zdrvgg.f

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

     $                      ISEED, B, LDA )
               IADD = KADD( KBZERO( JTYPE ) )
               IF( IADD.NE.0 .AND. IADD.LE.N )
     $            B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) )
*
               IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
*
*                 Include rotations
*
*                 Generate Q, Z as Householder transformations times
*                 a diagonal matrix.
*
                  DO 50 JC = 1, N - 1
                     DO 40 JR = JC, N
                        Q( JR, JC ) = ZLARND( 3, ISEED )
                        Z( JR, JC ) = ZLARND( 3, ISEED )
   40                CONTINUE
                     CALL ZLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1,
     $                            WORK( JC ) )
                     WORK( 2*N+JC ) = SIGN( ONE, DBLE( Q( JC, JC ) ) )
                     Q( JC, JC ) = CONE
                     CALL ZLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1,
     $                            WORK( N+JC ) )
                     WORK( 3*N+JC ) = SIGN( ONE, DBLE( Z( JC, JC ) ) )
                     Z( JC, JC ) = CONE
   50             CONTINUE
                  CTEMP = ZLARND( 3, ISEED )
                  Q( N, N ) = CONE
                  WORK( N ) = CZERO
                  WORK( 3*N ) = CTEMP / ABS( CTEMP )
                  CTEMP = ZLARND( 3, ISEED )
                  Z( N, N ) = CONE
                  WORK( 2*N ) = CZERO
                  WORK( 4*N ) = CTEMP / ABS( CTEMP )
*
*                 Apply the diagonal matrices
*
                  DO 70 JC = 1, N
                     DO 60 JR = 1, N
                        A( JR, JC ) = WORK( 2*N+JR )*
     $                                DCONJG( WORK( 3*N+JC ) )*
     $                                A( JR, JC )
                        B( JR, JC ) = WORK( 2*N+JR )*
     $                                DCONJG( WORK( 3*N+JC ) )*
     $                                B( JR, JC )
   60                CONTINUE
   70             CONTINUE
                  CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A,
     $                         LDA, WORK( 2*N+1 ), IINFO )
                  IF( IINFO.NE.0 )
     $               GO TO 100
                  CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ),
     $                         A, LDA, WORK( 2*N+1 ), IINFO )
                  IF( IINFO.NE.0 )
     $               GO TO 100
                  CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B,
     $                         LDA, WORK( 2*N+1 ), IINFO )
                  IF( IINFO.NE.0 )
     $               GO TO 100
                  CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ),
     $                         B, LDA, WORK( 2*N+1 ), IINFO )
                  IF( IINFO.NE.0 )
     $               GO TO 100
               END IF
            ELSE
*
*              Random matrices
*
               DO 90 JC = 1, N
                  DO 80 JR = 1, N
                     A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
     $                             ZLARND( 4, ISEED )
                     B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
     $                             ZLARND( 4, ISEED )
   80             CONTINUE
   90          CONTINUE
            END IF
*
  100       CONTINUE
*
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               RETURN
            END IF
*
  110       CONTINUE
*
*           Call ZGEGS to compute H, T, Q, Z, alpha, and beta.
*
            CALL ZLACPY( ' ', N, N, A, LDA, S, LDA )
            CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
            NTEST = 1
            RESULT( 1 ) = ULPINV
