ddrgev.f

来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 806 行 · 第 1/3 页

                  Q( N, N ) = ONE
                  WORK( N ) = ZERO
                  WORK( 3*N ) = SIGN( ONE, DLARND( 2, ISEED ) )
                  Z( N, N ) = ONE
                  WORK( 2*N ) = ZERO
                  WORK( 4*N ) = SIGN( ONE, DLARND( 2, ISEED ) )
*
*                 Apply the diagonal matrices
*
                  DO 60 JC = 1, N
                     DO 50 JR = 1, N
                        A( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )*
     $                                A( JR, JC )
                        B( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )*
     $                                B( JR, JC )
   50                CONTINUE
   60             CONTINUE
                  CALL DORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A,
     $                         LDA, WORK( 2*N+1 ), IERR )
                  IF( IERR.NE.0 )
     $               GO TO 90
                  CALL DORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ),
     $                         A, LDA, WORK( 2*N+1 ), IERR )
                  IF( IERR.NE.0 )
     $               GO TO 90
                  CALL DORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B,
     $                         LDA, WORK( 2*N+1 ), IERR )
                  IF( IERR.NE.0 )
     $               GO TO 90
                  CALL DORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ),
     $                         B, LDA, WORK( 2*N+1 ), IERR )
                  IF( IERR.NE.0 )
     $               GO TO 90
               END IF
            ELSE
*
*              Random matrices
*
               DO 80 JC = 1, N
                  DO 70 JR = 1, N
                     A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
     $                             DLARND( 2, ISEED )
                     B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
     $                             DLARND( 2, ISEED )
   70             CONTINUE
   80          CONTINUE
            END IF
*
   90       CONTINUE
*
            IF( IERR.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'Generator', IERR, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IERR )
               RETURN
            END IF
*
  100       CONTINUE
*
            DO 110 I = 1, 7
               RESULT( I ) = -ONE
  110       CONTINUE
*
*           Call DGGEV to compute eigenvalues and eigenvectors.
*
            CALL DLACPY( ' ', N, N, A, LDA, S, LDA )
            CALL DLACPY( ' ', N, N, B, LDA, T, LDA )
            CALL DGGEV( 'V', 'V', N, S, LDA, T, LDA, ALPHAR, ALPHAI,
     $                  BETA, Q, LDQ, Z, LDQ, WORK, LWORK, IERR )
            IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
               RESULT( 1 ) = ULPINV
               WRITE( NOUNIT, FMT = 9999 )'DGGEV1', IERR, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IERR )
               GO TO 190
            END IF
*
*           Do the tests (1) and (2)
*
            CALL DGET52( .TRUE., N, A, LDA, B, LDA, Q, LDQ, ALPHAR,
     $                   ALPHAI, BETA, WORK, RESULT( 1 ) )
            IF( RESULT( 2 ).GT.THRESH ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Left', 'DGGEV1',
     $            RESULT( 2 ), N, JTYPE, IOLDSD
            END IF
*
*           Do the tests (3) and (4)
*
            CALL DGET52( .FALSE., N, A, LDA, B, LDA, Z, LDQ, ALPHAR,
     $                   ALPHAI, BETA, WORK, RESULT( 3 ) )
            IF( RESULT( 4 ).GT.THRESH ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Right', 'DGGEV1',
     $            RESULT( 4 ), N, JTYPE, IOLDSD
            END IF
*
*           Do the test (5)
*
            CALL DLACPY( ' ', N, N, A, LDA, S, LDA )
            CALL DLACPY( ' ', N, N, B, LDA, T, LDA )
            CALL DGGEV( 'N', 'N', N, S, LDA, T, LDA, ALPHR1, ALPHI1,
     $                  BETA1, Q, LDQ, Z, LDQ, WORK, LWORK, IERR )
            IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
               RESULT( 1 ) = ULPINV
               WRITE( NOUNIT, FMT = 9999 )'DGGEV2', IERR, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IERR )
               GO TO 190
            END IF
*
            DO 120 J = 1, N
               IF( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE.
     $             ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) )RESULT( 5 )
     $              = ULPINV
  120       CONTINUE
*
*           Do the test (6): Compute eigenvalues and left eigenvectors,
*           and test them
*
            CALL DLACPY( ' ', N, N, A, LDA, S, LDA )
            CALL DLACPY( ' ', N, N, B, LDA, T, LDA )
            CALL DGGEV( 'V', 'N', N, S, LDA, T, LDA, ALPHR1, ALPHI1,
     $                  BETA1, QE, LDQE, Z, LDQ, WORK, LWORK, IERR )
            IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
               RESULT( 1 ) = ULPINV
               WRITE( NOUNIT, FMT = 9999 )'DGGEV3', IERR, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IERR )
               GO TO 190
            END IF
*
            DO 130 J = 1, N
               IF( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE.
     $             ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) )RESULT( 6 )
     $              = ULPINV
  130       CONTINUE
*
            DO 150 J = 1, N
               DO 140 JC = 1, N
                  IF( Q( J, JC ).NE.QE( J, JC ) )
     $               RESULT( 6 ) = ULPINV
  140          CONTINUE
  150       CONTINUE
