dchkhs.f

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

            CALL DHSEQR( 'S', 'N', N, ILO, IHI, T2, LDA, WR1, WI1, UZ,
     $                   LDU, WORK, NWORK, IINFO )
            IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'DHSEQR(S)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               GO TO 250
            END IF
*
*           Eigenvalues (WR1,WI1), Schur Form (T1), and Schur vectors
*           (UZ)
*
            CALL DLACPY( ' ', N, N, H, LDA, T1, LDA )
            CALL DLACPY( ' ', N, N, U, LDU, UZ, LDA )
*
            CALL DHSEQR( 'S', 'V', N, ILO, IHI, T1, LDA, WR1, WI1, UZ,
     $                   LDU, WORK, NWORK, IINFO )
            IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'DHSEQR(V)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               GO TO 250
            END IF
*
*           Compute Z = U' UZ
*
            CALL DGEMM( 'T', 'N', N, N, N, ONE, U, LDU, UZ, LDU, ZERO,
     $                  Z, LDU )
            NTEST = 8
*
*           Do Tests 3: | H - Z T Z' | / ( |H| n ulp )
*                and 4: | I - Z Z' | / ( n ulp )
*
            CALL DHST01( N, ILO, IHI, H, LDA, T1, LDA, Z, LDU, WORK,
     $                   NWORK, RESULT( 3 ) )
*
*           Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp )
*                and 6: | I - UZ (UZ)' | / ( n ulp )
*
            CALL DHST01( N, ILO, IHI, A, LDA, T1, LDA, UZ, LDU, WORK,
     $                   NWORK, RESULT( 5 ) )
*
*           Do Test 7: | T2 - T1 | / ( |T| n ulp )
*
            CALL DGET10( N, N, T2, LDA, T1, LDA, WORK, RESULT( 7 ) )