*
            CALL ZGEGS( 'V', 'V', N, S, LDA, T, LDA, ALPHA1, BETA1, Q,
     $                  LDQ, Z, LDQ, WORK, LWORK, RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZGEGS', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               GO TO 130
            END IF
*
            NTEST = 4
*
*           Do tests 1--4
*
            CALL ZGET51( 1, N, A, LDA, S, LDA, Q, LDQ, Z, LDQ, WORK,
     $                   RWORK, RESULT( 1 ) )
            CALL ZGET51( 1, N, B, LDA, T, LDA, Q, LDQ, Z, LDQ, WORK,
     $                   RWORK, RESULT( 2 ) )
            CALL ZGET51( 3, N, B, LDA, T, LDA, Q, LDQ, Q, LDQ, WORK,
     $                   RWORK, RESULT( 3 ) )
            CALL ZGET51( 3, N, B, LDA, T, LDA, Z, LDQ, Z, LDQ, WORK,
     $                   RWORK, RESULT( 4 ) )
*
*           Do test 5: compare eigenvalues with diagonals.
*
            TEMP1 = ZERO
*
            DO 120 J = 1, N
               TEMP2 = ( ABS1( ALPHA1( J )-S( J, J ) ) /
     $                 MAX( SAFMIN, ABS1( ALPHA1( J ) ), ABS1( S( J,
     $                 J ) ) )+ABS1( BETA1( J )-T( J, J ) ) /
     $                 MAX( SAFMIN, ABS1( BETA1( J ) ), ABS1( T( J,
     $                 J ) ) ) ) / ULP
               TEMP1 = MAX( TEMP1, TEMP2 )
  120       CONTINUE
            RESULT( 5 ) = TEMP1
*
*           Call ZGEGV to compute S2, T2, VL, and VR, do tests.
*
*           Eigenvalues and Eigenvectors
*
            CALL ZLACPY( ' ', N, N, A, LDA, S2, LDA )
            CALL ZLACPY( ' ', N, N, B, LDA, T2, LDA )
            NTEST = 6
            RESULT( 6 ) = ULPINV
*
            CALL ZGEGV( 'V', 'V', N, S2, LDA, T2, LDA, ALPHA2, BETA2,
     $                  VL, LDQ, VR, LDQ, WORK, LWORK, RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZGEGV', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               GO TO 130
            END IF
*
            NTEST = 7
*
*           Do Tests 6 and 7
*
            CALL ZGET52( .TRUE., N, A, LDA, B, LDA, VL, LDQ, ALPHA2,
     $                   BETA2, WORK, RWORK, DUMMA( 1 ) )
            RESULT( 6 ) = DUMMA( 1 )
            IF( DUMMA( 2 ).GT.THRSHN ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Left', 'ZGEGV', DUMMA( 2 ),
     $            N, JTYPE, IOLDSD
            END IF
*
            CALL ZGET52( .FALSE., N, A, LDA, B, LDA, VR, LDQ, ALPHA2,
     $                   BETA2, WORK, RWORK, DUMMA( 1 ) )
            RESULT( 7 ) = DUMMA( 1 )
            IF( DUMMA( 2 ).GT.THRESH ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Right', 'ZGEGV', DUMMA( 2 ),
     $            N, JTYPE, IOLDSD
            END IF
*
*           End of Loop -- Check for RESULT(j) > THRESH
*
  130       CONTINUE
*
            NTESTT = NTESTT + NTEST
*
*           Print out tests which fail.
*
            DO 140 JR = 1, NTEST
               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 )'ZGG'
*
*                    Matrix types
*
                     WRITE( NOUNIT, FMT = 9996 )
                     WRITE( NOUNIT, FMT = 9995 )
                     WRITE( NOUNIT, FMT = 9994 )'Unitary'
*
*                    Tests performed
*
                     WRITE( NOUNIT, FMT = 9993 )'unitary', '*',
     $                  'conjugate transpose', ( '*', J = 1, 5 )
*
                  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
  140       CONTINUE
*
  150    CONTINUE
  160 CONTINUE
*
*     Summary
*
      CALL ALASVM( 'ZGG', NOUNIT, NERRS, NTESTT, 0 )
      RETURN
*
 9999 FORMAT( ' ZDRVGG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
*
 9998 FORMAT( ' ZDRVGG: ', A, ' Eigenvectors from ', A, ' incorrectly ',
     $      'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
     $      'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5,
     $      ')' )
*
 9997 FORMAT( / 1X, A3,
     $      ' -- Complex Generalized eigenvalue problem driver' )
*
 9996 FORMAT( ' Matrix types (see ZDRVGG 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:  (S is Schur, T is triangular, ',
     $      'Q and Z are ', A, ',', / 20X,
     $      'l and r are the appropriate left and right', / 19X,
     $      'eigenvectors, resp., a is alpha, b is beta, and', / 19X, A,
     $      ' means ', A, '.)', / ' 1 = | A - Q S Z', A,
     $      ' | / ( |A| n ulp )      2 = | B - Q T Z', A,
     $      ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A,
     $      ' | / ( n ulp )             4 = | I - ZZ', A,
     $      ' | / ( n ulp )', /
     $      ' 5 = difference between (alpha,beta) and diagonals of',
     $      ' (S,T)', / ' 6 = max | ( b A - a B )', A,
     $      ' l | / const.   7 = max | ( b A - a B ) r | / const.',
     $      / 1X )
 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
     $      4( I4, ',' ), ' result ', I3, ' is', 0P, F8.2 )
 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
     $      4( I4, ',' ), ' result ', I3, ' is', 1P, D10.3 )
*
*     End of ZDRVGG
*
      END