*
*           DO the test (7): Compute eigenvalues and right eigenvectors,
*           and test them
*
            CALL DLACPY( ' ', N, N, A, LDA, S, LDA )
            CALL DLACPY( ' ', N, N, B, LDA, T, LDA )
            CALL DGGEV( 'N', 'V', N, S, LDA, T, LDA, ALPHR1, ALPHI1,
     $                  BETA1, Q, LDQ, QE, LDQE, WORK, LWORK, IERR )
            IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
               RESULT( 1 ) = ULPINV
               WRITE( NOUNIT, FMT = 9999 )'DGGEV4', IERR, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IERR )
               GO TO 190
            END IF
*
            DO 160 J = 1, N
               IF( ALPHAR( J ).NE.ALPHR1( J ) .OR. ALPHAI( J ).NE.
     $             ALPHI1( J ) .OR. BETA( J ).NE.BETA1( J ) )RESULT( 7 )
     $              = ULPINV
  160       CONTINUE
*
            DO 180 J = 1, N
               DO 170 JC = 1, N
                  IF( Z( J, JC ).NE.QE( J, JC ) )
     $               RESULT( 7 ) = ULPINV
  170          CONTINUE
  180       CONTINUE
*
*           End of Loop -- Check for RESULT(j) > THRESH
*
  190       CONTINUE
*
            NTESTT = NTESTT + 7
*
*           Print out tests which fail.
*
            DO 200 JR = 1, 7
               IF( RESULT( JR ).GE.THRESH ) THEN
*
*                 If this is the first test to fail,
*                 print a header to the data file.
*
                  IF( NERRS.EQ.0 ) THEN
                     WRITE( NOUNIT, FMT = 9997 )'DGV'
*
*                    Matrix types
*
                     WRITE( NOUNIT, FMT = 9996 )
                     WRITE( NOUNIT, FMT = 9995 )
                     WRITE( NOUNIT, FMT = 9994 )'Orthogonal'
*
*                    Tests performed
*
                     WRITE( NOUNIT, FMT = 9993 )
*
                  END IF
                  NERRS = NERRS + 1
                  IF( RESULT( JR ).LT.10000.0D0 ) THEN
                     WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR,
     $                  RESULT( JR )
                  ELSE
                     WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
     $                  RESULT( JR )
                  END IF
               END IF
  200       CONTINUE
*
  210    CONTINUE
  220 CONTINUE
*
*     Summary
*
      CALL ALASVM( 'DGV', NOUNIT, NERRS, NTESTT, 0 )
*
      WORK( 1 ) = MAXWRK
*
      RETURN
*
 9999 FORMAT( ' DDRGEV: ', A, ' returned INFO=', I6, '.', / 3X, 'N=',
     $      I6, ', JTYPE=', I6, ', ISEED=(', 4( I4, ',' ), I5, ')' )
*
 9998 FORMAT( ' DDRGEV: ', A, ' Eigenvectors from ', A, ' incorrectly ',
     $      'normalized.', / ' Bits of error=', 0P, G10.3, ',', 3X,
     $      'N=', I4, ', JTYPE=', I3, ', ISEED=(', 4( I4, ',' ), I5,
     $      ')' )
*
 9997 FORMAT( / 1X, A3, ' -- Real Generalized eigenvalue problem driver'
     $       )
*
 9996 FORMAT( ' Matrix types (see DDRGEV for details): ' )
*
 9995 FORMAT( ' Special Matrices:', 23X,
     $      '(J''=transposed Jordan block)',
     $      / '   1=(0,0)  2=(I,0)  3=(0,I)  4=(I,I)  5=(J'',J'')  ',
     $      '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices:  ( ',
     $      'D=diag(0,1,2,...) )', / '   7=(D,I)   9=(large*D, small*I',
     $      ')  11=(large*I, small*D)  13=(large*D, large*I)', /
     $      '   8=(I,D)  10=(small*D, large*I)  12=(small*I, large*D) ',
     $      ' 14=(small*D, small*I)', / '  15=(D, reversed D)' )
 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
     $      / '  16=Transposed Jordan Blocks             19=geometric ',
     $      'alpha, beta=0,1', / '  17=arithm. alpha&beta             ',
     $      '      20=arithmetic alpha, beta=0,1', / '  18=clustered ',
     $      'alpha, beta=0,1            21=random alpha, beta=0,1',
     $      / ' Large & Small Matrices:', / '  22=(large, small)   ',
     $      '23=(small,large)    24=(small,small)    25=(large,large)',
     $      / '  26=random O(1) matrices.' )
*
 9993 FORMAT( / ' Tests performed:    ',
     $      / ' 1 = max | ( b A - a B )''*l | / const.,',
     $      / ' 2 = | |VR(i)| - 1 | / ulp,',
     $      / ' 3 = max | ( b A - a B )*r | / const.',
     $      / ' 4 = | |VL(i)| - 1 | / ulp,',
     $      / ' 5 = 0 if W same no matter if r or l computed,',
     $      / ' 6 = 0 if l same no matter if l computed,',
     $      / ' 7 = 0 if r same no matter if r computed,', / 1X )
 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
     $      4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
     $      4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 )
*
*     End of DDRGEV
*
      END