Add numerical tests for TRECV3 (Reference-LAPACK 682)

This commit is contained in:
Martin Kroeker 2023-03-20 10:04:05 +01:00 committed by GitHub
parent 2a83ec1f79
commit 147e2fbf87
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 297 additions and 24 deletions

View File

@ -21,7 +21,7 @@
* .. Array Arguments .. * .. Array Arguments ..
* LOGICAL DOTYPE( * ), SELECT( * ) * LOGICAL DOTYPE( * ), SELECT( * )
* INTEGER ISEED( 4 ), IWORK( * ), NN( * ) * INTEGER ISEED( 4 ), IWORK( * ), NN( * )
* REAL RESULT( 14 ), RWORK( * ) * REAL RESULT( 16 ), RWORK( * )
* COMPLEX A( LDA, * ), EVECTL( LDU, * ), * COMPLEX A( LDA, * ), EVECTL( LDU, * ),
* $ EVECTR( LDU, * ), EVECTX( LDU, * ), * $ EVECTR( LDU, * ), EVECTX( LDU, * ),
* $ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ), * $ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ),
@ -64,10 +64,15 @@
*> eigenvectors of H. Y is lower triangular, and X is *> eigenvectors of H. Y is lower triangular, and X is
*> upper triangular. *> upper triangular.
*> *>
*> CTREVC3 computes left and right eigenvector matrices
*> from a Schur matrix T and backtransforms them with Z
*> to eigenvector matrices L and R for A. L and R are
*> GE matrices.
*>
*> When CCHKHS is called, a number of matrix "sizes" ("n's") and a *> When CCHKHS is called, a number of matrix "sizes" ("n's") and a
*> number of matrix "types" are specified. For each size ("n") *> number of matrix "types" are specified. For each size ("n")
*> and each type of matrix, one matrix will be generated and used *> and each type of matrix, one matrix will be generated and used
*> to test the nonsymmetric eigenroutines. For each matrix, 14 *> to test the nonsymmetric eigenroutines. For each matrix, 16
*> tests will be performed: *> tests will be performed:
*> *>
*> (1) | A - U H U**H | / ( |A| n ulp ) *> (1) | A - U H U**H | / ( |A| n ulp )
@ -98,6 +103,10 @@
*> *>
*> (14) | Y**H A - W**H Y | / ( |A| |Y| ulp ) *> (14) | Y**H A - W**H Y | / ( |A| |Y| ulp )
*> *>
*> (15) | AR - RW | / ( |A| |R| ulp )
*>
*> (16) | LA - WL | / ( |A| |L| ulp )
*>
*> The "sizes" are specified by an array NN(1:NSIZES); the value of *> The "sizes" are specified by an array NN(1:NSIZES); the value of
*> each element NN(j) specifies one size. *> each element NN(j) specifies one size.
*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
@ -331,7 +340,7 @@
*> Workspace. Could be equivalenced to IWORK, but not RWORK. *> Workspace. Could be equivalenced to IWORK, but not RWORK.
*> Modified. *> Modified.
*> *>
*> RESULT - REAL array, dimension (14) *> RESULT - REAL array, dimension (16)
*> The values computed by the fourteen tests described above. *> The values computed by the fourteen tests described above.
*> The values are currently limited to 1/ulp, to avoid *> The values are currently limited to 1/ulp, to avoid
*> overflow. *> overflow.
@ -421,7 +430,7 @@
* .. Array Arguments .. * .. Array Arguments ..
LOGICAL DOTYPE( * ), SELECT( * ) LOGICAL DOTYPE( * ), SELECT( * )
INTEGER ISEED( 4 ), IWORK( * ), NN( * ) INTEGER ISEED( 4 ), IWORK( * ), NN( * )
REAL RESULT( 14 ), RWORK( * ) REAL RESULT( 16 ), RWORK( * )
COMPLEX A( LDA, * ), EVECTL( LDU, * ), COMPLEX A( LDA, * ), EVECTL( LDU, * ),
$ EVECTR( LDU, * ), EVECTX( LDU, * ), $ EVECTR( LDU, * ), EVECTX( LDU, * ),
$ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ), $ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ),
@ -463,8 +472,8 @@
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL CCOPY, CGEHRD, CGEMM, CGET10, CGET22, CHSEIN, EXTERNAL CCOPY, CGEHRD, CGEMM, CGET10, CGET22, CHSEIN,
$ CHSEQR, CHST01, CLACPY, CLASET, CLATME, CLATMR, $ CHSEQR, CHST01, CLACPY, CLASET, CLATME, CLATMR,
$ CLATMS, CTREVC, CUNGHR, CUNMHR, SLABAD, SLAFTS, $ CLATMS, CTREVC, CTREVC3, CUNGHR, CUNMHR,
$ SLASUM, XERBLA $ SLABAD, SLAFTS, SLASUM, XERBLA
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, REAL, SQRT INTRINSIC ABS, MAX, MIN, REAL, SQRT
@ -1067,6 +1076,66 @@
$ RESULT( 14 ) = DUMMA( 3 )*ANINV $ RESULT( 14 ) = DUMMA( 3 )*ANINV
END IF END IF
* *
* Compute Left and Right Eigenvectors of A
*
* Compute a Right eigenvector matrix:
*
NTEST = 15
RESULT( 15 ) = ULPINV
*
CALL CLACPY( ' ', N, N, UZ, LDU, EVECTR, LDU )
*
CALL CTREVC3( 'Right', 'Back', SELECT, N, T1, LDA, CDUMMA,
$ LDU, EVECTR, LDU, N, IN, WORK, NWORK, RWORK,
$ N, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CTREVC3(R,B)', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
GO TO 250
END IF
*
* Test 15: | AR - RW | / ( |A| |R| ulp )
*
* (from Schur decomposition)
*
CALL CGET22( 'N', 'N', 'N', N, A, LDA, EVECTR, LDU, W1,
$ WORK, RWORK, DUMMA( 1 ) )
RESULT( 15 ) = DUMMA( 1 )
IF( DUMMA( 2 ).GT.THRESH ) THEN
WRITE( NOUNIT, FMT = 9998 )'Right', 'CTREVC3',
$ DUMMA( 2 ), N, JTYPE, IOLDSD
END IF
*
* Compute a Left eigenvector matrix:
*
NTEST = 16
RESULT( 16 ) = ULPINV
*
CALL CLACPY( ' ', N, N, UZ, LDU, EVECTL, LDU )
*
CALL CTREVC3( 'Left', 'Back', SELECT, N, T1, LDA, EVECTL,
$ LDU, CDUMMA, LDU, N, IN, WORK, NWORK, RWORK,
$ N, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'CTREVC3(L,B)', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
GO TO 250
END IF
*
* Test 16: | LA - WL | / ( |A| |L| ulp )
*
* (from Schur decomposition)
*
CALL CGET22( 'Conj', 'N', 'Conj', N, A, LDA, EVECTL, LDU,
$ W1, WORK, RWORK, DUMMA( 3 ) )
RESULT( 16 ) = DUMMA( 3 )
IF( DUMMA( 4 ).GT.THRESH ) THEN
WRITE( NOUNIT, FMT = 9998 )'Left', 'CTREVC3', DUMMA( 4 ),
$ N, JTYPE, IOLDSD
END IF
*
* End of Loop -- Check for RESULT(j) > THRESH * End of Loop -- Check for RESULT(j) > THRESH
* *
240 CONTINUE 240 CONTINUE

View File

@ -23,7 +23,7 @@
* INTEGER ISEED( 4 ), IWORK( * ), NN( * ) * INTEGER ISEED( 4 ), IWORK( * ), NN( * )
* DOUBLE PRECISION A( LDA, * ), EVECTL( LDU, * ), * DOUBLE PRECISION A( LDA, * ), EVECTL( LDU, * ),
* $ EVECTR( LDU, * ), EVECTX( LDU, * ), * $ EVECTR( LDU, * ), EVECTX( LDU, * ),
* $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ), * $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 16 ),
* $ T1( LDA, * ), T2( LDA, * ), TAU( * ), * $ T1( LDA, * ), T2( LDA, * ), TAU( * ),
* $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ), * $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ),
* $ WI1( * ), WI2( * ), WI3( * ), WORK( * ), * $ WI1( * ), WI2( * ), WI3( * ), WORK( * ),
@ -49,15 +49,21 @@
*> T is "quasi-triangular", and the eigenvalue vector W. *> T is "quasi-triangular", and the eigenvalue vector W.
*> *>
*> DTREVC computes the left and right eigenvector matrices *> DTREVC computes the left and right eigenvector matrices
*> L and R for T. *> L and R for T. L is lower quasi-triangular, and R is
*> upper quasi-triangular.
*> *>
*> DHSEIN computes the left and right eigenvector matrices *> DHSEIN computes the left and right eigenvector matrices
*> Y and X for H, using inverse iteration. *> Y and X for H, using inverse iteration.
*> *>
*> DTREVC3 computes left and right eigenvector matrices
*> from a Schur matrix T and backtransforms them with Z
*> to eigenvector matrices L and R for A. L and R are
*> GE matrices.
*>
*> When DCHKHS is called, a number of matrix "sizes" ("n's") and a *> When DCHKHS is called, a number of matrix "sizes" ("n's") and a
*> number of matrix "types" are specified. For each size ("n") *> number of matrix "types" are specified. For each size ("n")
*> and each type of matrix, one matrix will be generated and used *> and each type of matrix, one matrix will be generated and used
*> to test the nonsymmetric eigenroutines. For each matrix, 14 *> to test the nonsymmetric eigenroutines. For each matrix, 16
*> tests will be performed: *> tests will be performed:
*> *>
*> (1) | A - U H U**T | / ( |A| n ulp ) *> (1) | A - U H U**T | / ( |A| n ulp )
@ -88,6 +94,10 @@
*> *>
*> (14) | Y**H A - W**H Y | / ( |A| |Y| ulp ) *> (14) | Y**H A - W**H Y | / ( |A| |Y| ulp )
*> *>
*> (15) | AR - RW | / ( |A| |R| ulp )
*>
*> (16) | LA - WL | / ( |A| |L| ulp )
*>
*> The "sizes" are specified by an array NN(1:NSIZES); the value of *> The "sizes" are specified by an array NN(1:NSIZES); the value of
*> each element NN(j) specifies one size. *> each element NN(j) specifies one size.
*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
@ -331,7 +341,7 @@
*> Workspace. *> Workspace.
*> Modified. *> Modified.
*> *>
*> RESULT - DOUBLE PRECISION array, dimension (14) *> RESULT - DOUBLE PRECISION array, dimension (16)
*> The values computed by the fourteen tests described above. *> The values computed by the fourteen tests described above.
*> The values are currently limited to 1/ulp, to avoid *> The values are currently limited to 1/ulp, to avoid
*> overflow. *> overflow.
@ -423,7 +433,7 @@
INTEGER ISEED( 4 ), IWORK( * ), NN( * ) INTEGER ISEED( 4 ), IWORK( * ), NN( * )
DOUBLE PRECISION A( LDA, * ), EVECTL( LDU, * ), DOUBLE PRECISION A( LDA, * ), EVECTL( LDU, * ),
$ EVECTR( LDU, * ), EVECTX( LDU, * ), $ EVECTR( LDU, * ), EVECTX( LDU, * ),
$ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ), $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 16 ),
$ T1( LDA, * ), T2( LDA, * ), TAU( * ), $ T1( LDA, * ), T2( LDA, * ), TAU( * ),
$ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ), $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ),
$ WI1( * ), WI2( * ), WI3( * ), WORK( * ), $ WI1( * ), WI2( * ), WI3( * ), WORK( * ),
@ -461,7 +471,7 @@
EXTERNAL DCOPY, DGEHRD, DGEMM, DGET10, DGET22, DHSEIN, EXTERNAL DCOPY, DGEHRD, DGEMM, DGET10, DGET22, DHSEIN,
$ DHSEQR, DHST01, DLABAD, DLACPY, DLAFTS, DLASET, $ DHSEQR, DHST01, DLABAD, DLACPY, DLAFTS, DLASET,
$ DLASUM, DLATME, DLATMR, DLATMS, DORGHR, DORMHR, $ DLASUM, DLATME, DLATMR, DLATMS, DORGHR, DORMHR,
$ DTREVC, XERBLA $ DTREVC, DTREVC3, XERBLA
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, MIN, SQRT INTRINSIC ABS, DBLE, MAX, MIN, SQRT
@ -561,7 +571,7 @@
* *
* Initialize RESULT * Initialize RESULT
* *
DO 30 J = 1, 14 DO 30 J = 1, 16
RESULT( J ) = ZERO RESULT( J ) = ZERO
30 CONTINUE 30 CONTINUE
* *
@ -1108,6 +1118,64 @@
$ RESULT( 14 ) = DUMMA( 3 )*ANINV $ RESULT( 14 ) = DUMMA( 3 )*ANINV
END IF END IF
* *
* Compute Left and Right Eigenvectors of A
*
* Compute a Right eigenvector matrix:
*
NTEST = 15
RESULT( 15 ) = ULPINV
*
CALL DLACPY( ' ', N, N, UZ, LDU, EVECTR, LDU )
*
CALL DTREVC3( 'Right', 'Back', SELECT, N, T1, LDA, DUMMA,
$ LDU, EVECTR, LDU, N, IN, WORK, NWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'DTREVC3(R,B)', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
GO TO 250
END IF
*
* Test 15: | AR - RW | / ( |A| |R| ulp )
*
* (from Schur decomposition)
*
CALL DGET22( 'N', 'N', 'N', N, A, LDA, EVECTR, LDU, WR1,
$ WI1, WORK, DUMMA( 1 ) )
RESULT( 15 ) = DUMMA( 1 )
IF( DUMMA( 2 ).GT.THRESH ) THEN
WRITE( NOUNIT, FMT = 9998 )'Right', 'DTREVC3',
$ DUMMA( 2 ), N, JTYPE, IOLDSD
END IF
*
* Compute a Left eigenvector matrix:
*
NTEST = 16
RESULT( 16 ) = ULPINV
*
CALL DLACPY( ' ', N, N, UZ, LDU, EVECTL, LDU )
*
CALL DTREVC3( 'Left', 'Back', SELECT, N, T1, LDA, EVECTL,
$ LDU, DUMMA, LDU, N, IN, WORK, NWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'DTREVC3(L,B)', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
GO TO 250
END IF
*
* Test 16: | LA - WL | / ( |A| |L| ulp )
*
* (from Schur decomposition)
*
CALL DGET22( 'Trans', 'N', 'Conj', N, A, LDA, EVECTL, LDU,
$ WR1, WI1, WORK, DUMMA( 3 ) )
RESULT( 16 ) = DUMMA( 3 )
IF( DUMMA( 4 ).GT.THRESH ) THEN
WRITE( NOUNIT, FMT = 9998 )'Left', 'DTREVC3', DUMMA( 4 ),
$ N, JTYPE, IOLDSD
END IF
*
* End of Loop -- Check for RESULT(j) > THRESH * End of Loop -- Check for RESULT(j) > THRESH
* *
250 CONTINUE 250 CONTINUE

View File

@ -23,7 +23,7 @@
* INTEGER ISEED( 4 ), IWORK( * ), NN( * ) * INTEGER ISEED( 4 ), IWORK( * ), NN( * )
* REAL A( LDA, * ), EVECTL( LDU, * ), * REAL A( LDA, * ), EVECTL( LDU, * ),
* $ EVECTR( LDU, * ), EVECTX( LDU, * ), * $ EVECTR( LDU, * ), EVECTX( LDU, * ),
* $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ), * $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 16 ),
* $ T1( LDA, * ), T2( LDA, * ), TAU( * ), * $ T1( LDA, * ), T2( LDA, * ), TAU( * ),
* $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ), * $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ),
* $ WI1( * ), WI2( * ), WI3( * ), WORK( * ), * $ WI1( * ), WI2( * ), WI3( * ), WORK( * ),
@ -54,10 +54,15 @@
*> SHSEIN computes the left and right eigenvector matrices *> SHSEIN computes the left and right eigenvector matrices
*> Y and X for H, using inverse iteration. *> Y and X for H, using inverse iteration.
*> *>
*> STREVC3 computes left and right eigenvector matrices
*> from a Schur matrix T and backtransforms them with Z
*> to eigenvector matrices L and R for A. L and R are
*> GE matrices.
*>
*> When SCHKHS is called, a number of matrix "sizes" ("n's") and a *> When SCHKHS is called, a number of matrix "sizes" ("n's") and a
*> number of matrix "types" are specified. For each size ("n") *> number of matrix "types" are specified. For each size ("n")
*> and each type of matrix, one matrix will be generated and used *> and each type of matrix, one matrix will be generated and used
*> to test the nonsymmetric eigenroutines. For each matrix, 14 *> to test the nonsymmetric eigenroutines. For each matrix, 16
*> tests will be performed: *> tests will be performed:
*> *>
*> (1) | A - U H U**T | / ( |A| n ulp ) *> (1) | A - U H U**T | / ( |A| n ulp )
@ -88,6 +93,10 @@
*> *>
*> (14) | Y**H A - W**H Y | / ( |A| |Y| ulp ) *> (14) | Y**H A - W**H Y | / ( |A| |Y| ulp )
*> *>
*> (15) | AR - RW | / ( |A| |R| ulp )
*>
*> (16) | LA - WL | / ( |A| |L| ulp )
*>
*> The "sizes" are specified by an array NN(1:NSIZES); the value of *> The "sizes" are specified by an array NN(1:NSIZES); the value of
*> each element NN(j) specifies one size. *> each element NN(j) specifies one size.
*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
@ -331,7 +340,7 @@
*> Workspace. *> Workspace.
*> Modified. *> Modified.
*> *>
*> RESULT - REAL array, dimension (14) *> RESULT - REAL array, dimension (16)
*> The values computed by the fourteen tests described above. *> The values computed by the fourteen tests described above.
*> The values are currently limited to 1/ulp, to avoid *> The values are currently limited to 1/ulp, to avoid
*> overflow. *> overflow.
@ -423,7 +432,7 @@
INTEGER ISEED( 4 ), IWORK( * ), NN( * ) INTEGER ISEED( 4 ), IWORK( * ), NN( * )
REAL A( LDA, * ), EVECTL( LDU, * ), REAL A( LDA, * ), EVECTL( LDU, * ),
$ EVECTR( LDU, * ), EVECTX( LDU, * ), $ EVECTR( LDU, * ), EVECTX( LDU, * ),
$ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ), $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 16 ),
$ T1( LDA, * ), T2( LDA, * ), TAU( * ), $ T1( LDA, * ), T2( LDA, * ), TAU( * ),
$ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ), $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ),
$ WI1( * ), WI2( * ), WI3( * ), WORK( * ), $ WI1( * ), WI2( * ), WI3( * ), WORK( * ),
@ -461,7 +470,7 @@
EXTERNAL SCOPY, SGEHRD, SGEMM, SGET10, SGET22, SHSEIN, EXTERNAL SCOPY, SGEHRD, SGEMM, SGET10, SGET22, SHSEIN,
$ SHSEQR, SHST01, SLABAD, SLACPY, SLAFTS, SLASET, $ SHSEQR, SHST01, SLABAD, SLACPY, SLAFTS, SLASET,
$ SLASUM, SLATME, SLATMR, SLATMS, SORGHR, SORMHR, $ SLASUM, SLATME, SLATMR, SLATMS, SORGHR, SORMHR,
$ STREVC, XERBLA $ STREVC, STREVC3, XERBLA
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, REAL, SQRT INTRINSIC ABS, MAX, MIN, REAL, SQRT
@ -561,7 +570,7 @@
* *
* Initialize RESULT * Initialize RESULT
* *
DO 30 J = 1, 14 DO 30 J = 1, 16
RESULT( J ) = ZERO RESULT( J ) = ZERO
30 CONTINUE 30 CONTINUE
* *
@ -1108,6 +1117,64 @@
$ RESULT( 14 ) = DUMMA( 3 )*ANINV $ RESULT( 14 ) = DUMMA( 3 )*ANINV
END IF END IF
* *
* Compute Left and Right Eigenvectors of A
*
* Compute a Right eigenvector matrix:
*
NTEST = 15
RESULT( 15 ) = ULPINV
*
CALL SLACPY( ' ', N, N, UZ, LDU, EVECTR, LDU )
*
CALL STREVC3( 'Right', 'Back', SELECT, N, T1, LDA, DUMMA,
$ LDU, EVECTR, LDU, N, IN, WORK, NWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'STREVC3(R,B)', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
GO TO 250
END IF
*
* Test 15: | AR - RW | / ( |A| |R| ulp )
*
* (from Schur decomposition)
*
CALL SGET22( 'N', 'N', 'N', N, A, LDA, EVECTR, LDU, WR1,
$ WI1, WORK, DUMMA( 1 ) )
RESULT( 15 ) = DUMMA( 1 )
IF( DUMMA( 2 ).GT.THRESH ) THEN
WRITE( NOUNIT, FMT = 9998 )'Right', 'STREVC3',
$ DUMMA( 2 ), N, JTYPE, IOLDSD
END IF
*
* Compute a Left eigenvector matrix:
*
NTEST = 16
RESULT( 16 ) = ULPINV
*
CALL SLACPY( ' ', N, N, UZ, LDU, EVECTL, LDU )
*
CALL STREVC3( 'Left', 'Back', SELECT, N, T1, LDA, EVECTL,
$ LDU, DUMMA, LDU, N, IN, WORK, NWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'STREVC3(L,B)', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
GO TO 250
END IF
*
* Test 16: | LA - WL | / ( |A| |L| ulp )
*
* (from Schur decomposition)
*
CALL SGET22( 'Trans', 'N', 'Conj', N, A, LDA, EVECTL, LDU,
$ WR1, WI1, WORK, DUMMA( 3 ) )
RESULT( 16 ) = DUMMA( 3 )
IF( DUMMA( 4 ).GT.THRESH ) THEN
WRITE( NOUNIT, FMT = 9998 )'Left', 'STREVC3', DUMMA( 4 ),
$ N, JTYPE, IOLDSD
END IF
*
* End of Loop -- Check for RESULT(j) > THRESH * End of Loop -- Check for RESULT(j) > THRESH
* *
250 CONTINUE 250 CONTINUE

View File

@ -21,7 +21,7 @@
* .. Array Arguments .. * .. Array Arguments ..
* LOGICAL DOTYPE( * ), SELECT( * ) * LOGICAL DOTYPE( * ), SELECT( * )
* INTEGER ISEED( 4 ), IWORK( * ), NN( * ) * INTEGER ISEED( 4 ), IWORK( * ), NN( * )
* DOUBLE PRECISION RESULT( 14 ), RWORK( * ) * DOUBLE PRECISION RESULT( 16 ), RWORK( * )
* COMPLEX*16 A( LDA, * ), EVECTL( LDU, * ), * COMPLEX*16 A( LDA, * ), EVECTL( LDU, * ),
* $ EVECTR( LDU, * ), EVECTX( LDU, * ), * $ EVECTR( LDU, * ), EVECTX( LDU, * ),
* $ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ), * $ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ),
@ -64,10 +64,15 @@
*> eigenvectors of H. Y is lower triangular, and X is *> eigenvectors of H. Y is lower triangular, and X is
*> upper triangular. *> upper triangular.
*> *>
*> ZTREVC3 computes left and right eigenvector matrices
*> from a Schur matrix T and backtransforms them with Z
*> to eigenvector matrices L and R for A. L and R are
*> GE matrices.
*>
*> When ZCHKHS is called, a number of matrix "sizes" ("n's") and a *> When ZCHKHS is called, a number of matrix "sizes" ("n's") and a
*> number of matrix "types" are specified. For each size ("n") *> number of matrix "types" are specified. For each size ("n")
*> and each type of matrix, one matrix will be generated and used *> and each type of matrix, one matrix will be generated and used
*> to test the nonsymmetric eigenroutines. For each matrix, 14 *> to test the nonsymmetric eigenroutines. For each matrix, 16
*> tests will be performed: *> tests will be performed:
*> *>
*> (1) | A - U H U**H | / ( |A| n ulp ) *> (1) | A - U H U**H | / ( |A| n ulp )
@ -98,6 +103,10 @@
*> *>
*> (14) | Y**H A - W**H Y | / ( |A| |Y| ulp ) *> (14) | Y**H A - W**H Y | / ( |A| |Y| ulp )
*> *>
*> (15) | AR - RW | / ( |A| |R| ulp )
*>
*> (16) | LA - WL | / ( |A| |L| ulp )
*>
*> The "sizes" are specified by an array NN(1:NSIZES); the value of *> The "sizes" are specified by an array NN(1:NSIZES); the value of
*> each element NN(j) specifies one size. *> each element NN(j) specifies one size.
*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
@ -331,7 +340,7 @@
*> Workspace. Could be equivalenced to IWORK, but not RWORK. *> Workspace. Could be equivalenced to IWORK, but not RWORK.
*> Modified. *> Modified.
*> *>
*> RESULT - DOUBLE PRECISION array, dimension (14) *> RESULT - DOUBLE PRECISION array, dimension (16)
*> The values computed by the fourteen tests described above. *> The values computed by the fourteen tests described above.
*> The values are currently limited to 1/ulp, to avoid *> The values are currently limited to 1/ulp, to avoid
*> overflow. *> overflow.
@ -421,7 +430,7 @@
* .. Array Arguments .. * .. Array Arguments ..
LOGICAL DOTYPE( * ), SELECT( * ) LOGICAL DOTYPE( * ), SELECT( * )
INTEGER ISEED( 4 ), IWORK( * ), NN( * ) INTEGER ISEED( 4 ), IWORK( * ), NN( * )
DOUBLE PRECISION RESULT( 14 ), RWORK( * ) DOUBLE PRECISION RESULT( 16 ), RWORK( * )
COMPLEX*16 A( LDA, * ), EVECTL( LDU, * ), COMPLEX*16 A( LDA, * ), EVECTL( LDU, * ),
$ EVECTR( LDU, * ), EVECTX( LDU, * ), $ EVECTR( LDU, * ), EVECTX( LDU, * ),
$ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ), $ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ),
@ -464,7 +473,7 @@
EXTERNAL DLABAD, DLAFTS, DLASUM, XERBLA, ZCOPY, ZGEHRD, EXTERNAL DLABAD, DLAFTS, DLASUM, XERBLA, ZCOPY, ZGEHRD,
$ ZGEMM, ZGET10, ZGET22, ZHSEIN, ZHSEQR, ZHST01, $ ZGEMM, ZGET10, ZGET22, ZHSEIN, ZHSEQR, ZHST01,
$ ZLACPY, ZLASET, ZLATME, ZLATMR, ZLATMS, ZTREVC, $ ZLACPY, ZLASET, ZLATME, ZLATMR, ZLATMS, ZTREVC,
$ ZUNGHR, ZUNMHR $ ZTREVC3, ZUNGHR, ZUNMHR
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, MIN, SQRT INTRINSIC ABS, DBLE, MAX, MIN, SQRT
@ -1067,6 +1076,66 @@
$ RESULT( 14 ) = DUMMA( 3 )*ANINV $ RESULT( 14 ) = DUMMA( 3 )*ANINV
END IF END IF
* *
* Compute Left and Right Eigenvectors of A
*
* Compute a Right eigenvector matrix:
*
NTEST = 15
RESULT( 15 ) = ULPINV
*
CALL ZLACPY( ' ', N, N, UZ, LDU, EVECTR, LDU )
*
CALL ZTREVC3( 'Right', 'Back', SELECT, N, T1, LDA, CDUMMA,
$ LDU, EVECTR, LDU, N, IN, WORK, NWORK, RWORK,
$ N, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZTREVC3(R,B)', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
GO TO 250
END IF
*
* Test 15: | AR - RW | / ( |A| |R| ulp )
*
* (from Schur decomposition)
*
CALL ZGET22( 'N', 'N', 'N', N, A, LDA, EVECTR, LDU, W1,
$ WORK, RWORK, DUMMA( 1 ) )
RESULT( 15 ) = DUMMA( 1 )
IF( DUMMA( 2 ).GT.THRESH ) THEN
WRITE( NOUNIT, FMT = 9998 )'Right', 'ZTREVC3',
$ DUMMA( 2 ), N, JTYPE, IOLDSD
END IF
*
* Compute a Left eigenvector matrix:
*
NTEST = 16
RESULT( 16 ) = ULPINV
*
CALL ZLACPY( ' ', N, N, UZ, LDU, EVECTL, LDU )
*
CALL ZTREVC3( 'Left', 'Back', SELECT, N, T1, LDA, EVECTL,
$ LDU, CDUMMA, LDU, N, IN, WORK, NWORK, RWORK,
$ N, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUNIT, FMT = 9999 )'ZTREVC3(L,B)', IINFO, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
GO TO 250
END IF
*
* Test 16: | LA - WL | / ( |A| |L| ulp )
*
* (from Schur decomposition)
*
CALL ZGET22( 'Conj', 'N', 'Conj', N, A, LDA, EVECTL, LDU,
$ W1, WORK, RWORK, DUMMA( 3 ) )
RESULT( 16 ) = DUMMA( 3 )
IF( DUMMA( 4 ).GT.THRESH ) THEN
WRITE( NOUNIT, FMT = 9998 )'Left', 'ZTREVC3', DUMMA( 4 ),
$ N, JTYPE, IOLDSD
END IF
*
* End of Loop -- Check for RESULT(j) > THRESH * End of Loop -- Check for RESULT(j) > THRESH
* *
240 CONTINUE 240 CONTINUE