Merge pull request #3909 from martin-frbg/lapack796

Fix ill-conditioned matrix in LIN testsuite test_rfp (LAPACK PR 796)
This commit is contained in:
Martin Kroeker 2023-02-15 12:56:47 +01:00 committed by GitHub
commit 10be02c896
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 168 additions and 50 deletions

View File

@ -156,9 +156,10 @@
REAL RESULT( NTESTS ) REAL RESULT( NTESTS )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME
REAL SLAMCH, CLANGE REAL SLAMCH, CLANGE
COMPLEX CLARND COMPLEX CLARND
EXTERNAL SLAMCH, CLARND, CLANGE EXTERNAL SLAMCH, CLARND, CLANGE, LSAME
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL CTRTTF, CGEQRF, CGEQLF, CTFSM, CTRSM EXTERNAL CTRTTF, CGEQRF, CGEQLF, CTFSM, CTRSM
@ -222,9 +223,9 @@
* *
DO 100 IALPHA = 1, 3 DO 100 IALPHA = 1, 3
* *
IF ( IALPHA.EQ. 1) THEN IF ( IALPHA.EQ.1 ) THEN
ALPHA = ZERO ALPHA = ZERO
ELSE IF ( IALPHA.EQ. 2) THEN ELSE IF ( IALPHA.EQ.2 ) THEN
ALPHA = ONE ALPHA = ONE
ELSE ELSE
ALPHA = CLARND( 4, ISEED ) ALPHA = CLARND( 4, ISEED )
@ -263,7 +264,7 @@
* *
DO J = 1, NA DO J = 1, NA
DO I = 1, NA DO I = 1, NA
A( I, J) = CLARND( 4, ISEED ) A( I, J ) = CLARND( 4, ISEED )
END DO END DO
END DO END DO
* *
@ -276,6 +277,20 @@
CALL CGEQRF( NA, NA, A, LDA, TAU, CALL CGEQRF( NA, NA, A, LDA, TAU,
+ C_WORK_CGEQRF, LDA, + C_WORK_CGEQRF, LDA,
+ INFO ) + INFO )
*
* Forcing main diagonal of test matrix to
* be unit makes it ill-conditioned for
* some test cases
*
IF ( LSAME( DIAG, 'U' ) ) THEN
DO J = 1, NA
DO I = 1, J
A( I, J ) = A( I, J ) /
+ ( 2.0 * A( J, J ) )
END DO
END DO
END IF
*
ELSE ELSE
* *
* The case IUPLO.EQ.2 is when SIDE.EQ.'L' * The case IUPLO.EQ.2 is when SIDE.EQ.'L'
@ -285,6 +300,20 @@
CALL CGELQF( NA, NA, A, LDA, TAU, CALL CGELQF( NA, NA, A, LDA, TAU,
+ C_WORK_CGEQRF, LDA, + C_WORK_CGEQRF, LDA,
+ INFO ) + INFO )
*
* Forcing main diagonal of test matrix to
* be unit makes it ill-conditioned for
* some test cases
*
IF ( LSAME( DIAG, 'U' ) ) THEN
DO I = 1, NA
DO J = 1, I
A( I, J ) = A( I, J ) /
+ ( 2.0 * A( I, I ) )
END DO
END DO
END IF
*
END IF END IF
* *
* After the QR factorization, the diagonal * After the QR factorization, the diagonal
@ -293,7 +322,8 @@
* value 1.0E+00. * value 1.0E+00.
* *
DO J = 1, NA DO J = 1, NA
A( J, J) = A(J,J) * CLARND( 5, ISEED ) A( J, J ) = A( J, J ) *
+ CLARND( 5, ISEED )
END DO END DO
* *
* Store a copy of A in RFP format (in ARF). * Store a copy of A in RFP format (in ARF).
@ -307,8 +337,8 @@
* *
DO J = 1, N DO J = 1, N
DO I = 1, M DO I = 1, M
B1( I, J) = CLARND( 4, ISEED ) B1( I, J ) = CLARND( 4, ISEED )
B2( I, J) = B1( I, J) B2( I, J ) = B1( I, J )
END DO END DO
END DO END DO
* *
@ -331,24 +361,24 @@
* *
DO J = 1, N DO J = 1, N
DO I = 1, M DO I = 1, M
B1( I, J) = B2( I, J ) - B1( I, J ) B1( I, J ) = B2( I, J ) - B1( I, J )
END DO END DO
END DO END DO
* *
RESULT(1) = CLANGE( 'I', M, N, B1, LDA, RESULT( 1 ) = CLANGE( 'I', M, N, B1, LDA,
+ S_WORK_CLANGE ) + S_WORK_CLANGE )
* *
RESULT(1) = RESULT(1) / SQRT( EPS ) RESULT( 1 ) = RESULT( 1 ) / SQRT( EPS )
+ / MAX ( MAX( M, N), 1 ) + / MAX ( MAX( M, N ), 1 )
* *
IF( RESULT(1).GE.THRESH ) THEN IF( RESULT( 1 ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 ) THEN IF( NFAIL.EQ.0 ) THEN
WRITE( NOUT, * ) WRITE( NOUT, * )
WRITE( NOUT, FMT = 9999 ) WRITE( NOUT, FMT = 9999 )
END IF END IF
WRITE( NOUT, FMT = 9997 ) 'CTFSM', WRITE( NOUT, FMT = 9997 ) 'CTFSM',
+ CFORM, SIDE, UPLO, TRANS, DIAG, M, + CFORM, SIDE, UPLO, TRANS, DIAG, M,
+ N, RESULT(1) + N, RESULT( 1 )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
* *

View File

@ -153,8 +153,9 @@
DOUBLE PRECISION RESULT( NTESTS ) DOUBLE PRECISION RESULT( NTESTS )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME
DOUBLE PRECISION DLAMCH, DLANGE, DLARND DOUBLE PRECISION DLAMCH, DLANGE, DLARND
EXTERNAL DLAMCH, DLANGE, DLARND EXTERNAL DLAMCH, DLANGE, DLARND, LSAME
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL DTRTTF, DGEQRF, DGEQLF, DTFSM, DTRSM EXTERNAL DTRTTF, DGEQRF, DGEQLF, DTFSM, DTRSM
@ -218,9 +219,9 @@
* *
DO 100 IALPHA = 1, 3 DO 100 IALPHA = 1, 3
* *
IF ( IALPHA.EQ. 1) THEN IF ( IALPHA.EQ.1 ) THEN
ALPHA = ZERO ALPHA = ZERO
ELSE IF ( IALPHA.EQ. 2) THEN ELSE IF ( IALPHA.EQ.2 ) THEN
ALPHA = ONE ALPHA = ONE
ELSE ELSE
ALPHA = DLARND( 2, ISEED ) ALPHA = DLARND( 2, ISEED )
@ -259,7 +260,7 @@
* *
DO J = 1, NA DO J = 1, NA
DO I = 1, NA DO I = 1, NA
A( I, J) = DLARND( 2, ISEED ) A( I, J ) = DLARND( 2, ISEED )
END DO END DO
END DO END DO
* *
@ -272,6 +273,20 @@
CALL DGEQRF( NA, NA, A, LDA, TAU, CALL DGEQRF( NA, NA, A, LDA, TAU,
+ D_WORK_DGEQRF, LDA, + D_WORK_DGEQRF, LDA,
+ INFO ) + INFO )
*
* Forcing main diagonal of test matrix to
* be unit makes it ill-conditioned for
* some test cases
*
IF ( LSAME( DIAG, 'U' ) ) THEN
DO J = 1, NA
DO I = 1, J
A( I, J ) = A( I, J ) /
+ ( 2.0 * A( J, J ) )
END DO
END DO
END IF
*
ELSE ELSE
* *
* The case IUPLO.EQ.2 is when SIDE.EQ.'L' * The case IUPLO.EQ.2 is when SIDE.EQ.'L'
@ -281,6 +296,20 @@
CALL DGELQF( NA, NA, A, LDA, TAU, CALL DGELQF( NA, NA, A, LDA, TAU,
+ D_WORK_DGEQRF, LDA, + D_WORK_DGEQRF, LDA,
+ INFO ) + INFO )
*
* Forcing main diagonal of test matrix to
* be unit makes it ill-conditioned for
* some test cases
*
IF ( LSAME( DIAG, 'U' ) ) THEN
DO I = 1, NA
DO J = 1, I
A( I, J ) = A( I, J ) /
+ ( 2.0 * A( I, I ) )
END DO
END DO
END IF
*
END IF END IF
* *
* Store a copy of A in RFP format (in ARF). * Store a copy of A in RFP format (in ARF).
@ -294,8 +323,8 @@
* *
DO J = 1, N DO J = 1, N
DO I = 1, M DO I = 1, M
B1( I, J) = DLARND( 2, ISEED ) B1( I, J ) = DLARND( 2, ISEED )
B2( I, J) = B1( I, J) B2( I, J ) = B1( I, J )
END DO END DO
END DO END DO
* *
@ -318,24 +347,24 @@
* *
DO J = 1, N DO J = 1, N
DO I = 1, M DO I = 1, M
B1( I, J) = B2( I, J ) - B1( I, J ) B1( I, J ) = B2( I, J ) - B1( I, J )
END DO END DO
END DO END DO
* *
RESULT(1) = DLANGE( 'I', M, N, B1, LDA, RESULT( 1 ) = DLANGE( 'I', M, N, B1, LDA,
+ D_WORK_DLANGE ) + D_WORK_DLANGE )
* *
RESULT(1) = RESULT(1) / SQRT( EPS ) RESULT( 1 ) = RESULT( 1 ) / SQRT( EPS )
+ / MAX ( MAX( M, N), 1 ) + / MAX ( MAX( M, N ), 1 )
* *
IF( RESULT(1).GE.THRESH ) THEN IF( RESULT( 1 ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 ) THEN IF( NFAIL.EQ.0 ) THEN
WRITE( NOUT, * ) WRITE( NOUT, * )
WRITE( NOUT, FMT = 9999 ) WRITE( NOUT, FMT = 9999 )
END IF END IF
WRITE( NOUT, FMT = 9997 ) 'DTFSM', WRITE( NOUT, FMT = 9997 ) 'DTFSM',
+ CFORM, SIDE, UPLO, TRANS, DIAG, M, + CFORM, SIDE, UPLO, TRANS, DIAG, M,
+ N, RESULT(1) + N, RESULT( 1 )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
* *

View File

@ -153,8 +153,9 @@
REAL RESULT( NTESTS ) REAL RESULT( NTESTS )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME
REAL SLAMCH, SLANGE, SLARND REAL SLAMCH, SLANGE, SLARND
EXTERNAL SLAMCH, SLANGE, SLARND EXTERNAL SLAMCH, SLANGE, SLARND, LSAME
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL STRTTF, SGEQRF, SGEQLF, STFSM, STRSM EXTERNAL STRTTF, SGEQRF, SGEQLF, STFSM, STRSM
@ -218,9 +219,9 @@
* *
DO 100 IALPHA = 1, 3 DO 100 IALPHA = 1, 3
* *
IF ( IALPHA.EQ. 1) THEN IF ( IALPHA.EQ.1 ) THEN
ALPHA = ZERO ALPHA = ZERO
ELSE IF ( IALPHA.EQ. 2) THEN ELSE IF ( IALPHA.EQ.2 ) THEN
ALPHA = ONE ALPHA = ONE
ELSE ELSE
ALPHA = SLARND( 2, ISEED ) ALPHA = SLARND( 2, ISEED )
@ -259,7 +260,7 @@
* *
DO J = 1, NA DO J = 1, NA
DO I = 1, NA DO I = 1, NA
A( I, J) = SLARND( 2, ISEED ) A( I, J ) = SLARND( 2, ISEED )
END DO END DO
END DO END DO
* *
@ -272,6 +273,20 @@
CALL SGEQRF( NA, NA, A, LDA, TAU, CALL SGEQRF( NA, NA, A, LDA, TAU,
+ S_WORK_SGEQRF, LDA, + S_WORK_SGEQRF, LDA,
+ INFO ) + INFO )
*
* Forcing main diagonal of test matrix to
* be unit makes it ill-conditioned for
* some test cases
*
IF ( LSAME( DIAG, 'U' ) ) THEN
DO J = 1, NA
DO I = 1, J
A( I, J ) = A( I, J ) /
+ ( 2.0 * A( J, J ) )
END DO
END DO
END IF
*
ELSE ELSE
* *
* The case IUPLO.EQ.2 is when SIDE.EQ.'L' * The case IUPLO.EQ.2 is when SIDE.EQ.'L'
@ -281,6 +296,20 @@
CALL SGELQF( NA, NA, A, LDA, TAU, CALL SGELQF( NA, NA, A, LDA, TAU,
+ S_WORK_SGEQRF, LDA, + S_WORK_SGEQRF, LDA,
+ INFO ) + INFO )
*
* Forcing main diagonal of test matrix to
* be unit makes it ill-conditioned for
* some test cases
*
IF ( LSAME( DIAG, 'U' ) ) THEN
DO I = 1, NA
DO J = 1, I
A( I, J ) = A( I, J ) /
+ ( 2.0 * A( I, I ) )
END DO
END DO
END IF
*
END IF END IF
* *
* Store a copy of A in RFP format (in ARF). * Store a copy of A in RFP format (in ARF).
@ -294,8 +323,8 @@
* *
DO J = 1, N DO J = 1, N
DO I = 1, M DO I = 1, M
B1( I, J) = SLARND( 2, ISEED ) B1( I, J ) = SLARND( 2, ISEED )
B2( I, J) = B1( I, J) B2( I, J ) = B1( I, J )
END DO END DO
END DO END DO
* *
@ -318,24 +347,24 @@
* *
DO J = 1, N DO J = 1, N
DO I = 1, M DO I = 1, M
B1( I, J) = B2( I, J ) - B1( I, J ) B1( I, J ) = B2( I, J ) - B1( I, J )
END DO END DO
END DO END DO
* *
RESULT(1) = SLANGE( 'I', M, N, B1, LDA, RESULT( 1 ) = SLANGE( 'I', M, N, B1, LDA,
+ S_WORK_SLANGE ) + S_WORK_SLANGE )
* *
RESULT(1) = RESULT(1) / SQRT( EPS ) RESULT( 1 ) = RESULT( 1 ) / SQRT( EPS )
+ / MAX ( MAX( M, N), 1 ) + / MAX ( MAX( M, N ), 1 )
* *
IF( RESULT(1).GE.THRESH ) THEN IF( RESULT( 1 ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 ) THEN IF( NFAIL.EQ.0 ) THEN
WRITE( NOUT, * ) WRITE( NOUT, * )
WRITE( NOUT, FMT = 9999 ) WRITE( NOUT, FMT = 9999 )
END IF END IF
WRITE( NOUT, FMT = 9997 ) 'STFSM', WRITE( NOUT, FMT = 9997 ) 'STFSM',
+ CFORM, SIDE, UPLO, TRANS, DIAG, M, + CFORM, SIDE, UPLO, TRANS, DIAG, M,
+ N, RESULT(1) + N, RESULT( 1 )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
* *

View File

@ -156,9 +156,10 @@
DOUBLE PRECISION RESULT( NTESTS ) DOUBLE PRECISION RESULT( NTESTS )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME
DOUBLE PRECISION DLAMCH, ZLANGE DOUBLE PRECISION DLAMCH, ZLANGE
COMPLEX*16 ZLARND COMPLEX*16 ZLARND
EXTERNAL DLAMCH, ZLARND, ZLANGE EXTERNAL DLAMCH, ZLARND, ZLANGE, LSAME
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZTRTTF, ZGEQRF, ZGEQLF, ZTFSM, ZTRSM EXTERNAL ZTRTTF, ZGEQRF, ZGEQLF, ZTFSM, ZTRSM
@ -222,9 +223,9 @@
* *
DO 100 IALPHA = 1, 3 DO 100 IALPHA = 1, 3
* *
IF ( IALPHA.EQ. 1) THEN IF ( IALPHA.EQ.1 ) THEN
ALPHA = ZERO ALPHA = ZERO
ELSE IF ( IALPHA.EQ. 2) THEN ELSE IF ( IALPHA.EQ.2 ) THEN
ALPHA = ONE ALPHA = ONE
ELSE ELSE
ALPHA = ZLARND( 4, ISEED ) ALPHA = ZLARND( 4, ISEED )
@ -263,7 +264,7 @@
* *
DO J = 1, NA DO J = 1, NA
DO I = 1, NA DO I = 1, NA
A( I, J) = ZLARND( 4, ISEED ) A( I, J ) = ZLARND( 4, ISEED )
END DO END DO
END DO END DO
* *
@ -276,6 +277,20 @@
CALL ZGEQRF( NA, NA, A, LDA, TAU, CALL ZGEQRF( NA, NA, A, LDA, TAU,
+ Z_WORK_ZGEQRF, LDA, + Z_WORK_ZGEQRF, LDA,
+ INFO ) + INFO )
*
* Forcing main diagonal of test matrix to
* be unit makes it ill-conditioned for
* some test cases
*
IF ( LSAME( DIAG, 'U' ) ) THEN
DO J = 1, NA
DO I = 1, J
A( I, J ) = A( I, J ) /
+ ( 2.0 * A( J, J ) )
END DO
END DO
END IF
*
ELSE ELSE
* *
* The case IUPLO.EQ.2 is when SIDE.EQ.'L' * The case IUPLO.EQ.2 is when SIDE.EQ.'L'
@ -285,6 +300,20 @@
CALL ZGELQF( NA, NA, A, LDA, TAU, CALL ZGELQF( NA, NA, A, LDA, TAU,
+ Z_WORK_ZGEQRF, LDA, + Z_WORK_ZGEQRF, LDA,
+ INFO ) + INFO )
*
* Forcing main diagonal of test matrix to
* be unit makes it ill-conditioned for
* some test cases
*
IF ( LSAME( DIAG, 'U' ) ) THEN
DO I = 1, NA
DO J = 1, I
A( I, J ) = A( I, J ) /
+ ( 2.0 * A( I, I ) )
END DO
END DO
END IF
*
END IF END IF
* *
* After the QR factorization, the diagonal * After the QR factorization, the diagonal
@ -293,7 +322,8 @@
* value 1.0E+00. * value 1.0E+00.
* *
DO J = 1, NA DO J = 1, NA
A( J, J) = A(J,J) * ZLARND( 5, ISEED ) A( J, J ) = A( J, J ) *
+ ZLARND( 5, ISEED )
END DO END DO
* *
* Store a copy of A in RFP format (in ARF). * Store a copy of A in RFP format (in ARF).
@ -307,8 +337,8 @@
* *
DO J = 1, N DO J = 1, N
DO I = 1, M DO I = 1, M
B1( I, J) = ZLARND( 4, ISEED ) B1( I, J ) = ZLARND( 4, ISEED )
B2( I, J) = B1( I, J) B2( I, J ) = B1( I, J )
END DO END DO
END DO END DO
* *
@ -331,24 +361,24 @@
* *
DO J = 1, N DO J = 1, N
DO I = 1, M DO I = 1, M
B1( I, J) = B2( I, J ) - B1( I, J ) B1( I, J ) = B2( I, J ) - B1( I, J )
END DO END DO
END DO END DO
* *
RESULT(1) = ZLANGE( 'I', M, N, B1, LDA, RESULT( 1 ) = ZLANGE( 'I', M, N, B1, LDA,
+ D_WORK_ZLANGE ) + D_WORK_ZLANGE )
* *
RESULT(1) = RESULT(1) / SQRT( EPS ) RESULT( 1 ) = RESULT( 1 ) / SQRT( EPS )
+ / MAX ( MAX( M, N), 1 ) + / MAX ( MAX( M, N ), 1 )
* *
IF( RESULT(1).GE.THRESH ) THEN IF( RESULT( 1 ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 ) THEN IF( NFAIL.EQ.0 ) THEN
WRITE( NOUT, * ) WRITE( NOUT, * )
WRITE( NOUT, FMT = 9999 ) WRITE( NOUT, FMT = 9999 )
END IF END IF
WRITE( NOUT, FMT = 9997 ) 'ZTFSM', WRITE( NOUT, FMT = 9997 ) 'ZTFSM',
+ CFORM, SIDE, UPLO, TRANS, DIAG, M, + CFORM, SIDE, UPLO, TRANS, DIAG, M,
+ N, RESULT(1) + N, RESULT( 1 )
NFAIL = NFAIL + 1 NFAIL = NFAIL + 1
END IF END IF
* *