cchkgg.f

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

*
            CALL CGET51( 1, N, A, LDA, H, LDA, U, LDU, V, LDU, WORK,
     $                   RWORK, RESULT( 1 ) )
            CALL CGET51( 1, N, B, LDA, T, LDA, U, LDU, V, LDU, WORK,
     $                   RWORK, RESULT( 2 ) )
            CALL CGET51( 3, N, B, LDA, T, LDA, U, LDU, U, LDU, WORK,
     $                   RWORK, RESULT( 3 ) )
            CALL CGET51( 3, N, B, LDA, T, LDA, V, LDU, V, LDU, WORK,
     $                   RWORK, RESULT( 4 ) )
*
*           Call CHGEQZ to compute S1, P1, S2, P2, Q, and Z, do tests.
*
*           Compute T1 and UZ
*
*           Eigenvalues only
*
            CALL CLACPY( ' ', N, N, H, LDA, S2, LDA )
            CALL CLACPY( ' ', N, N, T, LDA, P2, LDA )
            NTEST = 5
            RESULT( 5 ) = ULPINV
*
            CALL CHGEQZ( 'E', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
     $                   ALPHA3, BETA3, Q, LDU, Z, LDU, WORK, LWORK,
     $                   RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(E)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               GO TO 210
            END IF
*
*           Eigenvalues and Full Schur Form
*
            CALL CLACPY( ' ', N, N, H, LDA, S2, LDA )
            CALL CLACPY( ' ', N, N, T, LDA, P2, LDA )
*
            CALL CHGEQZ( 'S', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
     $                   ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
     $                   RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(S)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               GO TO 210
            END IF
*
*           Eigenvalues, Schur Form, and Schur Vectors
*
            CALL CLACPY( ' ', N, N, H, LDA, S1, LDA )
            CALL CLACPY( ' ', N, N, T, LDA, P1, LDA )
*
            CALL CHGEQZ( 'S', 'I', 'I', N, 1, N, S1, LDA, P1, LDA,
     $                   ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
     $                   RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(V)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               GO TO 210
            END IF
*
            NTEST = 8
*
*           Do Tests 5--8
*
            CALL CGET51( 1, N, H, LDA, S1, LDA, Q, LDU, Z, LDU, WORK,
     $                   RWORK, RESULT( 5 ) )
            CALL CGET51( 1, N, T, LDA, P1, LDA, Q, LDU, Z, LDU, WORK,
     $                   RWORK, RESULT( 6 ) )
            CALL CGET51( 3, N, T, LDA, P1, LDA, Q, LDU, Q, LDU, WORK,
     $                   RWORK, RESULT( 7 ) )
            CALL CGET51( 3, N, T, LDA, P1, LDA, Z, LDU, Z, LDU, WORK,
     $                   RWORK, RESULT( 8 ) )
*
*           Compute the Left and Right Eigenvectors of (S1,P1)
*
*           9: Compute the left eigenvector Matrix without
*              back transforming:
*
            NTEST = 9
            RESULT( 9 ) = ULPINV
*
*           To test "SELECT" option, compute half of the eigenvectors
*           in one call, and half in another
*
            I1 = N / 2
            DO 120 J = 1, I1
               LLWORK( J ) = .TRUE.
  120       CONTINUE
            DO 130 J = I1 + 1, N
               LLWORK( J ) = .FALSE.
  130       CONTINUE
*
            CALL CTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
     $                   LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,S1)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 210
            END IF
*
            I1 = IN
            DO 140 J = 1, I1
               LLWORK( J ) = .FALSE.
  140       CONTINUE
            DO 150 J = I1 + 1, N
               LLWORK( J ) = .TRUE.
  150       CONTINUE
*
            CALL CTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA,
     $                   EVECTL( 1, I1+1 ), LDU, CDUMMA, LDU, N, IN,
     $                   WORK, RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,S2)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 210
            END IF
*
            CALL CGET52( .TRUE., N, S1, LDA, P1, LDA, EVECTL, LDU,
     $                   ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
            RESULT( 9 ) = DUMMA( 1 )
            IF( DUMMA( 2 ).GT.THRSHN ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Left', 'CTGEVC(HOWMNY=S)',
     $            DUMMA( 2 ), N, JTYPE, IOLDSD
            END IF
*
*           10: Compute the left eigenvector Matrix with
*               back transforming:
*
            NTEST = 10
            RESULT( 10 ) = ULPINV
            CALL CLACPY( 'F', N, N, Q, LDU, EVECTL, LDU )
            CALL CTGEVC( 'L', 'B', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
     $                   LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,B)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 210
            END IF
*
            CALL CGET52( .TRUE., N, H, LDA, T, LDA, EVECTL, LDU, ALPHA1,
     $                   BETA1, WORK, RWORK, DUMMA( 1 ) )
            RESULT( 10 ) = DUMMA( 1 )
            IF( DUMMA( 2 ).GT.THRSHN ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Left', 'CTGEVC(HOWMNY=B)',
     $            DUMMA( 2 ), N, JTYPE, IOLDSD
            END IF
*
*           11: Compute the right eigenvector Matrix without
*               back transforming:
*
            NTEST = 11
            RESULT( 11 ) = ULPINV
*
*           To test "SELECT" option, compute half of the eigenvectors
*           in one call, and half in another
*
            I1 = N / 2
            DO 160 J = 1, I1
               LLWORK( J ) = .TRUE.
  160       CONTINUE
            DO 170 J = I1 + 1, N
               LLWORK( J ) = .FALSE.
  170       CONTINUE
*
            CALL CTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
     $                   LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,S1)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 210
            END IF
*
            I1 = IN
            DO 180 J = 1, I1
               LLWORK( J ) = .FALSE.
  180       CONTINUE
            DO 190 J = I1 + 1, N
               LLWORK( J ) = .TRUE.
  190       CONTINUE
*
            CALL CTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
     $                   LDU, EVECTR( 1, I1+1 ), LDU, N, IN, WORK,
     $                   RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,S2)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 210
            END IF
*
            CALL CGET52( .FALSE., N, S1, LDA, P1, LDA, EVECTR, LDU,
     $                   ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
            RESULT( 11 ) = DUMMA( 1 )
            IF( DUMMA( 2 ).GT.THRESH ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Right', 'CTGEVC(HOWMNY=S)',
     $            DUMMA( 2 ), N, JTYPE, IOLDSD
            END IF
*
*           12: Compute the right eigenvector Matrix with
*               back transforming:
*
            NTEST = 12
            RESULT( 12 ) = ULPINV
            CALL CLACPY( 'F', N, N, Z, LDU, EVECTR, LDU )
            CALL CTGEVC( 'R', 'B', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
     $                   LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,B)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 210
            END IF
*
            CALL CGET52( .FALSE., N, H, LDA, T, LDA, EVECTR, LDU,
     $                   ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
            RESULT( 12 ) = DUMMA( 1 )
            IF( DUMMA( 2 ).GT.THRESH ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Right', 'CTGEVC(HOWMNY=B)',
     $            DUMMA( 2 ), N, JTYPE, IOLDSD
            END IF
*
*           Tests 13--15 are done only on request
*
            IF( TSTDIF ) THEN
*
*              Do Tests 13--14
*
               CALL CGET51( 2, N, S1, LDA, S2, LDA, Q, LDU, Z, LDU,
     $                      WORK, RWORK, RESULT( 13 ) )
               CALL CGET51( 2, N, P1, LDA, P2, LDA, Q, LDU, Z, LDU,
     $                      WORK, RWORK, RESULT( 14 ) )
*
*              Do Test 15
*
               TEMP1 = ZERO
               TEMP2 = ZERO
               DO 200 J = 1, N
                  TEMP1 = MAX( TEMP1, ABS( ALPHA1( J )-ALPHA3( J ) ) )
                  TEMP2 = MAX( TEMP2, ABS( BETA1( J )-BETA3( J ) ) )
  200          CONTINUE
*
               TEMP1 = TEMP1 / MAX( SAFMIN, ULP*MAX( TEMP1, ANORM ) )
               TEMP2 = TEMP2 / MAX( SAFMIN, ULP*MAX( TEMP2, BNORM ) )
               RESULT( 15 ) = MAX( TEMP1, TEMP2 )
               NTEST = 15
            ELSE
               RESULT( 13 ) = ZERO
               RESULT( 14 ) = ZERO
               RESULT( 15 ) = ZERO
               NTEST = 12
            END IF
*
*           End of Loop -- Check for RESULT(j) > THRESH
*
  210       CONTINUE
*
            NTESTT = NTESTT + NTEST
*
*           Print out tests which fail.
*
            DO 220 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 )'CGG'
*
*                    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, 10 )
*
                  END IF
                  NERRS = NERRS + 1
                  IF( RESULT( JR ).LT.10000.0 ) 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
  220       CONTINUE
*
  230    CONTINUE
  240 CONTINUE
*
*     Summary
*
      CALL SLASUM( 'CGG', NOUNIT, NERRS, NTESTT )
      RETURN
*
 9999 FORMAT( ' CCHKGG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
*
 9998 FORMAT( ' CCHKGG: ', 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' )
*
 9996 FORMAT( ' Matrix types (see CCHKGG 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:   (H is Hessenberg, S is Schur, B, ',
     $      'T, P are triangular,', / 20X, 'U, V, Q, and Z are ', A,
     $      ', l and r are the', / 20X,
     $      'appropriate left and right eigenvectors, resp., a is',
     $      / 20X, 'alpha, b is beta, and ', A, ' means ', A, '.)',
     $      / ' 1 = | A - U H V', A,
     $      ' | / ( |A| n ulp )      2 = | B - U T V', A,
     $      ' | / ( |B| n ulp )', / ' 3 = | I - UU', A,
     $      ' | / ( n ulp )             4 = | I - VV', A,
     $      ' | / ( n ulp )', / ' 5 = | H - Q S Z', A,
     $      ' | / ( |H| n ulp )', 6X, '6 = | T - Q P Z', A,
     $      ' | / ( |T| n ulp )', / ' 7 = | I - QQ', A,
     $      ' | / ( n ulp )             8 = | I - ZZ', A,
     $      ' | / ( n ulp )', / ' 9 = max | ( b S - a P )', A,
     $      ' l | / const.  10 = max | ( b H - a T )', A,
     $      ' l | / const.', /
     $      ' 11= max | ( b S - a P ) r | / const.   12 = max | ( b H',
     $      ' - a T ) r | / const.', / 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, E10.3 )
*
*     End of CCHKGG
*
      END