From ecca781af87fc8419bf2ab0354a6a6ea1d254379 Mon Sep 17 00:00:00 2001 From: Angelika Schwarz Date: Sun, 26 Jun 2022 18:06:54 +0100 Subject: [PATCH 1/2] Closes #217 --- SRC/csyswapr.f | 43 +++++++++++++------------------------------ SRC/dsyswapr.f | 47 +++++++++++++++-------------------------------- SRC/ssyswapr.f | 47 +++++++++++++++-------------------------------- SRC/zsyswapr.f | 47 +++++++++++++++-------------------------------- 4 files changed, 58 insertions(+), 126 deletions(-) diff --git a/SRC/csyswapr.f b/SRC/csyswapr.f index 185d819225..04004f3c11 100644 --- a/SRC/csyswapr.f +++ b/SRC/csyswapr.f @@ -58,15 +58,13 @@ *> \param[in,out] A *> \verbatim *> A is COMPLEX array, dimension (LDA,N) -*> On entry, the NB diagonal matrix D and the multipliers -*> used to obtain the factor U or L as computed by CSYTRF. -*> -*> On exit, if INFO = 0, the (symmetric) inverse of the original -*> matrix. If UPLO = 'U', the upper triangular part of the -*> inverse is formed and the part of A below the diagonal is not -*> referenced; if UPLO = 'L' the lower triangular part of the -*> inverse is formed and the part of A above the diagonal is -*> not referenced. +*> On entry, the N-by-N matrix A. On exit, the permuted matrix +*> where the rows I1 and I2 and columns I1 and I2 are interchanged. +*> If UPLO = 'U', the interchanges are applied to the upper +*> triangular part and the strictly lower triangular part of A is +*> not referenced; if UPLO = 'L', the interchanges are applied to +*> the lower triangular part and the part of A above the diagonal +*> is not referenced. *> \endverbatim *> *> \param[in] LDA @@ -116,7 +114,6 @@ SUBROUTINE CSYSWAPR( UPLO, N, A, LDA, I1, I2) * .. * .. Local Scalars .. LOGICAL UPPER - INTEGER I COMPLEX TMP * * .. External Functions .. @@ -143,19 +140,12 @@ SUBROUTINE CSYSWAPR( UPLO, N, A, LDA, I1, I2) A(I1,I1)=A(I2,I2) A(I2,I2)=TMP * - DO I=1,I2-I1-1 - TMP=A(I1,I1+I) - A(I1,I1+I)=A(I1+I,I2) - A(I1+I,I2)=TMP - END DO + CALL CSWAP( I2-I1-1, A(I1,I1+1), LDA, A(I1+1,I2), 1 ) * * third swap * - swap row I1 and I2 from I2+1 to N - DO I=I2+1,N - TMP=A(I1,I) - A(I1,I)=A(I2,I) - A(I2,I)=TMP - END DO + IF ( I2.LT.N ) + $ CALL CSWAP( N-I2, A(I1,I2+1), LDA, A(I2,I2+1), LDA ) * ELSE * @@ -171,19 +161,12 @@ SUBROUTINE CSYSWAPR( UPLO, N, A, LDA, I1, I2) A(I1,I1)=A(I2,I2) A(I2,I2)=TMP * - DO I=1,I2-I1-1 - TMP=A(I1+I,I1) - A(I1+I,I1)=A(I2,I1+I) - A(I2,I1+I)=TMP - END DO + CALL CSWAP( I2-I1-1, A(I1+1,I1), 1, A(I2,I1+1), LDA ) * * third swap * - swap col I1 and I2 from I2+1 to N - DO I=I2+1,N - TMP=A(I,I1) - A(I,I1)=A(I,I2) - A(I,I2)=TMP - END DO + IF ( I2.LT.N ) + $ CALL CSWAP( N-I2, A(I2+1,I1), 1, A(I2+1,I2), 1 ) * ENDIF END SUBROUTINE CSYSWAPR diff --git a/SRC/dsyswapr.f b/SRC/dsyswapr.f index c60ccbefc3..93f6195f22 100644 --- a/SRC/dsyswapr.f +++ b/SRC/dsyswapr.f @@ -57,16 +57,14 @@ *> *> \param[in,out] A *> \verbatim -*> A is DOUBLE PRECISION array, dimension (LDA,N) -*> On entry, the NB diagonal matrix D and the multipliers -*> used to obtain the factor U or L as computed by DSYTRF. -*> -*> On exit, if INFO = 0, the (symmetric) inverse of the original -*> matrix. If UPLO = 'U', the upper triangular part of the -*> inverse is formed and the part of A below the diagonal is not -*> referenced; if UPLO = 'L' the lower triangular part of the -*> inverse is formed and the part of A above the diagonal is -*> not referenced. +*> A is DOUBLE PRECISION array, dimension (LDA,*) +*> On entry, the N-by-N matrix A. On exit, the permuted matrix +*> where the rows I1 and I2 and columns I1 and I2 are interchanged. +*> If UPLO = 'U', the interchanges are applied to the upper +*> triangular part and the strictly lower triangular part of A is +*> not referenced; if UPLO = 'L', the interchanges are applied to +*> the lower triangular part and the part of A above the diagonal +*> is not referenced. *> \endverbatim *> *> \param[in] LDA @@ -109,14 +107,13 @@ SUBROUTINE DSYSWAPR( UPLO, N, A, LDA, I1, I2) INTEGER I1, I2, LDA, N * .. * .. Array Arguments .. - DOUBLE PRECISION A( LDA, N ) + DOUBLE PRECISION A( LDA, * ) * * ===================================================================== * * .. * .. Local Scalars .. LOGICAL UPPER - INTEGER I DOUBLE PRECISION TMP * * .. External Functions .. @@ -143,19 +140,12 @@ SUBROUTINE DSYSWAPR( UPLO, N, A, LDA, I1, I2) A(I1,I1)=A(I2,I2) A(I2,I2)=TMP * - DO I=1,I2-I1-1 - TMP=A(I1,I1+I) - A(I1,I1+I)=A(I1+I,I2) - A(I1+I,I2)=TMP - END DO + CALL DSWAP( I2-I1-1, A(I1,I1+1), LDA, A(I1+1,I2), 1 ) * * third swap * - swap row I1 and I2 from I2+1 to N - DO I=I2+1,N - TMP=A(I1,I) - A(I1,I)=A(I2,I) - A(I2,I)=TMP - END DO + IF ( I2.LT.N ) + $ CALL DSWAP( N-I2, A(I1,I2+1), LDA, A(I2,I2+1), LDA ) * ELSE * @@ -171,19 +161,12 @@ SUBROUTINE DSYSWAPR( UPLO, N, A, LDA, I1, I2) A(I1,I1)=A(I2,I2) A(I2,I2)=TMP * - DO I=1,I2-I1-1 - TMP=A(I1+I,I1) - A(I1+I,I1)=A(I2,I1+I) - A(I2,I1+I)=TMP - END DO + CALL DSWAP( I2-I1-1, A(I1+1,I1), 1, A(I2,I1+1), LDA ) * * third swap * - swap col I1 and I2 from I2+1 to N - DO I=I2+1,N - TMP=A(I,I1) - A(I,I1)=A(I,I2) - A(I,I2)=TMP - END DO + IF ( I2.LT.N ) + $ CALL DSWAP( N-I2, A(I2+1,I1), 1, A(I2+1,I2), 1 ) * ENDIF END SUBROUTINE DSYSWAPR diff --git a/SRC/ssyswapr.f b/SRC/ssyswapr.f index 5e4265d7a8..e1ab5a22ac 100644 --- a/SRC/ssyswapr.f +++ b/SRC/ssyswapr.f @@ -57,16 +57,14 @@ *> *> \param[in,out] A *> \verbatim -*> A is REAL array, dimension (LDA,N) -*> On entry, the NB diagonal matrix D and the multipliers -*> used to obtain the factor U or L as computed by SSYTRF. -*> -*> On exit, if INFO = 0, the (symmetric) inverse of the original -*> matrix. If UPLO = 'U', the upper triangular part of the -*> inverse is formed and the part of A below the diagonal is not -*> referenced; if UPLO = 'L' the lower triangular part of the -*> inverse is formed and the part of A above the diagonal is -*> not referenced. +*> A is REAL array, dimension (LDA,*) +*> On entry, the N-by-N matrix A. On exit, the permuted matrix +*> where the rows I1 and I2 and columns I1 and I2 are interchanged. +*> If UPLO = 'U', the interchanges are applied to the upper +*> triangular part and the strictly lower triangular part of A is +*> not referenced; if UPLO = 'L', the interchanges are applied to +*> the lower triangular part and the part of A above the diagonal +*> is not referenced. *> \endverbatim *> *> \param[in] LDA @@ -109,14 +107,13 @@ SUBROUTINE SSYSWAPR( UPLO, N, A, LDA, I1, I2) INTEGER I1, I2, LDA, N * .. * .. Array Arguments .. - REAL A( LDA, N ) + REAL A( LDA, * ) * * ===================================================================== * * .. * .. Local Scalars .. LOGICAL UPPER - INTEGER I REAL TMP * * .. External Functions .. @@ -143,19 +140,12 @@ SUBROUTINE SSYSWAPR( UPLO, N, A, LDA, I1, I2) A(I1,I1)=A(I2,I2) A(I2,I2)=TMP * - DO I=1,I2-I1-1 - TMP=A(I1,I1+I) - A(I1,I1+I)=A(I1+I,I2) - A(I1+I,I2)=TMP - END DO + CALL SSWAP( I2-I1-1, A(I1,I1+1), LDA, A(I1+1,I2), 1 ) * * third swap * - swap row I1 and I2 from I2+1 to N - DO I=I2+1,N - TMP=A(I1,I) - A(I1,I)=A(I2,I) - A(I2,I)=TMP - END DO + IF ( I2.LT.N ) + $ CALL SSWAP( N-I2, A(I1,I2+1), LDA, A(I2,I2+1), LDA ) * ELSE * @@ -171,19 +161,12 @@ SUBROUTINE SSYSWAPR( UPLO, N, A, LDA, I1, I2) A(I1,I1)=A(I2,I2) A(I2,I2)=TMP * - DO I=1,I2-I1-1 - TMP=A(I1+I,I1) - A(I1+I,I1)=A(I2,I1+I) - A(I2,I1+I)=TMP - END DO + CALL SSWAP( I2-I1-1, A(I1+1,I1), 1, A(I2,I1+1), LDA ) * * third swap * - swap col I1 and I2 from I2+1 to N - DO I=I2+1,N - TMP=A(I,I1) - A(I,I1)=A(I,I2) - A(I,I2)=TMP - END DO + IF ( I2.LT.N ) + $ CALL SSWAP( N-I2, A(I2+1,I1), 1, A(I2+1,I2), 1 ) * ENDIF END SUBROUTINE SSYSWAPR diff --git a/SRC/zsyswapr.f b/SRC/zsyswapr.f index 1f1a878574..eb3c98c34b 100644 --- a/SRC/zsyswapr.f +++ b/SRC/zsyswapr.f @@ -57,16 +57,14 @@ *> *> \param[in,out] A *> \verbatim -*> A is COMPLEX*16 array, dimension (LDA,N) -*> On entry, the NB diagonal matrix D and the multipliers -*> used to obtain the factor U or L as computed by ZSYTRF. -*> -*> On exit, if INFO = 0, the (symmetric) inverse of the original -*> matrix. If UPLO = 'U', the upper triangular part of the -*> inverse is formed and the part of A below the diagonal is not -*> referenced; if UPLO = 'L' the lower triangular part of the -*> inverse is formed and the part of A above the diagonal is -*> not referenced. +*> A is COMPLEX*16 array, dimension (LDA,*) +*> On entry, the N-by-N matrix A. On exit, the permuted matrix +*> where the rows I1 and I2 and columns I1 and I2 are interchanged. +*> If UPLO = 'U', the interchanges are applied to the upper +*> triangular part and the strictly lower triangular part of A is +*> not referenced; if UPLO = 'L', the interchanges are applied to +*> the lower triangular part and the part of A above the diagonal +*> is not referenced. *> \endverbatim *> *> \param[in] LDA @@ -109,14 +107,13 @@ SUBROUTINE ZSYSWAPR( UPLO, N, A, LDA, I1, I2) INTEGER I1, I2, LDA, N * .. * .. Array Arguments .. - COMPLEX*16 A( LDA, N ) + COMPLEX*16 A( LDA, * ) * * ===================================================================== * * .. * .. Local Scalars .. LOGICAL UPPER - INTEGER I COMPLEX*16 TMP * * .. External Functions .. @@ -143,19 +140,12 @@ SUBROUTINE ZSYSWAPR( UPLO, N, A, LDA, I1, I2) A(I1,I1)=A(I2,I2) A(I2,I2)=TMP * - DO I=1,I2-I1-1 - TMP=A(I1,I1+I) - A(I1,I1+I)=A(I1+I,I2) - A(I1+I,I2)=TMP - END DO + CALL ZSWAP( I2-I1-1, A(I1,I1+1), LDA, A(I1+1,I2), 1 ) * * third swap * - swap row I1 and I2 from I2+1 to N - DO I=I2+1,N - TMP=A(I1,I) - A(I1,I)=A(I2,I) - A(I2,I)=TMP - END DO + IF ( I2.LT.N ) + $ CALL ZSWAP( N-I2, A(I1,I2+1), LDA, A(I2,I2+1), LDA ) * ELSE * @@ -171,19 +161,12 @@ SUBROUTINE ZSYSWAPR( UPLO, N, A, LDA, I1, I2) A(I1,I1)=A(I2,I2) A(I2,I2)=TMP * - DO I=1,I2-I1-1 - TMP=A(I1+I,I1) - A(I1+I,I1)=A(I2,I1+I) - A(I2,I1+I)=TMP - END DO + CALL ZSWAP( I2-I1-1, A(I1+1,I1), 1, A(I2,I1+1), LDA ) * * third swap * - swap col I1 and I2 from I2+1 to N - DO I=I2+1,N - TMP=A(I,I1) - A(I,I1)=A(I,I2) - A(I,I2)=TMP - END DO + IF ( I2.LT.N ) + $ CALL ZSWAP( N-I2, A(I2+1,I1), 1, A(I2+1,I2), 1 ) * ENDIF END SUBROUTINE ZSYSWAPR From c8a5cf510e5c3b16441dec8a9c20f7ff3a679efa Mon Sep 17 00:00:00 2001 From: Angelika Schwarz Date: Sat, 25 Jun 2022 20:22:24 +0100 Subject: [PATCH 2/2] Fix out-of-bounds write in [ds]get40 The test driver allocates a scalar for INFO, but the test writes to 3 entries. Revise INFO allocation & propagation: * Allocate sufficient space for the two INFOs * Instead of discarding INFO computed in [ds]get40, return INFO to test driver * Fix documentation of input/output arguments [ds]get31: Fix typo in docs --- TESTING/EIG/dchkec.f | 8 ++++---- TESTING/EIG/dget31.f | 2 +- TESTING/EIG/dget40.f | 23 ++++++++++++----------- TESTING/EIG/schkec.f | 8 ++++---- TESTING/EIG/sget31.f | 2 +- TESTING/EIG/sget40.f | 25 +++++++++++++------------ 6 files changed, 35 insertions(+), 33 deletions(-) diff --git a/TESTING/EIG/dchkec.f b/TESTING/EIG/dchkec.f index 8549618842..fbdf924c8c 100644 --- a/TESTING/EIG/dchkec.f +++ b/TESTING/EIG/dchkec.f @@ -92,14 +92,14 @@ SUBROUTINE DCHKEC( THRESH, TSTERR, NIN, NOUT ) INTEGER KLAEXC, KLALN2, KLANV2, KLAQTR, KLASY2, KTREXC, $ KTRSEN, KTRSNA, KTRSYL, LLAEXC, LLALN2, LLANV2, $ LLAQTR, LLASY2, LTREXC, LTRSYL, NLANV2, NLAQTR, - $ NLASY2, NTESTS, NTRSYL, KTGEXC, NTGEXC, LTGEXC + $ NLASY2, NTESTS, NTRSYL, KTGEXC, LTGEXC DOUBLE PRECISION EPS, RLAEXC, RLALN2, RLANV2, RLAQTR, RLASY2, $ RTREXC, RTRSYL, SFMIN, RTGEXC * .. * .. Local Arrays .. INTEGER LTRSEN( 3 ), LTRSNA( 3 ), NLAEXC( 2 ), - $ NLALN2( 2 ), NTREXC( 3 ), NTRSEN( 3 ), - $ NTRSNA( 3 ) + $ NLALN2( 2 ), NTGEXC( 2 ), NTREXC( 3 ), + $ NTRSEN( 3 ), NTRSNA( 3 ) DOUBLE PRECISION RTRSEN( 3 ), RTRSNA( 3 ) * .. * .. External Subroutines .. @@ -227,7 +227,7 @@ SUBROUTINE DCHKEC( THRESH, TSTERR, NIN, NOUT ) 9987 FORMAT( ' Routines pass computational tests if test ratio is les', $ 's than', F8.2, / / ) 9986 FORMAT( ' Error in DTGEXC: RMAX =', D12.3, / ' LMAX = ', I8, ' N', - $ 'INFO=', I8, ' KNT=', I8 ) + $ 'INFO=', 2I8, ' KNT=', I8 ) * * End of DCHKEC * diff --git a/TESTING/EIG/dget31.f b/TESTING/EIG/dget31.f index b952a91d24..d1e83114bb 100644 --- a/TESTING/EIG/dget31.f +++ b/TESTING/EIG/dget31.f @@ -65,7 +65,7 @@ *> *> \param[out] NINFO *> \verbatim -*> NINFO is INTEGER array, dimension (3) +*> NINFO is INTEGER array, dimension (2) *> NINFO(1) = number of examples with INFO less than 0 *> NINFO(2) = number of examples with INFO greater than 0 *> \endverbatim diff --git a/TESTING/EIG/dget40.f b/TESTING/EIG/dget40.f index 4d087f5240..0c3926022a 100644 --- a/TESTING/EIG/dget40.f +++ b/TESTING/EIG/dget40.f @@ -15,7 +15,7 @@ * DOUBLE PRECISION RMAX * .. * .. Array Arguments .. -* INTEGER NINFO( 3 ) +* INTEGER NINFO( 2 ) * * *> \par Purpose: @@ -53,8 +53,9 @@ *> *> \param[out] NINFO *> \verbatim -*> NINFO is INTEGER(3) -*> Number of examples where INFO is nonzero. +*> NINFO is INTEGER array, dimension (2) +*> NINFO( 1 ) = DTGEXC without accumulation returned INFO nonzero +*> NINFO( 2 ) = DTGEXC with accumulation returned INFO nonzero *> \endverbatim *> *> \param[out] KNT @@ -63,9 +64,10 @@ *> Total number of examples tested. *> \endverbatim *> -*> \param[out] NIN +*> \param[in] NIN *> \verbatim -*> NINFO is INTEGER +*> NIN is INTEGER +*> Input logical unit number. *> \endverbatim * * Authors: @@ -90,7 +92,7 @@ SUBROUTINE DGET40( RMAX, LMAX, NINFO, KNT, NIN ) DOUBLE PRECISION RMAX * .. * .. Array Arguments .. - INTEGER NINFO( 3 ) + INTEGER NINFO( 2 ) * .. * * ===================================================================== @@ -103,7 +105,7 @@ SUBROUTINE DGET40( RMAX, LMAX, NINFO, KNT, NIN ) * .. * .. Local Scalars .. INTEGER I, IFST, IFST1, IFST2, IFSTSV, ILST, ILST1, - $ ILST2, ILSTSV, INFO1, INFO2, J, LOC, N + $ ILST2, ILSTSV, J, LOC, N DOUBLE PRECISION EPS, RES * .. * .. Local Arrays .. @@ -130,7 +132,6 @@ SUBROUTINE DGET40( RMAX, LMAX, NINFO, KNT, NIN ) KNT = 0 NINFO( 1 ) = 0 NINFO( 2 ) = 0 - NINFO( 3 ) = 0 * * Read input data until N=0 * @@ -164,7 +165,7 @@ SUBROUTINE DGET40( RMAX, LMAX, NINFO, KNT, NIN ) CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDT ) CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDT ) CALL DTGEXC( .FALSE., .FALSE., N, T1, LDT, S1, LDT, Q, LDT, - $ Z, LDT, IFST1, ILST1, WORK, LWORK, INFO1 ) + $ Z, LDT, IFST1, ILST1, WORK, LWORK, NINFO ( 1 ) ) DO 40 I = 1, N DO 30 J = 1, N IF( I.EQ.J .AND. Q( I, J ).NE.ONE ) @@ -183,7 +184,7 @@ SUBROUTINE DGET40( RMAX, LMAX, NINFO, KNT, NIN ) CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDT ) CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDT ) CALL DTGEXC( .TRUE., .TRUE., N, T2, LDT, S2, LDT, Q, LDT, - $ Z, LDT, IFST2, ILST2, WORK, LWORK, INFO2 ) + $ Z, LDT, IFST2, ILST2, WORK, LWORK, NINFO ( 2 ) ) * * Compare T1 with T2 and S1 with S2 * @@ -199,7 +200,7 @@ SUBROUTINE DGET40( RMAX, LMAX, NINFO, KNT, NIN ) $ RES = RES + ONE / EPS IF( ILST1.NE.ILST2 ) $ RES = RES + ONE / EPS - IF( INFO1.NE.INFO2 ) + IF( NINFO( 1 ).NE.NINFO( 2 ) ) $ RES = RES + ONE / EPS * * Test orthogonality of Q and Z and backward error on T2 and S2 diff --git a/TESTING/EIG/schkec.f b/TESTING/EIG/schkec.f index e6123e1ad9..f742c5b36e 100644 --- a/TESTING/EIG/schkec.f +++ b/TESTING/EIG/schkec.f @@ -92,14 +92,14 @@ SUBROUTINE SCHKEC( THRESH, TSTERR, NIN, NOUT ) INTEGER KLAEXC, KLALN2, KLANV2, KLAQTR, KLASY2, KTREXC, $ KTRSEN, KTRSNA, KTRSYL, LLAEXC, LLALN2, LLANV2, $ LLAQTR, LLASY2, LTREXC, LTRSYL, NLANV2, NLAQTR, - $ NLASY2, NTESTS, NTRSYL, KTGEXC, NTGEXC, LTGEXC + $ NLASY2, NTESTS, NTRSYL, KTGEXC, LTGEXC REAL EPS, RLAEXC, RLALN2, RLANV2, RLAQTR, RLASY2, $ RTREXC, RTRSYL, SFMIN, RTGEXC * .. * .. Local Arrays .. INTEGER LTRSEN( 3 ), LTRSNA( 3 ), NLAEXC( 2 ), - $ NLALN2( 2 ), NTREXC( 3 ), NTRSEN( 3 ), - $ NTRSNA( 3 ) + $ NLALN2( 2 ), NTGEXC( 2 ), NTREXC( 3 ), + $ NTRSEN( 3 ), NTRSNA( 3 ) REAL RTRSEN( 3 ), RTRSNA( 3 ) * .. * .. External Subroutines .. @@ -227,7 +227,7 @@ SUBROUTINE SCHKEC( THRESH, TSTERR, NIN, NOUT ) 9987 FORMAT( ' Routines pass computational tests if test ratio is les', $ 's than', F8.2, / / ) 9986 FORMAT( ' Error in STGEXC: RMAX =', E12.3, / ' LMAX = ', I8, ' N', - $ 'INFO=', I8, ' KNT=', I8 ) + $ 'INFO=', 2I8, ' KNT=', I8 ) * * End of SCHKEC * diff --git a/TESTING/EIG/sget31.f b/TESTING/EIG/sget31.f index e6214007be..fd2d980d3b 100644 --- a/TESTING/EIG/sget31.f +++ b/TESTING/EIG/sget31.f @@ -65,7 +65,7 @@ *> *> \param[out] NINFO *> \verbatim -*> NINFO is INTEGER array, dimension (3) +*> NINFO is INTEGER array, dimension (2) *> NINFO(1) = number of examples with INFO less than 0 *> NINFO(2) = number of examples with INFO greater than 0 *> \endverbatim diff --git a/TESTING/EIG/sget40.f b/TESTING/EIG/sget40.f index 1fd3d400dd..15a25fbf1a 100644 --- a/TESTING/EIG/sget40.f +++ b/TESTING/EIG/sget40.f @@ -12,10 +12,10 @@ * * .. Scalar Arguments .. * INTEGER KNT, LMAX, NIN -* REAL RMAX +* REAL RMAX * .. * .. Array Arguments .. -* INTEGER NINFO( 3 ) +* INTEGER NINFO( 2 ) * * *> \par Purpose: @@ -53,8 +53,9 @@ *> *> \param[out] NINFO *> \verbatim -*> NINFO is INTEGER -*> Number of examples where INFO is nonzero. +*> NINFO is INTEGER array, dimension (2) +*> NINFO( 1 ) = STGEXC without accumulation returned INFO nonzero +*> NINFO( 2 ) = STGEXC with accumulation returned INFO nonzero *> \endverbatim *> *> \param[out] KNT @@ -63,9 +64,10 @@ *> Total number of examples tested. *> \endverbatim *> -*> \param[out] NIN +*> \param[in] NIN *> \verbatim -*> NINFO is INTEGER +*> NIN is INTEGER +*> Input logical unit number. *> \endverbatim * * Authors: @@ -90,7 +92,7 @@ SUBROUTINE SGET40( RMAX, LMAX, NINFO, KNT, NIN ) REAL RMAX * .. * .. Array Arguments .. - INTEGER NINFO( 3 ) + INTEGER NINFO( 2 ) * .. * * ===================================================================== @@ -103,7 +105,7 @@ SUBROUTINE SGET40( RMAX, LMAX, NINFO, KNT, NIN ) * .. * .. Local Scalars .. INTEGER I, IFST, IFST1, IFST2, IFSTSV, ILST, ILST1, - $ ILST2, ILSTSV, INFO1, INFO2, J, LOC, N + $ ILST2, ILSTSV, J, LOC, N REAL EPS, RES * .. * .. Local Arrays .. @@ -130,7 +132,6 @@ SUBROUTINE SGET40( RMAX, LMAX, NINFO, KNT, NIN ) KNT = 0 NINFO( 1 ) = 0 NINFO( 2 ) = 0 - NINFO( 3 ) = 0 * * Read input data until N=0 * @@ -164,7 +165,7 @@ SUBROUTINE SGET40( RMAX, LMAX, NINFO, KNT, NIN ) CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDT ) CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDT ) CALL STGEXC( .FALSE., .FALSE., N, T1, LDT, S1, LDT, Q, LDT, - $ Z, LDT, IFST1, ILST1, WORK, LWORK, INFO1 ) + $ Z, LDT, IFST1, ILST1, WORK, LWORK, NINFO( 1 ) ) DO 40 I = 1, N DO 30 J = 1, N IF( I.EQ.J .AND. Q( I, J ).NE.ONE ) @@ -183,7 +184,7 @@ SUBROUTINE SGET40( RMAX, LMAX, NINFO, KNT, NIN ) CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDT ) CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDT ) CALL STGEXC( .TRUE., .TRUE., N, T2, LDT, S2, LDT, Q, LDT, - $ Z, LDT, IFST2, ILST2, WORK, LWORK, INFO2 ) + $ Z, LDT, IFST2, ILST2, WORK, LWORK, NINFO( 2 ) ) * * Compare T1 with T2 and S1 with S2 * @@ -199,7 +200,7 @@ SUBROUTINE SGET40( RMAX, LMAX, NINFO, KNT, NIN ) $ RES = RES + ONE / EPS IF( ILST1.NE.ILST2 ) $ RES = RES + ONE / EPS - IF( INFO1.NE.INFO2 ) + IF( NINFO( 1 ).NE.NINFO( 2 ) ) $ RES = RES + ONE / EPS * * Test orthogonality of Q and Z and backward error on T2 and S2