zchkhs.f

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

            NTEST = 2
*
            CALL ZHST01( N, ILO, IHI, A, LDA, H, LDA, U, LDU, WORK,
     $                   NWORK, RWORK, RESULT( 1 ) )
*
*           Call ZHSEQR to compute T1, T2 and Z, do tests.
*
*           Eigenvalues only (W3)
*
            CALL ZLACPY( ' ', N, N, H, LDA, T2, LDA )
            NTEST = 3
            RESULT( 3 ) = ULPINV
*
            CALL ZHSEQR( 'E', 'N', N, ILO, IHI, T2, LDA, W3, UZ, LDU,
     $                   WORK, NWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(E)', IINFO, N, JTYPE,
     $            IOLDSD
               IF( IINFO.LE.N+2 ) THEN
                  INFO = ABS( IINFO )
                  GO TO 240
               END IF
            END IF
*
*           Eigenvalues (W1) and Full Schur Form (T2)
*
            CALL ZLACPY( ' ', N, N, H, LDA, T2, LDA )
*
            CALL ZHSEQR( 'S', 'N', N, ILO, IHI, T2, LDA, W1, UZ, LDU,
     $                   WORK, NWORK, IINFO )
            IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(S)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               GO TO 240
            END IF
*
*           Eigenvalues (W1), Schur Form (T1), and Schur Vectors (UZ)
*
            CALL ZLACPY( ' ', N, N, H, LDA, T1, LDA )
            CALL ZLACPY( ' ', N, N, U, LDU, UZ, LDU )
*
            CALL ZHSEQR( 'S', 'V', N, ILO, IHI, T1, LDA, W1, UZ, LDU,
     $                   WORK, NWORK, IINFO )
            IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZHSEQR(V)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               GO TO 240
            END IF
*
*           Compute Z = U' UZ
*
            CALL ZGEMM( 'C', 'N', N, N, N, CONE, U, LDU, UZ, LDU, CZERO,
     $                  Z, LDU )
            NTEST = 8
*
*           Do Tests 3: | H - Z T Z' | / ( |H| n ulp )
*                and 4: | I - Z Z' | / ( n ulp )
*
            CALL ZHST01( N, ILO, IHI, H, LDA, T1, LDA, Z, LDU, WORK,
     $                   NWORK, RWORK, RESULT( 3 ) )
*
*           Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp )
*                and 6: | I - UZ (UZ)' | / ( n ulp )
*
            CALL ZHST01( N, ILO, IHI, A, LDA, T1, LDA, UZ, LDU, WORK,
     $                   NWORK, RWORK, RESULT( 5 ) )
*
*           Do Test 7: | T2 - T1 | / ( |T| n ulp )
*
            CALL ZGET10( N, N, T2, LDA, T1, LDA, WORK, RWORK,
     $                   RESULT( 7 ) )
*
*           Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp )
*
            TEMP1 = ZERO
            TEMP2 = ZERO
            DO 130 J = 1, N
               TEMP1 = MAX( TEMP1, ABS( W1( J ) ), ABS( W3( J ) ) )
               TEMP2 = MAX( TEMP2, ABS( W1( J )-W3( J ) ) )
  130       CONTINUE