*
*           Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp )
*
            TEMP1 = ZERO
            TEMP2 = ZERO
            DO 130 J = 1, N
               TEMP1 = MAX( TEMP1, ABS( WR1( J ) )+ABS( WI1( J ) ),
     $                 ABS( WR3( J ) )+ABS( WI3( J ) ) )
               TEMP2 = MAX( TEMP2, ABS( WR1( J )-WR3( J ) )+
     $                 ABS( WR1( J )-WR3( 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 last max(N/4,1) real, max(N/4,1) complex eigenvectors
*
            NSELC = 0
            NSELR = 0
            J = N
  140       CONTINUE
            IF( WI1( J ).EQ.ZERO ) THEN
               IF( NSELR.LT.MAX( N / 4, 1 ) ) THEN
                  NSELR = NSELR + 1
                  SELECT( J ) = .TRUE.
               ELSE
                  SELECT( J ) = .FALSE.
               END IF
               J = J - 1
            ELSE
               IF( NSELC.LT.MAX( N / 4, 1 ) ) THEN
                  NSELC = NSELC + 1
                  SELECT( J ) = .TRUE.
                  SELECT( J-1 ) = .FALSE.
               ELSE
                  SELECT( J ) = .FALSE.
                  SELECT( J-1 ) = .FALSE.
               END IF
               J = J - 2
            END IF
            IF( J.GT.0 )
     $         GO TO 140
*
            CALL DTREVC( 'Right', 'All', SELECT, N, T1, LDA, DUMMA, LDU,
     $                   EVECTR, LDU, N, IN, WORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'DTREVC(R,A)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 250
            END IF
*
*           Test 9:  | TR - RW | / ( |T| |R| ulp )
*
            CALL DGET22( 'N', 'N', 'N', N, T1, LDA, EVECTR, LDU, WR1,
     $                   WI1, WORK, DUMMA( 1 ) )
            RESULT( 9 ) = DUMMA( 1 )
            IF( DUMMA( 2 ).GT.THRESH ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Right', 'DTREVC',
     $            DUMMA( 2 ), N, JTYPE, IOLDSD
            END IF
*
*           Compute selected right eigenvectors and confirm that
*           they agree with previous right eigenvectors
*
            CALL DTREVC( 'Right', 'Some', SELECT, N, T1, LDA, DUMMA,
     $                   LDU, EVECTL, LDU, N, IN, WORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'DTREVC(R,S)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 250
            END IF
*
            K = 1
            MATCH = .TRUE.
            DO 170 J = 1, N
               IF( SELECT( J ) .AND. WI1( J ).EQ.ZERO ) THEN
                  DO 150 JJ = 1, N
                     IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) ) THEN
                        MATCH = .FALSE.
                        GO TO 180
                     END IF
  150             CONTINUE
                  K = K + 1
               ELSE IF( SELECT( J ) .AND. WI1( J ).NE.ZERO ) THEN
                  DO 160 JJ = 1, N
                     IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) .OR.
     $                   EVECTR( JJ, J+1 ).NE.EVECTL( JJ, K+1 ) ) THEN
                        MATCH = .FALSE.
                        GO TO 180
                     END IF
  160             CONTINUE
                  K = K + 2
               END IF
  170       CONTINUE
  180       CONTINUE
            IF( .NOT.MATCH )
     $         WRITE( NOUNIT, FMT = 9997 )'Right', 'DTREVC', N, JTYPE,
     $         IOLDSD
*
*           Compute the Left eigenvector Matrix:
*
            NTEST = 10
            RESULT( 10 ) = ULPINV
            CALL DTREVC( 'Left', 'All', SELECT, N, T1, LDA, EVECTL, LDU,
     $                   DUMMA, LDU, N, IN, WORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'DTREVC(L,A)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 250
            END IF
*
*           Test 10:  | LT - WL | / ( |T| |L| ulp )
*
            CALL DGET22( 'Trans', 'N', 'Conj', N, T1, LDA, EVECTL, LDU,
     $                   WR1, WI1, WORK, DUMMA( 3 ) )
            RESULT( 10 ) = DUMMA( 3 )
            IF( DUMMA( 4 ).GT.THRESH ) THEN
               WRITE( NOUNIT, FMT = 9998 )'Left', 'DTREVC', DUMMA( 4 ),
     $            N, JTYPE, IOLDSD
            END IF
*
*           Compute selected left eigenvectors and confirm that
*           they agree with previous left eigenvectors
*
            CALL DTREVC( 'Left', 'Some', SELECT, N, T1, LDA, EVECTR,
     $                   LDU, DUMMA, LDU, N, IN, WORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'DTREVC(L,S)', IINFO, N,
     $            JTYPE, IOLDSD
               INFO = ABS( IINFO )
               GO TO 250
            END IF
*
            K = 1
            MATCH = .TRUE.
            DO 210 J = 1, N
               IF( SELECT( J ) .AND. WI1( J ).EQ.ZERO ) THEN
                  DO 190 JJ = 1, N
                     IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) ) THEN
                        MATCH = .FALSE.
                        GO TO 220
                     END IF
  190             CONTINUE
                  K = K + 1
               ELSE IF( SELECT( J ) .AND. WI1( J ).NE.ZERO ) THEN
                  DO 200 JJ = 1, N
                     IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) .OR.
     $                   EVECTL( JJ, J+1 ).NE.EVECTR( JJ, K+1 ) ) THEN
                        MATCH = .FALSE.
                        GO TO 220
                     END IF
  200             CONTINUE
                  K = K + 2
               END IF
  210       CONTINUE
  220       CONTINUE
            IF( .NOT.MATCH )
     $         WRITE( NOUNIT, FMT = 9997 )'Left', 'DTREVC', N, JTYPE,
     $         IOLDSD
*
*           Call DHSEIN for Right eigenvectors of H, do test 11
*
            NTEST = 11
            RESULT( 11 ) = ULPINV
            DO 230 J = 1, N
               SELECT( J ) = .TRUE.
  230       CONTINUE
