Skip to content

Add numerical tests for trevc3 #682

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Mar 18, 2023
Merged

Add numerical tests for trevc3 #682

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
22d1721
Add missing numerical tests for TREVC3
angsch Jun 19, 2022
9f9295f
standardize style in laqr5
angsch Jul 3, 2022
0a6cd43
Rewrite [ds]hgeqz to use FMA with Householder reflectors
angsch Jul 5, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
Add missing numerical tests for TREVC3 
At least some tests, though there are still code paths
that are not covered
* input sizes defined in nep.in are small
* RWORK in [CZ]TREVC3 is de factor defined as N-vector
  from the input file and limits the blocked computation
  • Loading branch information
@angsch
angsch committed Jul 6, 2022
commit 22d172188ce2c1cab440ea4568b5a061a46e557b
81 changes: 75 additions & 6 deletions TESTING/EIG/cchkhs.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
* .. Array Arguments ..
* LOGICAL DOTYPE( * ), SELECT( * )
* INTEGER ISEED( 4 ), IWORK( * ), NN( * )
* REAL RESULT( 14 ), RWORK( * )
* REAL RESULT( 16 ), RWORK( * )
* COMPLEX A( LDA, * ), EVECTL( LDU, * ),
* $ EVECTR( LDU, * ), EVECTX( LDU, * ),
* $ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ),
Expand Down Expand Up @@ -64,10 +64,15 @@
*> eigenvectors of H. Y is lower triangular, and X is
*> 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
*> number of matrix "types" are specified. For each size ("n")
*> 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:
*>
*> (1) | A - U H U**H | / ( |A| n ulp )
Expand Down Expand Up @@ -98,6 +103,10 @@
*>
*> (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
*> each element NN(j) specifies one size.
*> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
Expand Down Expand Up @@ -331,7 +340,7 @@
*> Workspace. Could be equivalenced to IWORK, but not RWORK.
*> Modified.
*>
*> RESULT - REAL array, dimension (14)
*> RESULT - REAL array, dimension (16)
*> The values computed by the fourteen tests described above.
*> The values are currently limited to 1/ulp, to avoid
*> overflow.
Expand Down Expand Up @@ -421,7 +430,7 @@ SUBROUTINE CCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
* .. Array Arguments ..
LOGICAL DOTYPE( * ), SELECT( * )
INTEGER ISEED( 4 ), IWORK( * ), NN( * )
REAL RESULT( 14 ), RWORK( * )
REAL RESULT( 16 ), RWORK( * )
COMPLEX A( LDA, * ), EVECTL( LDU, * ),
$ EVECTR( LDU, * ), EVECTX( LDU, * ),
$ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ),
Expand Down Expand Up @@ -463,8 +472,8 @@ SUBROUTINE CCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
* .. External Subroutines ..
EXTERNAL CCOPY, CGEHRD, CGEMM, CGET10, CGET22, CHSEIN,
$ CHSEQR, CHST01, CLACPY, CLASET, CLATME, CLATMR,
$ CLATMS, CTREVC, CUNGHR, CUNMHR, SLABAD, SLAFTS,
$ SLASUM, XERBLA
$ CLATMS, CTREVC, CTREVC3, CUNGHR, CUNMHR,
$ SLABAD, SLAFTS, SLASUM, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, REAL, SQRT
Expand Down Expand Up @@ -1067,6 +1076,66 @@ SUBROUTINE CCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
$ RESULT( 14 ) = DUMMA( 3 )*ANINV
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
*
240 CONTINUE
Expand Down
82 changes: 75 additions & 7 deletions TESTING/EIG/dchkhs.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
* INTEGER ISEED( 4 ), IWORK( * ), NN( * )
* DOUBLE PRECISION A( LDA, * ), EVECTL( LDU, * ),
* $ EVECTR( LDU, * ), EVECTX( LDU, * ),
* $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ),
* $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 16 ),
* $ T1( LDA, * ), T2( LDA, * ), TAU( * ),
* $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ),
* $ WI1( * ), WI2( * ), WI3( * ), WORK( * ),
Expand All @@ -49,15 +49,21 @@
*> T is "quasi-triangular", and the eigenvalue vector W.
*>
*> 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
*> 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
*> number of matrix "types" are specified. For each size ("n")
*> 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:
*>
*> (1) | A - U H U**T | / ( |A| n ulp )
Expand Down Expand Up @@ -88,6 +94,10 @@
*>
*> (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
*> each element NN(j) specifies one size.
*> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
Expand Down Expand Up @@ -331,7 +341,7 @@
*> Workspace.
*> Modified.
*>
*> RESULT - DOUBLE PRECISION array, dimension (14)
*> RESULT - DOUBLE PRECISION array, dimension (16)
*> The values computed by the fourteen tests described above.
*> The values are currently limited to 1/ulp, to avoid
*> overflow.
Expand Down Expand Up @@ -423,7 +433,7 @@ SUBROUTINE DCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
INTEGER ISEED( 4 ), IWORK( * ), NN( * )
DOUBLE PRECISION A( LDA, * ), EVECTL( LDU, * ),
$ EVECTR( LDU, * ), EVECTX( LDU, * ),
$ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ),
$ EVECTY( LDU, * ), H( LDA, * ), RESULT( 16 ),
$ T1( LDA, * ), T2( LDA, * ), TAU( * ),
$ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ),
$ WI1( * ), WI2( * ), WI3( * ), WORK( * ),
Expand Down Expand Up @@ -461,7 +471,7 @@ SUBROUTINE DCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
EXTERNAL DCOPY, DGEHRD, DGEMM, DGET10, DGET22, DHSEIN,
$ DHSEQR, DHST01, DLABAD, DLACPY, DLAFTS, DLASET,
$ DLASUM, DLATME, DLATMR, DLATMS, DORGHR, DORMHR,
$ DTREVC, XERBLA
$ DTREVC, DTREVC3, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, MIN, SQRT
Expand Down Expand Up @@ -561,7 +571,7 @@ SUBROUTINE DCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
*
* Initialize RESULT
*
DO 30 J = 1, 14
DO 30 J = 1, 16
RESULT( J ) = ZERO
30 CONTINUE
*
Expand Down Expand Up @@ -1108,6 +1118,64 @@ SUBROUTINE DCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
$ RESULT( 14 ) = DUMMA( 3 )*ANINV
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
*
250 CONTINUE
Expand Down
Loading