*
            RESULT( 8 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
*
*           Compute the Left and Right Eigenvectors of T
*
*           Compute the Right eigenvector Matrix:
*
            NTEST = 9
            RESULT( 9 ) = ULPINV
*
*           Select every other eigenvector
*
            DO 140 J = 1, N
               SELECT( J ) = .FALSE.
  140       CONTINUE
            DO 150 J = 1, N, 2
               SELECT( J ) = .TRUE.
  150       CONTINUE
            CALL ZTREVC( 'Right', 'All', SELECT, N, T1, LDA, CDUMMA,
     $                   LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZTREVC(R,A)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 240
            END IF
*
*           Test 9:  | TR - RW | / ( |T| |R| ulp )
*
            CALL ZGET22( 'N', 'N', 'N', N, T1, LDA, EVECTR, LDU, W1,
     $                   WORK, RWORK, DUMMA( 1 ) )
            RESULT( 9 ) = DUMMA( 1 )
            IF( DUMMA( 2 ).GT.THRESH ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Right', 'ZTREVC',
     $            DUMMA( 2 ), N, JTYPE, IOLDSD
            END IF
*
*           Compute selected right eigenvectors and confirm that
*           they agree with previous right eigenvectors
*
            CALL ZTREVC( 'Right', 'Some', SELECT, N, T1, LDA, CDUMMA,
     $                   LDU, EVECTL, LDU, N, IN, WORK, RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZTREVC(R,S)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 240
            END IF
*
            K = 1
            MATCH = .TRUE.
            DO 170 J = 1, N
               IF( SELECT( J ) ) THEN
                  DO 160 JJ = 1, N
                     IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) ) THEN
                        MATCH = .FALSE.
                        GO TO 180
                     END IF
  160             CONTINUE
                  K = K + 1
               END IF
  170       CONTINUE
  180       CONTINUE
            IF( .NOT.MATCH )
     $         WRITE( NOUNIT, FMT = 9997 )'Right', 'ZTREVC', N, JTYPE,
     $         IOLDSD
*
*           Compute the Left eigenvector Matrix:
*
            NTEST = 10
            RESULT( 10 ) = ULPINV
            CALL ZTREVC( 'Left', 'All', SELECT, N, T1, LDA, EVECTL, LDU,
     $                   CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZTREVC(L,A)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 240
            END IF
*
*           Test 10:  | LT - WL | / ( |T| |L| ulp )
*
            CALL ZGET22( 'C', 'N', 'C', N, T1, LDA, EVECTL, LDU, W1,
     $                   WORK, RWORK, DUMMA( 3 ) )
            RESULT( 10 ) = DUMMA( 3 )
            IF( DUMMA( 4 ).GT.THRESH ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Left', 'ZTREVC', DUMMA( 4 ),
     $            N, JTYPE, IOLDSD
            END IF
*
*           Compute selected left eigenvectors and confirm that
*           they agree with previous left eigenvectors
*
            CALL ZTREVC( 'Left', 'Some', SELECT, N, T1, LDA, EVECTR,
     $                   LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZTREVC(L,S)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 240
            END IF
*
            K = 1
            MATCH = .TRUE.
            DO 200 J = 1, N
               IF( SELECT( J ) ) THEN
                  DO 190 JJ = 1, N
                     IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) ) THEN
                        MATCH = .FALSE.
                        GO TO 210
                     END IF
  190             CONTINUE
                  K = K + 1
               END IF
  200       CONTINUE
  210       CONTINUE
            IF( .NOT.MATCH )
     $         WRITE( NOUNIT, FMT = 9997 )'Left', 'ZTREVC', N, JTYPE,
     $         IOLDSD
*
*           Call ZHSEIN for Right eigenvectors of H, do test 11
*
            NTEST = 11
            RESULT( 11 ) = ULPINV
            DO 220 J = 1, N
               SELECT( J ) = .TRUE.
  220       CONTINUE
*
            CALL ZHSEIN( 'Right', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3,
     $                   CDUMMA, LDU, EVECTX, LDU, N1, IN, WORK, RWORK,
     $                   IWORK, IWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZHSEIN(R)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 )
     $            GO TO 240
            ELSE
*
*              Test 11:  | HX - XW | / ( |H| |X| ulp )
*
*                        (from inverse iteration)
*
               CALL ZGET22( 'N', 'N', 'N', N, H, LDA, EVECTX, LDU, W3,
     $                      WORK, RWORK, DUMMA( 1 ) )
               IF( DUMMA( 1 ).LT.ULPINV )
     $            RESULT( 11 ) = DUMMA( 1 )*ANINV
               IF( DUMMA( 2 ).GT.THRESH ) THEN
                  WRITE( NOUNIT, FMT = 9998 )'Right', 'ZHSEIN',
     $               DUMMA( 2 ), N, JTYPE, IOLDSD
               END IF
            END IF
*
*           Call ZHSEIN for Left eigenvectors of H, do test 12
*
            NTEST = 12
            RESULT( 12 ) = ULPINV
            DO 230 J = 1, N
               SELECT( J ) = .TRUE.
  230       CONTINUE
*
            CALL ZHSEIN( 'Left', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3,
     $                   EVECTY, LDU, CDUMMA, LDU, N1, IN, WORK, RWORK,
     $                   IWORK, IWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZHSEIN(L)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 )
     $            GO TO 240
            ELSE
*
*              Test 12:  | YH - WY | / ( |H| |Y| ulp )
*
*                        (from inverse iteration)
*
               CALL ZGET22( 'C', 'N', 'C', N, H, LDA, EVECTY, LDU, W3,
     $                      WORK, RWORK, DUMMA( 3 ) )
               IF( DUMMA( 3 ).LT.ULPINV )
     $            RESULT( 12 ) = DUMMA( 3 )*ANINV
               IF( DUMMA( 4 ).GT.THRESH ) THEN
                  WRITE( NOUNIT, FMT = 9998 )'Left', 'ZHSEIN',
     $               DUMMA( 4 ), N, JTYPE, IOLDSD
               END IF
            END IF
*
*           Call ZUNMHR for Right eigenvectors of A, do test 13
*
            NTEST = 13
            RESULT( 13 ) = ULPINV
*
            CALL ZUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU,
     $                   LDU, TAU, EVECTX, LDU, WORK, NWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZUNMHR(L)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 )
     $            GO TO 240
            ELSE
*
*              Test 13:  | AX - XW | / ( |A| |X| ulp )
*
*                        (from inverse iteration)
*
               CALL ZGET22( 'N', 'N', 'N', N, A, LDA, EVECTX, LDU, W3,
     $                      WORK, RWORK, DUMMA( 1 ) )
               IF( DUMMA( 1 ).LT.ULPINV )
     $            RESULT( 13 ) = DUMMA( 1 )*ANINV
            END IF
*
*           Call ZUNMHR for Left eigenvectors of A, do test 14
*
            NTEST = 14
            RESULT( 14 ) = ULPINV
*
            CALL ZUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU,
     $                   LDU, TAU, EVECTY, LDU, WORK, NWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'ZUNMHR(L)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 )
     $            GO TO 240
            ELSE
*
*              Test 14:  | YA - WY | / ( |A| |Y| ulp )
*
*                        (from inverse iteration)
*
               CALL ZGET22( 'C', 'N', 'C', N, A, LDA, EVECTY, LDU, W3,
     $                      WORK, RWORK, DUMMA( 3 ) )
               IF( DUMMA( 3 ).LT.ULPINV )
     $            RESULT( 14 ) = DUMMA( 3 )*ANINV
            END IF
*
*           End of Loop -- Check for RESULT(j) > THRESH
*
  240       CONTINUE
*
            NTESTT = NTESTT + NTEST
            CALL DLAFTS( 'ZHS', N, N, JTYPE, NTEST, RESULT, IOLDSD,
     $                   THRESH, NOUNIT, NERRS )
*
  250    CONTINUE
  260 CONTINUE
*
*     Summary
*
      CALL DLASUM( 'ZHS', NOUNIT, NERRS, NTESTT )
*
      RETURN
*
 9999 FORMAT( ' ZCHKHS: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
 9998 FORMAT( ' ZCHKHS: ', A, ' Eigenvectors from ', A, ' incorrectly ',
     $      'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
     $      'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5,
     $      ')' )
 9997 FORMAT( ' ZCHKHS: Selected ', A, ' Eigenvectors from ', A,
     $      ' do not match other eigenvectors ', 9X, 'N=', I6,
     $      ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
*
*     End of ZCHKHS
*
      END