*
            CALL DHSEIN( 'Right', 'Qr', 'Ninitv', SELECT, N, H, LDA,
     $                   WR3, WI3, DUMMA, LDU, EVECTX, LDU, N1, IN,
     $                   WORK, IWORK, IWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'DHSEIN(R)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 )
     $            GO TO 250
            ELSE
*
*              Test 11:  | HX - XW | / ( |H| |X| ulp )
*
*                        (from inverse iteration)
*
               CALL DGET22( 'N', 'N', 'N', N, H, LDA, EVECTX, LDU, WR3,
     $                      WI3, WORK, DUMMA( 1 ) )
               IF( DUMMA( 1 ).LT.ULPINV )
     $            RESULT( 11 ) = DUMMA( 1 )*ANINV
               IF( DUMMA( 2 ).GT.THRESH ) THEN
                  WRITE( NOUNIT, FMT = 9998 )'Right', 'DHSEIN',
     $               DUMMA( 2 ), N, JTYPE, IOLDSD
               END IF
            END IF
*
*           Call DHSEIN for Left eigenvectors of H, do test 12
*
            NTEST = 12
            RESULT( 12 ) = ULPINV
            DO 240 J = 1, N
               SELECT( J ) = .TRUE.
  240       CONTINUE
*
            CALL DHSEIN( 'Left', 'Qr', 'Ninitv', SELECT, N, H, LDA, WR3,
     $                   WI3, EVECTY, LDU, DUMMA, LDU, N1, IN, WORK,
     $                   IWORK, IWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'DHSEIN(L)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 )
     $            GO TO 250
            ELSE
*
*              Test 12:  | YH - WY | / ( |H| |Y| ulp )
*
*                        (from inverse iteration)
*
               CALL DGET22( 'C', 'N', 'C', N, H, LDA, EVECTY, LDU, WR3,
     $                      WI3, WORK, DUMMA( 3 ) )
               IF( DUMMA( 3 ).LT.ULPINV )
     $            RESULT( 12 ) = DUMMA( 3 )*ANINV
               IF( DUMMA( 4 ).GT.THRESH ) THEN
                  WRITE( NOUNIT, FMT = 9998 )'Left', 'DHSEIN',
     $               DUMMA( 4 ), N, JTYPE, IOLDSD
               END IF
            END IF
*
*           Call DORMHR for Right eigenvectors of A, do test 13
*
            NTEST = 13
            RESULT( 13 ) = ULPINV
*
            CALL DORMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU,
     $                   LDU, TAU, EVECTX, LDU, WORK, NWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'DORMHR(R)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 )
     $            GO TO 250
            ELSE
*
*              Test 13:  | AX - XW | / ( |A| |X| ulp )
*
*                        (from inverse iteration)
*
               CALL DGET22( 'N', 'N', 'N', N, A, LDA, EVECTX, LDU, WR3,
     $                      WI3, WORK, DUMMA( 1 ) )
               IF( DUMMA( 1 ).LT.ULPINV )
     $            RESULT( 13 ) = DUMMA( 1 )*ANINV
            END IF
*
*           Call DORMHR for Left eigenvectors of A, do test 14
*
            NTEST = 14
            RESULT( 14 ) = ULPINV
*
            CALL DORMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU,
     $                   LDU, TAU, EVECTY, LDU, WORK, NWORK, IINFO )
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'DORMHR(L)', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               IF( IINFO.LT.0 )
     $            GO TO 250
            ELSE
*
*              Test 14:  | YA - WY | / ( |A| |Y| ulp )
*
*                        (from inverse iteration)
*
               CALL DGET22( 'C', 'N', 'C', N, A, LDA, EVECTY, LDU, WR3,
     $                      WI3, WORK, DUMMA( 3 ) )
               IF( DUMMA( 3 ).LT.ULPINV )
     $            RESULT( 14 ) = DUMMA( 3 )*ANINV
            END IF
*
*           End of Loop -- Check for RESULT(j) > THRESH
*
  250       CONTINUE
*
            NTESTT = NTESTT + NTEST
            CALL DLAFTS( 'DHS', N, N, JTYPE, NTEST, RESULT, IOLDSD,
     $                   THRESH, NOUNIT, NERRS )
*
  260    CONTINUE
  270 CONTINUE
*
*     Summary
*
      CALL DLASUM( 'DHS', NOUNIT, NERRS, NTESTT )
*
      RETURN
*
 9999 FORMAT( ' DCHKHS: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
 9998 FORMAT( ' DCHKHS: ', A, ' Eigenvectors from ', A, ' incorrectly ',
     $      'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
     $      'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5,
     $      ')' )
 9997 FORMAT( ' DCHKHS: Selected ', A, ' Eigenvectors from ', A,
     $      ' do not match other eigenvectors ', 9X, 'N=', I6,
     $      ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
*
*     End of DCHKHS
*
      END