From ff7c5ca217e1d3a7a9c27b2ef97847fc4d84784a Mon Sep 17 00:00:00 2001 From: Igor Zhuravlov Date: Sun, 30 May 2021 17:52:28 +1000 Subject: [PATCH 1/4] comments fixed: RWORK dimension The RWORK dimension was understated. --- TESTING/LIN/cchkgb.f | 2 +- TESTING/LIN/cdrvgb.f | 2 +- TESTING/LIN/cdrvgbx.f | 2 +- TESTING/LIN/dchkgb.f | 2 +- TESTING/LIN/ddrvgb.f | 2 +- TESTING/LIN/ddrvgbx.f | 2 +- TESTING/LIN/schkgb.f | 2 +- TESTING/LIN/sdrvgb.f | 2 +- TESTING/LIN/sdrvgbx.f | 2 +- TESTING/LIN/zchkgb.f | 2 +- TESTING/LIN/zdrvgb.f | 2 +- TESTING/LIN/zdrvgbx.f | 2 +- 12 files changed, 12 insertions(+), 12 deletions(-) diff --git a/TESTING/LIN/cchkgb.f b/TESTING/LIN/cchkgb.f index 433b45f4ec..e8a12c8a50 100644 --- a/TESTING/LIN/cchkgb.f +++ b/TESTING/LIN/cchkgb.f @@ -160,7 +160,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is REAL array, dimension -*> (max(NMAX,2*NSMAX)) +*> (NMAX+2*NSMAX) *> \endverbatim *> *> \param[out] IWORK diff --git a/TESTING/LIN/cdrvgb.f b/TESTING/LIN/cdrvgb.f index b0a8be002f..867669523b 100644 --- a/TESTING/LIN/cdrvgb.f +++ b/TESTING/LIN/cdrvgb.f @@ -141,7 +141,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is REAL array, dimension -*> (max(NMAX,2*NRHS)) +*> (NMAX+2*NRHS) *> \endverbatim *> *> \param[out] IWORK diff --git a/TESTING/LIN/cdrvgbx.f b/TESTING/LIN/cdrvgbx.f index 747cb40adb..1f1f664dbe 100644 --- a/TESTING/LIN/cdrvgbx.f +++ b/TESTING/LIN/cdrvgbx.f @@ -144,7 +144,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is REAL array, dimension -*> (max(NMAX,2*NRHS)) +*> (max(2*NMAX,NMAX+2*NRHS)) *> \endverbatim *> *> \param[out] IWORK diff --git a/TESTING/LIN/dchkgb.f b/TESTING/LIN/dchkgb.f index c8a7efeb76..f4b7685648 100644 --- a/TESTING/LIN/dchkgb.f +++ b/TESTING/LIN/dchkgb.f @@ -160,7 +160,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is DOUBLE PRECISION array, dimension -*> (max(NMAX,2*NSMAX)) +*> (NMAX+2*NSMAX) *> \endverbatim *> *> \param[out] IWORK diff --git a/TESTING/LIN/ddrvgb.f b/TESTING/LIN/ddrvgb.f index ed0e18a850..802ff7354b 100644 --- a/TESTING/LIN/ddrvgb.f +++ b/TESTING/LIN/ddrvgb.f @@ -141,7 +141,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is DOUBLE PRECISION array, dimension -*> (max(NMAX,2*NRHS)) +*> (NMAX+2*NRHS) *> \endverbatim *> *> \param[out] IWORK diff --git a/TESTING/LIN/ddrvgbx.f b/TESTING/LIN/ddrvgbx.f index b8236703ad..949e39fd8a 100644 --- a/TESTING/LIN/ddrvgbx.f +++ b/TESTING/LIN/ddrvgbx.f @@ -144,7 +144,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is DOUBLE PRECISION array, dimension -*> (max(NMAX,2*NRHS)) +*> (max(2*NMAX,NMAX+2*NRHS)) *> \endverbatim *> *> \param[out] IWORK diff --git a/TESTING/LIN/schkgb.f b/TESTING/LIN/schkgb.f index 7eee800aa4..1418e1afbe 100644 --- a/TESTING/LIN/schkgb.f +++ b/TESTING/LIN/schkgb.f @@ -160,7 +160,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is REAL array, dimension -*> (max(NMAX,2*NSMAX)) +*> (NMAX+2*NSMAX) *> \endverbatim *> *> \param[out] IWORK diff --git a/TESTING/LIN/sdrvgb.f b/TESTING/LIN/sdrvgb.f index 8c65b6ae9b..1dc318b154 100644 --- a/TESTING/LIN/sdrvgb.f +++ b/TESTING/LIN/sdrvgb.f @@ -141,7 +141,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is REAL array, dimension -*> (max(NMAX,2*NRHS)) +*> (NMAX+2*NRHS) *> \endverbatim *> *> \param[out] IWORK diff --git a/TESTING/LIN/sdrvgbx.f b/TESTING/LIN/sdrvgbx.f index a38c3486a4..e8cb538fa9 100644 --- a/TESTING/LIN/sdrvgbx.f +++ b/TESTING/LIN/sdrvgbx.f @@ -144,7 +144,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is REAL array, dimension -*> (max(NMAX,2*NRHS)) +*> (max(2*NMAX,NMAX+2*NRHS)) *> \endverbatim *> *> \param[out] IWORK diff --git a/TESTING/LIN/zchkgb.f b/TESTING/LIN/zchkgb.f index 39de25197c..6e002bf187 100644 --- a/TESTING/LIN/zchkgb.f +++ b/TESTING/LIN/zchkgb.f @@ -160,7 +160,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is DOUBLE PRECISION array, dimension -*> (max(NMAX,2*NSMAX)) +*> (NMAX+2*NSMAX) *> \endverbatim *> *> \param[out] IWORK diff --git a/TESTING/LIN/zdrvgb.f b/TESTING/LIN/zdrvgb.f index 3431b4b31d..5ee8cd8195 100644 --- a/TESTING/LIN/zdrvgb.f +++ b/TESTING/LIN/zdrvgb.f @@ -141,7 +141,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is DOUBLE PRECISION array, dimension -*> (max(NMAX,2*NRHS)) +*> (NMAX+2*NRHS) *> \endverbatim *> *> \param[out] IWORK diff --git a/TESTING/LIN/zdrvgbx.f b/TESTING/LIN/zdrvgbx.f index cd0e077549..1a437daa69 100644 --- a/TESTING/LIN/zdrvgbx.f +++ b/TESTING/LIN/zdrvgbx.f @@ -144,7 +144,7 @@ *> \param[out] RWORK *> \verbatim *> RWORK is DOUBLE PRECISION array, dimension -*> (max(NMAX,2*NRHS)) +*> (max(2*NMAX,NMAX+2*NRHS)) *> \endverbatim *> *> \param[out] IWORK From 705bb8129640e478d8e0e0b2ee053e091388c18b Mon Sep 17 00:00:00 2001 From: Igor Zhuravlov Date: Sun, 30 May 2021 18:05:43 +1000 Subject: [PATCH 2/4] code fixed: check added to avoid subscript out of bounds error If M=1 or M Date: Sun, 30 May 2021 18:30:19 +1000 Subject: [PATCH 3/4] modifies xGBT02 and xGET02 to use 1-norm of op(A) not A (no.2 from PR #562) 1. Comments improved, unified and fixed. 2. Code changed. To unify with routines xGBT01, xGTT02, xTBT02, xTPT02, and xTRT02. --- TESTING/EIG/cget02.f | 25 ++++++++------ TESTING/EIG/dget02.f | 25 ++++++++------ TESTING/EIG/sget02.f | 25 ++++++++------ TESTING/EIG/zget02.f | 25 ++++++++------ TESTING/LIN/cchkgb.f | 4 +-- TESTING/LIN/cdrvgb.f | 4 ++- TESTING/LIN/cdrvgbx.f | 7 ++-- TESTING/LIN/cgbt02.f | 76 ++++++++++++++++++++++++++++++++----------- TESTING/LIN/cget02.f | 27 ++++++++------- TESTING/LIN/dchkgb.f | 4 +-- TESTING/LIN/ddrvgb.f | 4 ++- TESTING/LIN/ddrvgbx.f | 6 ++-- TESTING/LIN/dgbt02.f | 73 +++++++++++++++++++++++++++++------------ TESTING/LIN/dget02.f | 27 ++++++++------- TESTING/LIN/schkgb.f | 4 +-- TESTING/LIN/sdrvgb.f | 4 ++- TESTING/LIN/sdrvgbx.f | 6 ++-- TESTING/LIN/sgbt02.f | 73 +++++++++++++++++++++++++++++------------ TESTING/LIN/sget02.f | 27 ++++++++------- TESTING/LIN/zchkgb.f | 4 +-- TESTING/LIN/zdrvgb.f | 4 ++- TESTING/LIN/zdrvgbx.f | 7 ++-- TESTING/LIN/zgbt02.f | 76 ++++++++++++++++++++++++++++++++----------- TESTING/LIN/zget02.f | 27 ++++++++------- 24 files changed, 380 insertions(+), 184 deletions(-) diff --git a/TESTING/EIG/cget02.f b/TESTING/EIG/cget02.f index 373a41b456..3a8f6308ca 100644 --- a/TESTING/EIG/cget02.f +++ b/TESTING/EIG/cget02.f @@ -28,9 +28,10 @@ *> \verbatim *> *> CGET02 computes the residual for a solution of a system of linear -*> equations A*x = b or A'*x = b: -*> RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ), -*> where EPS is the machine epsilon. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A, A**T, or A**H, depending on TRANS, and EPS is the +*> machine epsilon. *> \endverbatim * * Arguments: @@ -40,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A^T*x = b, where A^T is the transpose of A -*> = 'C': A^H*x = b, where A^H is the conjugate transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose) *> \endverbatim *> *> \param[in] M @@ -114,7 +115,7 @@ *> \verbatim *> RESID is REAL *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -188,19 +189,23 @@ SUBROUTINE CGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = SLAMCH( 'Epsilon' ) - ANORM = CLANGE( '1', M, N, A, LDA, RWORK ) + IF( LSAME( TRANS, 'N' ) ) THEN + ANORM = CLANGE( '1', M, N, A, LDA, RWORK ) + ELSE + ANORM = CLANGE( 'I', M, N, A, LDA, RWORK ) + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN END IF * -* Compute B - A*X (or B - A'*X ) and store in B. +* Compute B - op(A)*X and store in B. * CALL CGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X, $ LDX, CONE, B, LDB ) * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ) . * RESID = ZERO DO 10 J = 1, NRHS diff --git a/TESTING/EIG/dget02.f b/TESTING/EIG/dget02.f index 1d131dd1a1..7be23a689b 100644 --- a/TESTING/EIG/dget02.f +++ b/TESTING/EIG/dget02.f @@ -28,9 +28,10 @@ *> \verbatim *> *> DGET02 computes the residual for a solution of a system of linear -*> equations A*x = b or A'*x = b: -*> RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ), -*> where EPS is the machine epsilon. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A or A**T, depending on TRANS, and EPS is the +*> machine epsilon. *> \endverbatim * * Arguments: @@ -40,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A'*x = b, where A' is the transpose of A -*> = 'C': A'*x = b, where A' is the transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose = Transpose) *> \endverbatim *> *> \param[in] M @@ -114,7 +115,7 @@ *> \verbatim *> RESID is DOUBLE PRECISION *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -186,19 +187,23 @@ SUBROUTINE DGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = DLAMCH( 'Epsilon' ) - ANORM = DLANGE( '1', M, N, A, LDA, RWORK ) + IF( LSAME( TRANS, 'N' ) ) THEN + ANORM = DLANGE( '1', M, N, A, LDA, RWORK ) + ELSE + ANORM = DLANGE( 'I', M, N, A, LDA, RWORK ) + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN END IF * -* Compute B - A*X (or B - A'*X ) and store in B. +* Compute B - op(A)*X and store in B. * CALL DGEMM( TRANS, 'No transpose', N1, NRHS, N2, -ONE, A, LDA, X, $ LDX, ONE, B, LDB ) * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ) . * RESID = ZERO DO 10 J = 1, NRHS diff --git a/TESTING/EIG/sget02.f b/TESTING/EIG/sget02.f index 5b6e00fc21..8c1e8ea69a 100644 --- a/TESTING/EIG/sget02.f +++ b/TESTING/EIG/sget02.f @@ -28,9 +28,10 @@ *> \verbatim *> *> SGET02 computes the residual for a solution of a system of linear -*> equations A*x = b or A'*x = b: -*> RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ), -*> where EPS is the machine epsilon. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A or A**T, depending on TRANS, and EPS is the +*> machine epsilon. *> \endverbatim * * Arguments: @@ -40,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A'*x = b, where A' is the transpose of A -*> = 'C': A'*x = b, where A' is the transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose = Transpose) *> \endverbatim *> *> \param[in] M @@ -114,7 +115,7 @@ *> \verbatim *> RESID is REAL *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -186,19 +187,23 @@ SUBROUTINE SGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = SLAMCH( 'Epsilon' ) - ANORM = SLANGE( '1', M, N, A, LDA, RWORK ) + IF( LSAME( TRANS, 'N' ) ) THEN + ANORM = SLANGE( '1', M, N, A, LDA, RWORK ) + ELSE + ANORM = SLANGE( 'I', M, N, A, LDA, RWORK ) + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN END IF * -* Compute B - A*X (or B - A'*X ) and store in B. +* Compute B - op(A)*X and store in B. * CALL SGEMM( TRANS, 'No transpose', N1, NRHS, N2, -ONE, A, LDA, X, $ LDX, ONE, B, LDB ) * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ) . * RESID = ZERO DO 10 J = 1, NRHS diff --git a/TESTING/EIG/zget02.f b/TESTING/EIG/zget02.f index 39cea4f786..c43704693f 100644 --- a/TESTING/EIG/zget02.f +++ b/TESTING/EIG/zget02.f @@ -28,9 +28,10 @@ *> \verbatim *> *> ZGET02 computes the residual for a solution of a system of linear -*> equations A*x = b or A'*x = b: -*> RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ), -*> where EPS is the machine epsilon. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A, A**T, or A**H, depending on TRANS, and EPS is the +*> machine epsilon. *> \endverbatim * * Arguments: @@ -40,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A^T*x = b, where A^T is the transpose of A -*> = 'C': A^H*x = b, where A^H is the conjugate transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose) *> \endverbatim *> *> \param[in] M @@ -114,7 +115,7 @@ *> \verbatim *> RESID is DOUBLE PRECISION *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -188,19 +189,23 @@ SUBROUTINE ZGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = DLAMCH( 'Epsilon' ) - ANORM = ZLANGE( '1', M, N, A, LDA, RWORK ) + IF( LSAME( TRANS, 'N' ) ) THEN + ANORM = ZLANGE( '1', M, N, A, LDA, RWORK ) + ELSE + ANORM = ZLANGE( 'I', M, N, A, LDA, RWORK ) + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN END IF * -* Compute B - A*X (or B - A'*X ) and store in B. +* Compute B - op(A)*X and store in B. * CALL ZGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X, $ LDX, CONE, B, LDB ) * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ) . * RESID = ZERO DO 10 J = 1, NRHS diff --git a/TESTING/LIN/cchkgb.f b/TESTING/LIN/cchkgb.f index e8a12c8a50..9bb6f29b79 100644 --- a/TESTING/LIN/cchkgb.f +++ b/TESTING/LIN/cchkgb.f @@ -563,7 +563,7 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, END IF * *+ TEST 2: -* Solve and compute residual for A * X = B. +* Solve and compute residual for op(A) * X = B. * SRNAMT = 'CLARHS' CALL CLARHS( PATH, XTYPE, ' ', TRANS, N, @@ -589,7 +589,7 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ WORK, LDB ) CALL CGBT02( TRANS, M, N, KL, KU, NRHS, A, $ LDA, X, LDB, WORK, LDB, - $ RESULT( 2 ) ) + $ RWORK, RESULT( 2 ) ) * *+ TEST 3: * Check solution from generated exact diff --git a/TESTING/LIN/cdrvgb.f b/TESTING/LIN/cdrvgb.f index 867669523b..389725eaed 100644 --- a/TESTING/LIN/cdrvgb.f +++ b/TESTING/LIN/cdrvgb.f @@ -582,7 +582,8 @@ SUBROUTINE CDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL CGBT02( 'No transpose', N, N, KL, $ KU, NRHS, A, LDA, X, LDB, - $ WORK, LDB, RESULT( 2 ) ) + $ WORK, LDB, RWORK, + $ RESULT( 2 ) ) * * Check solution from generated exact * solution. @@ -699,6 +700,7 @@ SUBROUTINE CDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL CGBT02( TRANS, N, N, KL, KU, NRHS, $ ASAV, LDA, X, LDB, WORK, LDB, + $ RWORK( 2*NRHS+1 ), $ RESULT( 2 ) ) * * Check solution from generated exact diff --git a/TESTING/LIN/cdrvgbx.f b/TESTING/LIN/cdrvgbx.f index 1f1f664dbe..68ff66e07e 100644 --- a/TESTING/LIN/cdrvgbx.f +++ b/TESTING/LIN/cdrvgbx.f @@ -590,7 +590,8 @@ SUBROUTINE CDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL CGBT02( 'No transpose', N, N, KL, $ KU, NRHS, A, LDA, X, LDB, - $ WORK, LDB, RESULT( 2 ) ) + $ WORK, LDB, RWORK, + $ RESULT( 2 ) ) * * Check solution from generated exact * solution. @@ -708,6 +709,7 @@ SUBROUTINE CDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL CGBT02( TRANS, N, N, KL, KU, NRHS, $ ASAV, LDA, X, LDB, WORK, LDB, + $ RWORK( 2*NRHS+1 ), $ RESULT( 2 ) ) * * Check solution from generated exact @@ -897,7 +899,8 @@ SUBROUTINE CDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, CALL CLACPY( 'Full', N, NRHS, BSAV, LDB, WORK, $ LDB ) CALL CGBT02( TRANS, N, N, KL, KU, NRHS, ASAV, - $ LDA, X, LDB, WORK, LDB, RESULT( 2 ) ) + $ LDA, X, LDB, WORK, LDB, RWORK, + $ RESULT( 2 ) ) * * Check solution from generated exact solution. * diff --git a/TESTING/LIN/cgbt02.f b/TESTING/LIN/cgbt02.f index e3f9d2b118..12fc57284b 100644 --- a/TESTING/LIN/cgbt02.f +++ b/TESTING/LIN/cgbt02.f @@ -9,7 +9,7 @@ * =========== * * SUBROUTINE CGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, -* LDB, RESID ) +* LDB, RWORK, RESID ) * * .. Scalar Arguments .. * CHARACTER TRANS @@ -17,6 +17,7 @@ * REAL RESID * .. * .. Array Arguments .. +* REAL RWORK( * ) * COMPLEX A( LDA, * ), B( LDB, * ), X( LDX, * ) * .. * @@ -27,9 +28,10 @@ *> \verbatim *> *> CGBT02 computes the residual for a solution of a banded system of -*> equations A*x = b or A'*x = b: -*> RESID = norm( B - A*X ) / ( norm(A) * norm(X) * EPS). -*> where EPS is the machine precision. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A, A**T, or A**H, depending on TRANS, and EPS is the +*> machine precision. *> \endverbatim * * Arguments: @@ -39,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A'*x = b, where A' is the transpose of A -*> = 'C': A'*x = b, where A' is the transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose) *> \endverbatim *> *> \param[in] M @@ -116,11 +118,18 @@ *> LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N). *> \endverbatim *> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (MAX(1,LRWORK)), +*> where LRWORK >= M when TRANS = 'T' or 'C'; otherwise, RWORK +*> is not referenced. +*> \endverbatim +* *> \param[out] RESID *> \verbatim *> RESID is REAL *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -135,7 +144,7 @@ * * ===================================================================== SUBROUTINE CGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, - $ LDB, RESID ) + $ LDB, RWORK, RESID ) * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -147,6 +156,7 @@ SUBROUTINE CGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, REAL RESID * .. * .. Array Arguments .. + REAL RWORK( * ) COMPLEX A( LDA, * ), B( LDB, * ), X( LDX, * ) * .. * @@ -161,6 +171,7 @@ SUBROUTINE CGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * .. Local Scalars .. INTEGER I1, I2, J, KD, N1 REAL ANORM, BNORM, EPS, XNORM + COMPLEX CDUM * .. * .. External Functions .. LOGICAL LSAME @@ -170,8 +181,14 @@ SUBROUTINE CGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * .. External Subroutines .. EXTERNAL CGBMV * .. +* .. Statement Functions .. + REAL CABS1 +* .. * .. Intrinsic Functions .. - INTRINSIC MAX, MIN + INTRINSIC ABS, AIMAG, MAX, MIN, REAL +* .. +* .. Statement Function definitions .. + CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) * .. * .. Executable Statements .. * @@ -185,14 +202,35 @@ SUBROUTINE CGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = SLAMCH( 'Epsilon' ) - KD = KU + 1 ANORM = ZERO - DO 10 J = 1, N - I1 = MAX( KD+1-J, 1 ) - I2 = MIN( KD+M-J, KL+KD ) - IF( I2.GE.I1 ) - $ ANORM = MAX( ANORM, SCASUM( I2-I1+1, A( I1, J ), 1 ) ) - 10 CONTINUE + IF( LSAME( TRANS, 'N' ) ) THEN +* +* Find norm1(A). +* + KD = KU + 1 + DO 10 J = 1, N + I1 = MAX( KD+1-J, 1 ) + I2 = MIN( KD+M-J, KL+KD ) + IF( I2.GE.I1 ) + $ ANORM = MAX( ANORM, SCASUM( I2-I1+1, A( I1, J ), 1 ) ) + 10 CONTINUE + ELSE +* +* Find normI(A). +* + DO 12 I1 = 1, M + RWORK( I1 ) = ZERO + 12 CONTINUE + DO 16 J = 1, N + KD = KU + 1 - J + DO 14 I1 = MAX( 1, J-KU ), MIN( M, J+KL ) + RWORK( I1 ) = RWORK( I1 ) + CABS1( A( KD+I1, J ) ) + 14 CONTINUE + 16 CONTINUE + DO 18 I1 = 1, M + ANORM = MAX( ANORM, RWORK( I1 ) ) + 18 CONTINUE + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN @@ -204,7 +242,7 @@ SUBROUTINE CGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, N1 = M END IF * -* Compute B - A*X (or B - A'*X ) +* Compute B - op(A)*X * DO 20 J = 1, NRHS CALL CGBMV( TRANS, M, N, KL, KU, -CONE, A, LDA, X( 1, J ), 1, @@ -212,7 +250,7 @@ SUBROUTINE CGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, 20 CONTINUE * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). * RESID = ZERO DO 30 J = 1, NRHS diff --git a/TESTING/LIN/cget02.f b/TESTING/LIN/cget02.f index 32f7652c63..805ad949ec 100644 --- a/TESTING/LIN/cget02.f +++ b/TESTING/LIN/cget02.f @@ -28,9 +28,10 @@ *> \verbatim *> *> CGET02 computes the residual for a solution of a system of linear -*> equations A*x = b or A'*x = b: -*> RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ), -*> where EPS is the machine epsilon. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A, A**T, or A**H, depending on TRANS, and EPS is the +*> machine epsilon. *> \endverbatim * * Arguments: @@ -40,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A^T*x = b, where A^T is the transpose of A -*> = 'C': A^H*x = b, where A^H is the conjugate transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose) *> \endverbatim *> *> \param[in] M @@ -95,7 +96,7 @@ *> B is COMPLEX array, dimension (LDB,NRHS) *> On entry, the right hand side vectors for the system of *> linear equations. -*> On exit, B is overwritten with the difference B - A*X. +*> On exit, B is overwritten with the difference B - op(A)*X. *> \endverbatim *> *> \param[in] LDB @@ -114,7 +115,7 @@ *> \verbatim *> RESID is REAL *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -188,19 +189,23 @@ SUBROUTINE CGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = SLAMCH( 'Epsilon' ) - ANORM = CLANGE( '1', M, N, A, LDA, RWORK ) + IF( LSAME( TRANS, 'N' ) ) THEN + ANORM = CLANGE( '1', M, N, A, LDA, RWORK ) + ELSE + ANORM = CLANGE( 'I', M, N, A, LDA, RWORK ) + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN END IF * -* Compute B - A*X (or B - A'*X ) and store in B. +* Compute B - op(A)*X and store in B. * CALL CGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X, $ LDX, CONE, B, LDB ) * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ) . * RESID = ZERO DO 10 J = 1, NRHS diff --git a/TESTING/LIN/dchkgb.f b/TESTING/LIN/dchkgb.f index f4b7685648..fa8dce5cff 100644 --- a/TESTING/LIN/dchkgb.f +++ b/TESTING/LIN/dchkgb.f @@ -563,7 +563,7 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, END IF * *+ TEST 2: -* Solve and compute residual for A * X = B. +* Solve and compute residual for op(A) * X = B. * SRNAMT = 'DLARHS' CALL DLARHS( PATH, XTYPE, ' ', TRANS, N, @@ -589,7 +589,7 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ WORK, LDB ) CALL DGBT02( TRANS, M, N, KL, KU, NRHS, A, $ LDA, X, LDB, WORK, LDB, - $ RESULT( 2 ) ) + $ RWORK, RESULT( 2 ) ) * *+ TEST 3: * Check solution from generated exact diff --git a/TESTING/LIN/ddrvgb.f b/TESTING/LIN/ddrvgb.f index 802ff7354b..ad323d5ec0 100644 --- a/TESTING/LIN/ddrvgb.f +++ b/TESTING/LIN/ddrvgb.f @@ -582,7 +582,8 @@ SUBROUTINE DDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL DGBT02( 'No transpose', N, N, KL, $ KU, NRHS, A, LDA, X, LDB, - $ WORK, LDB, RESULT( 2 ) ) + $ WORK, LDB, RWORK, + $ RESULT( 2 ) ) * * Check solution from generated exact * solution. @@ -699,6 +700,7 @@ SUBROUTINE DDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL DGBT02( TRANS, N, N, KL, KU, NRHS, $ ASAV, LDA, X, LDB, WORK, LDB, + $ RWORK( 2*NRHS+1 ), $ RESULT( 2 ) ) * * Check solution from generated exact diff --git a/TESTING/LIN/ddrvgbx.f b/TESTING/LIN/ddrvgbx.f index 949e39fd8a..8a9fdb645a 100644 --- a/TESTING/LIN/ddrvgbx.f +++ b/TESTING/LIN/ddrvgbx.f @@ -590,7 +590,8 @@ SUBROUTINE DDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL DGBT02( 'No transpose', N, N, KL, $ KU, NRHS, A, LDA, X, LDB, - $ WORK, LDB, RESULT( 2 ) ) + $ WORK, LDB, RWORK, + $ RESULT( 2 ) ) * * Check solution from generated exact * solution. @@ -707,6 +708,7 @@ SUBROUTINE DDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL DGBT02( TRANS, N, N, KL, KU, NRHS, $ ASAV, LDA, X, LDB, WORK, LDB, + $ RWORK( 2*NRHS+1 ), $ RESULT( 2 ) ) * * Check solution from generated exact @@ -892,7 +894,7 @@ SUBROUTINE DDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, CALL DLACPY( 'Full', N, NRHS, BSAV, LDB, WORK, $ LDB ) CALL DGBT02( TRANS, N, N, KL, KU, NRHS, ASAV, - $ LDA, X, LDB, WORK, LDB, + $ LDA, X, LDB, WORK, LDB, RWORK, $ RESULT( 2 ) ) * * Check solution from generated exact solution. diff --git a/TESTING/LIN/dgbt02.f b/TESTING/LIN/dgbt02.f index 7860e7f092..a5f669656f 100644 --- a/TESTING/LIN/dgbt02.f +++ b/TESTING/LIN/dgbt02.f @@ -9,7 +9,7 @@ * =========== * * SUBROUTINE DGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, -* LDB, RESID ) +* LDB, RWORK, RESID ) * * .. Scalar Arguments .. * CHARACTER TRANS @@ -17,7 +17,8 @@ * DOUBLE PRECISION RESID * .. * .. Array Arguments .. -* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * ) +* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * ), +* RWORK( * ) * .. * * @@ -27,9 +28,10 @@ *> \verbatim *> *> DGBT02 computes the residual for a solution of a banded system of -*> equations A*x = b or A'*x = b: -*> RESID = norm( B - A*X ) / ( norm(A) * norm(X) * EPS). -*> where EPS is the machine precision. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A or A**T, depending on TRANS, and EPS is the +*> machine precision. *> \endverbatim * * Arguments: @@ -39,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A'*x = b, where A' is the transpose of A -*> = 'C': A'*x = b, where A' is the transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose = Transpose) *> \endverbatim *> *> \param[in] M @@ -116,11 +118,18 @@ *> LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N). *> \endverbatim *> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)), +*> where LRWORK >= M when TRANS = 'T' or 'C'; otherwise, RWORK +*> is not referenced. +*> \endverbatim +* *> \param[out] RESID *> \verbatim *> RESID is DOUBLE PRECISION *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -135,7 +144,7 @@ * * ===================================================================== SUBROUTINE DGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, - $ LDB, RESID ) + $ LDB, RWORK, RESID ) * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -147,7 +156,8 @@ SUBROUTINE DGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, DOUBLE PRECISION RESID * .. * .. Array Arguments .. - DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * ) + DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * ), + $ RWORK( * ) * .. * * ===================================================================== @@ -169,7 +179,7 @@ SUBROUTINE DGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, EXTERNAL DGBMV * .. * .. Intrinsic Functions .. - INTRINSIC MAX, MIN + INTRINSIC ABS, MAX, MIN * .. * .. Executable Statements .. * @@ -183,14 +193,35 @@ SUBROUTINE DGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = DLAMCH( 'Epsilon' ) - KD = KU + 1 ANORM = ZERO - DO 10 J = 1, N - I1 = MAX( KD+1-J, 1 ) - I2 = MIN( KD+M-J, KL+KD ) - IF( I2.GE.I1 ) - $ ANORM = MAX( ANORM, DASUM( I2-I1+1, A( I1, J ), 1 ) ) - 10 CONTINUE + IF( LSAME( TRANS, 'N' ) ) THEN +* +* Find norm1(A). +* + KD = KU + 1 + DO 10 J = 1, N + I1 = MAX( KD+1-J, 1 ) + I2 = MIN( KD+M-J, KL+KD ) + IF( I2.GE.I1 ) + $ ANORM = MAX( ANORM, DASUM( I2-I1+1, A( I1, J ), 1 ) ) + 10 CONTINUE + ELSE +* +* Find normI(A). +* + DO 12 I1 = 1, M + RWORK( I1 ) = ZERO + 12 CONTINUE + DO 16 J = 1, N + KD = KU + 1 - J + DO 14 I1 = MAX( 1, J-KU ), MIN( M, J+KL ) + RWORK( I1 ) = RWORK( I1 ) + ABS( A( KD+I1, J ) ) + 14 CONTINUE + 16 CONTINUE + DO 18 I1 = 1, M + ANORM = MAX( ANORM, RWORK( I1 ) ) + 18 CONTINUE + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN @@ -202,7 +233,7 @@ SUBROUTINE DGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, N1 = M END IF * -* Compute B - A*X (or B - A'*X ) +* Compute B - op(A)*X * DO 20 J = 1, NRHS CALL DGBMV( TRANS, M, N, KL, KU, -ONE, A, LDA, X( 1, J ), 1, @@ -210,7 +241,7 @@ SUBROUTINE DGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, 20 CONTINUE * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). * RESID = ZERO DO 30 J = 1, NRHS diff --git a/TESTING/LIN/dget02.f b/TESTING/LIN/dget02.f index 02e54e7c1b..df42b0b6c5 100644 --- a/TESTING/LIN/dget02.f +++ b/TESTING/LIN/dget02.f @@ -28,9 +28,10 @@ *> \verbatim *> *> DGET02 computes the residual for a solution of a system of linear -*> equations A*x = b or A'*x = b: -*> RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ), -*> where EPS is the machine epsilon. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A or A**T, depending on TRANS, and EPS is the +*> machine epsilon. *> \endverbatim * * Arguments: @@ -40,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A'*x = b, where A' is the transpose of A -*> = 'C': A'*x = b, where A' is the transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose = Transpose) *> \endverbatim *> *> \param[in] M @@ -95,7 +96,7 @@ *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) *> On entry, the right hand side vectors for the system of *> linear equations. -*> On exit, B is overwritten with the difference B - A*X. +*> On exit, B is overwritten with the difference B - op(A)*X. *> \endverbatim *> *> \param[in] LDB @@ -114,7 +115,7 @@ *> \verbatim *> RESID is DOUBLE PRECISION *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -186,19 +187,23 @@ SUBROUTINE DGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = DLAMCH( 'Epsilon' ) - ANORM = DLANGE( '1', M, N, A, LDA, RWORK ) + IF( LSAME( TRANS, 'N' ) ) THEN + ANORM = DLANGE( '1', M, N, A, LDA, RWORK ) + ELSE + ANORM = DLANGE( 'I', M, N, A, LDA, RWORK ) + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN END IF * -* Compute B - A*X (or B - A'*X ) and store in B. +* Compute B - op(A)*X and store in B. * CALL DGEMM( TRANS, 'No transpose', N1, NRHS, N2, -ONE, A, LDA, X, $ LDX, ONE, B, LDB ) * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ) . * RESID = ZERO DO 10 J = 1, NRHS diff --git a/TESTING/LIN/schkgb.f b/TESTING/LIN/schkgb.f index 1418e1afbe..9ab8c0e018 100644 --- a/TESTING/LIN/schkgb.f +++ b/TESTING/LIN/schkgb.f @@ -563,7 +563,7 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, END IF * *+ TEST 2: -* Solve and compute residual for A * X = B. +* Solve and compute residual for op(A) * X = B. * SRNAMT = 'SLARHS' CALL SLARHS( PATH, XTYPE, ' ', TRANS, N, @@ -589,7 +589,7 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ WORK, LDB ) CALL SGBT02( TRANS, M, N, KL, KU, NRHS, A, $ LDA, X, LDB, WORK, LDB, - $ RESULT( 2 ) ) + $ RWORK, RESULT( 2 ) ) * *+ TEST 3: * Check solution from generated exact diff --git a/TESTING/LIN/sdrvgb.f b/TESTING/LIN/sdrvgb.f index 1dc318b154..d7c293ab6a 100644 --- a/TESTING/LIN/sdrvgb.f +++ b/TESTING/LIN/sdrvgb.f @@ -582,7 +582,8 @@ SUBROUTINE SDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL SGBT02( 'No transpose', N, N, KL, $ KU, NRHS, A, LDA, X, LDB, - $ WORK, LDB, RESULT( 2 ) ) + $ WORK, LDB, RWORK, + $ RESULT( 2 ) ) * * Check solution from generated exact * solution. @@ -699,6 +700,7 @@ SUBROUTINE SDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL SGBT02( TRANS, N, N, KL, KU, NRHS, $ ASAV, LDA, X, LDB, WORK, LDB, + $ RWORK( 2*NRHS+1 ), $ RESULT( 2 ) ) * * Check solution from generated exact diff --git a/TESTING/LIN/sdrvgbx.f b/TESTING/LIN/sdrvgbx.f index e8cb538fa9..342e25ac25 100644 --- a/TESTING/LIN/sdrvgbx.f +++ b/TESTING/LIN/sdrvgbx.f @@ -590,7 +590,8 @@ SUBROUTINE SDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL SGBT02( 'No transpose', N, N, KL, $ KU, NRHS, A, LDA, X, LDB, - $ WORK, LDB, RESULT( 2 ) ) + $ WORK, LDB, RWORK, + $ RESULT( 2 ) ) * * Check solution from generated exact * solution. @@ -707,6 +708,7 @@ SUBROUTINE SDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL SGBT02( TRANS, N, N, KL, KU, NRHS, $ ASAV, LDA, X, LDB, WORK, LDB, + $ RWORK( 2*NRHS+1 ), $ RESULT( 2 ) ) * * Check solution from generated exact @@ -893,7 +895,7 @@ SUBROUTINE SDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, CALL SLACPY( 'Full', N, NRHS, BSAV, LDB, WORK, $ LDB ) CALL SGBT02( TRANS, N, N, KL, KU, NRHS, ASAV, - $ LDA, X, LDB, WORK, LDB, + $ LDA, X, LDB, WORK, LDB, RWORK, $ RESULT( 2 ) ) * * Check solution from generated exact solution. diff --git a/TESTING/LIN/sgbt02.f b/TESTING/LIN/sgbt02.f index 14337e8032..b09b15bfd9 100644 --- a/TESTING/LIN/sgbt02.f +++ b/TESTING/LIN/sgbt02.f @@ -9,7 +9,7 @@ * =========== * * SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, -* LDB, RESID ) +* LDB, RWORK, RESID ) * * .. Scalar Arguments .. * CHARACTER TRANS @@ -17,7 +17,8 @@ * REAL RESID * .. * .. Array Arguments .. -* REAL A( LDA, * ), B( LDB, * ), X( LDX, * ) +* REAL A( LDA, * ), B( LDB, * ), X( LDX, * ), +* RWORK( * ) * .. * * @@ -27,9 +28,10 @@ *> \verbatim *> *> SGBT02 computes the residual for a solution of a banded system of -*> equations A*x = b or A'*x = b: -*> RESID = norm( B - A*X ) / ( norm(A) * norm(X) * EPS). -*> where EPS is the machine precision. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A or A**T, depending on TRANS, and EPS is the +*> machine precision. *> \endverbatim * * Arguments: @@ -39,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A'*x = b, where A' is the transpose of A -*> = 'C': A'*x = b, where A' is the transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose = Transpose) *> \endverbatim *> *> \param[in] M @@ -116,11 +118,18 @@ *> LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N). *> \endverbatim *> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (MAX(1,LRWORK)), +*> where LRWORK >= M when TRANS = 'T' or 'C'; otherwise, RWORK +*> is not referenced. +*> \endverbatim +* *> \param[out] RESID *> \verbatim *> RESID is REAL *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -135,7 +144,7 @@ * * ===================================================================== SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, - $ LDB, RESID ) + $ LDB, RWORK, RESID ) * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -147,7 +156,8 @@ SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, REAL RESID * .. * .. Array Arguments .. - REAL A( LDA, * ), B( LDB, * ), X( LDX, * ) + REAL A( LDA, * ), B( LDB, * ), X( LDX, * ), + $ RWORK( * ) * .. * * ===================================================================== @@ -169,7 +179,7 @@ SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, EXTERNAL SGBMV * .. * .. Intrinsic Functions .. - INTRINSIC MAX, MIN + INTRINSIC ABS, MAX, MIN * .. * .. Executable Statements .. * @@ -183,14 +193,35 @@ SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = SLAMCH( 'Epsilon' ) - KD = KU + 1 ANORM = ZERO - DO 10 J = 1, N - I1 = MAX( KD+1-J, 1 ) - I2 = MIN( KD+M-J, KL+KD ) - IF( I2.GE.I1 ) - $ ANORM = MAX( ANORM, SASUM( I2-I1+1, A( I1, J ), 1 ) ) - 10 CONTINUE + IF( LSAME( TRANS, 'N' ) ) THEN +* +* Find norm1(A). +* + KD = KU + 1 + DO 10 J = 1, N + I1 = MAX( KD+1-J, 1 ) + I2 = MIN( KD+M-J, KL+KD ) + IF( I2.GE.I1 ) + $ ANORM = MAX( ANORM, SASUM( I2-I1+1, A( I1, J ), 1 ) ) + 10 CONTINUE + ELSE +* +* Find normI(A). +* + DO 12 I1 = 1, M + RWORK( I1 ) = ZERO + 12 CONTINUE + DO 16 J = 1, N + KD = KU + 1 - J + DO 14 I1 = MAX( 1, J-KU ), MIN( M, J+KL ) + RWORK( I1 ) = RWORK( I1 ) + ABS( A( KD+I1, J ) ) + 14 CONTINUE + 16 CONTINUE + DO 18 I1 = 1, M + ANORM = MAX( ANORM, RWORK( I1 ) ) + 18 CONTINUE + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN @@ -202,7 +233,7 @@ SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, N1 = M END IF * -* Compute B - A*X (or B - A'*X ) +* Compute B - op(A)*X * DO 20 J = 1, NRHS CALL SGBMV( TRANS, M, N, KL, KU, -ONE, A, LDA, X( 1, J ), 1, @@ -210,7 +241,7 @@ SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, 20 CONTINUE * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). * RESID = ZERO DO 30 J = 1, NRHS diff --git a/TESTING/LIN/sget02.f b/TESTING/LIN/sget02.f index 817a9c049f..a04a2a1804 100644 --- a/TESTING/LIN/sget02.f +++ b/TESTING/LIN/sget02.f @@ -28,9 +28,10 @@ *> \verbatim *> *> SGET02 computes the residual for a solution of a system of linear -*> equations A*x = b or A'*x = b: -*> RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ), -*> where EPS is the machine epsilon. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A or A**T, depending on TRANS, and EPS is the +*> machine epsilon. *> \endverbatim * * Arguments: @@ -40,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A'*x = b, where A' is the transpose of A -*> = 'C': A'*x = b, where A' is the transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose = Transpose) *> \endverbatim *> *> \param[in] M @@ -95,7 +96,7 @@ *> B is REAL array, dimension (LDB,NRHS) *> On entry, the right hand side vectors for the system of *> linear equations. -*> On exit, B is overwritten with the difference B - A*X. +*> On exit, B is overwritten with the difference B - op(A)*X. *> \endverbatim *> *> \param[in] LDB @@ -114,7 +115,7 @@ *> \verbatim *> RESID is REAL *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -186,19 +187,23 @@ SUBROUTINE SGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = SLAMCH( 'Epsilon' ) - ANORM = SLANGE( '1', M, N, A, LDA, RWORK ) + IF( LSAME( TRANS, 'N' ) ) THEN + ANORM = SLANGE( '1', M, N, A, LDA, RWORK ) + ELSE + ANORM = SLANGE( 'I', M, N, A, LDA, RWORK ) + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN END IF * -* Compute B - A*X (or B - A'*X ) and store in B. +* Compute B - op(A)*X and store in B. * CALL SGEMM( TRANS, 'No transpose', N1, NRHS, N2, -ONE, A, LDA, X, $ LDX, ONE, B, LDB ) * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ) . * RESID = ZERO DO 10 J = 1, NRHS diff --git a/TESTING/LIN/zchkgb.f b/TESTING/LIN/zchkgb.f index 6e002bf187..d4d261d93a 100644 --- a/TESTING/LIN/zchkgb.f +++ b/TESTING/LIN/zchkgb.f @@ -563,7 +563,7 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, END IF * *+ TEST 2: -* Solve and compute residual for A * X = B. +* Solve and compute residual for op(A) * X = B. * SRNAMT = 'ZLARHS' CALL ZLARHS( PATH, XTYPE, ' ', TRANS, N, @@ -589,7 +589,7 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ WORK, LDB ) CALL ZGBT02( TRANS, M, N, KL, KU, NRHS, A, $ LDA, X, LDB, WORK, LDB, - $ RESULT( 2 ) ) + $ RWORK, RESULT( 2 ) ) * *+ TEST 3: * Check solution from generated exact diff --git a/TESTING/LIN/zdrvgb.f b/TESTING/LIN/zdrvgb.f index 5ee8cd8195..aa5712c3b5 100644 --- a/TESTING/LIN/zdrvgb.f +++ b/TESTING/LIN/zdrvgb.f @@ -582,7 +582,8 @@ SUBROUTINE ZDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL ZGBT02( 'No transpose', N, N, KL, $ KU, NRHS, A, LDA, X, LDB, - $ WORK, LDB, RESULT( 2 ) ) + $ WORK, LDB, RWORK, + $ RESULT( 2 ) ) * * Check solution from generated exact * solution. @@ -699,6 +700,7 @@ SUBROUTINE ZDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL ZGBT02( TRANS, N, N, KL, KU, NRHS, $ ASAV, LDA, X, LDB, WORK, LDB, + $ RWORK( 2*NRHS+1 ), $ RESULT( 2 ) ) * * Check solution from generated exact diff --git a/TESTING/LIN/zdrvgbx.f b/TESTING/LIN/zdrvgbx.f index 1a437daa69..9bf3560f20 100644 --- a/TESTING/LIN/zdrvgbx.f +++ b/TESTING/LIN/zdrvgbx.f @@ -590,7 +590,8 @@ SUBROUTINE ZDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL ZGBT02( 'No transpose', N, N, KL, $ KU, NRHS, A, LDA, X, LDB, - $ WORK, LDB, RESULT( 2 ) ) + $ WORK, LDB, RWORK, + $ RESULT( 2 ) ) * * Check solution from generated exact * solution. @@ -708,6 +709,7 @@ SUBROUTINE ZDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, $ WORK, LDB ) CALL ZGBT02( TRANS, N, N, KL, KU, NRHS, $ ASAV, LDA, X, LDB, WORK, LDB, + $ RWORK( 2*NRHS+1 ), $ RESULT( 2 ) ) * * Check solution from generated exact @@ -897,7 +899,8 @@ SUBROUTINE ZDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA, CALL ZLACPY( 'Full', N, NRHS, BSAV, LDB, WORK, $ LDB ) CALL ZGBT02( TRANS, N, N, KL, KU, NRHS, ASAV, - $ LDA, X, LDB, WORK, LDB, RESULT( 2 ) ) + $ LDA, X, LDB, WORK, LDB, RWORK, + $ RESULT( 2 ) ) * * Check solution from generated exact solution. * diff --git a/TESTING/LIN/zgbt02.f b/TESTING/LIN/zgbt02.f index a694e05665..21977ea64d 100644 --- a/TESTING/LIN/zgbt02.f +++ b/TESTING/LIN/zgbt02.f @@ -9,7 +9,7 @@ * =========== * * SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, -* LDB, RESID ) +* LDB, RWORK, RESID ) * * .. Scalar Arguments .. * CHARACTER TRANS @@ -17,6 +17,7 @@ * DOUBLE PRECISION RESID * .. * .. Array Arguments .. +* DOUBLE PRECISION RWORK( * ) * COMPLEX*16 A( LDA, * ), B( LDB, * ), X( LDX, * ) * .. * @@ -27,9 +28,10 @@ *> \verbatim *> *> ZGBT02 computes the residual for a solution of a banded system of -*> equations A*x = b or A'*x = b: -*> RESID = norm( B - A*X ) / ( norm(A) * norm(X) * EPS). -*> where EPS is the machine precision. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A, A**T, or A**H, depending on TRANS, and EPS is the +*> machine precision. *> \endverbatim * * Arguments: @@ -39,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A'*x = b, where A' is the transpose of A -*> = 'C': A'*x = b, where A' is the transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose) *> \endverbatim *> *> \param[in] M @@ -116,11 +118,18 @@ *> LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N). *> \endverbatim *> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)), +*> where LRWORK >= M when TRANS = 'T' or 'C'; otherwise, RWORK +*> is not referenced. +*> \endverbatim +* *> \param[out] RESID *> \verbatim *> RESID is DOUBLE PRECISION *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -135,7 +144,7 @@ * * ===================================================================== SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, - $ LDB, RESID ) + $ LDB, RWORK, RESID ) * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -147,6 +156,7 @@ SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, DOUBLE PRECISION RESID * .. * .. Array Arguments .. + DOUBLE PRECISION RWORK( * ) COMPLEX*16 A( LDA, * ), B( LDB, * ), X( LDX, * ) * .. * @@ -161,6 +171,7 @@ SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * .. Local Scalars .. INTEGER I1, I2, J, KD, N1 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM + COMPLEX*16 ZDUM * .. * .. External Functions .. LOGICAL LSAME @@ -170,8 +181,14 @@ SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * .. External Subroutines .. EXTERNAL ZGBMV * .. +* .. Statement Functions .. + DOUBLE PRECISION CABS1 +* .. * .. Intrinsic Functions .. - INTRINSIC MAX, MIN + INTRINSIC ABS, DBLE, DIMAG, MAX, MIN +* .. +* .. Statement Function definitions .. + CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) * .. * .. Executable Statements .. * @@ -185,14 +202,35 @@ SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = DLAMCH( 'Epsilon' ) - KD = KU + 1 ANORM = ZERO - DO 10 J = 1, N - I1 = MAX( KD+1-J, 1 ) - I2 = MIN( KD+M-J, KL+KD ) - IF( I2.GE.I1 ) - $ ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) ) - 10 CONTINUE + IF( LSAME( TRANS, 'N' ) ) THEN +* +* Find norm1(A). +* + KD = KU + 1 + DO 10 J = 1, N + I1 = MAX( KD+1-J, 1 ) + I2 = MIN( KD+M-J, KL+KD ) + IF( I2.GE.I1 ) + $ ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) ) + 10 CONTINUE + ELSE +* +* Find normI(A). +* + DO 12 I1 = 1, M + RWORK( I1 ) = ZERO + 12 CONTINUE + DO 16 J = 1, N + KD = KU + 1 - J + DO 14 I1 = MAX( 1, J-KU ), MIN( M, J+KL ) + RWORK( I1 ) = RWORK( I1 ) + CABS1( A( KD+I1, J ) ) + 14 CONTINUE + 16 CONTINUE + DO 18 I1 = 1, M + ANORM = MAX( ANORM, RWORK( I1 ) ) + 18 CONTINUE + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN @@ -204,7 +242,7 @@ SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, N1 = M END IF * -* Compute B - A*X (or B - A'*X ) +* Compute B - op(A)*X * DO 20 J = 1, NRHS CALL ZGBMV( TRANS, M, N, KL, KU, -CONE, A, LDA, X( 1, J ), 1, @@ -212,7 +250,7 @@ SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, 20 CONTINUE * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). * RESID = ZERO DO 30 J = 1, NRHS diff --git a/TESTING/LIN/zget02.f b/TESTING/LIN/zget02.f index 16e9ee0d92..1ffa387959 100644 --- a/TESTING/LIN/zget02.f +++ b/TESTING/LIN/zget02.f @@ -28,9 +28,10 @@ *> \verbatim *> *> ZGET02 computes the residual for a solution of a system of linear -*> equations A*x = b or A'*x = b: -*> RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ), -*> where EPS is the machine epsilon. +*> equations op(A)*X = B: +*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), +*> where op(A) = A, A**T, or A**H, depending on TRANS, and EPS is the +*> machine epsilon. *> \endverbatim * * Arguments: @@ -40,9 +41,9 @@ *> \verbatim *> TRANS is CHARACTER*1 *> Specifies the form of the system of equations: -*> = 'N': A *x = b -*> = 'T': A^T*x = b, where A^T is the transpose of A -*> = 'C': A^H*x = b, where A^H is the conjugate transpose of A +*> = 'N': A * X = B (No transpose) +*> = 'T': A**T * X = B (Transpose) +*> = 'C': A**H * X = B (Conjugate transpose) *> \endverbatim *> *> \param[in] M @@ -95,7 +96,7 @@ *> B is COMPLEX*16 array, dimension (LDB,NRHS) *> On entry, the right hand side vectors for the system of *> linear equations. -*> On exit, B is overwritten with the difference B - A*X. +*> On exit, B is overwritten with the difference B - op(A)*X. *> \endverbatim *> *> \param[in] LDB @@ -114,7 +115,7 @@ *> \verbatim *> RESID is DOUBLE PRECISION *> The maximum over the number of right hand sides of -*> norm(B - A*X) / ( norm(A) * norm(X) * EPS ). +*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). *> \endverbatim * * Authors: @@ -188,19 +189,23 @@ SUBROUTINE ZGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, * Exit with RESID = 1/EPS if ANORM = 0. * EPS = DLAMCH( 'Epsilon' ) - ANORM = ZLANGE( '1', M, N, A, LDA, RWORK ) + IF( LSAME( TRANS, 'N' ) ) THEN + ANORM = ZLANGE( '1', M, N, A, LDA, RWORK ) + ELSE + ANORM = ZLANGE( 'I', M, N, A, LDA, RWORK ) + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN END IF * -* Compute B - A*X (or B - A'*X ) and store in B. +* Compute B - op(A)*X and store in B. * CALL ZGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X, $ LDX, CONE, B, LDB ) * * Compute the maximum over the number of right hand sides of -* norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . +* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ) . * RESID = ZERO DO 10 J = 1, NRHS From 036ec87263ae6f300b43ec02ce7d0dc7037ea364 Mon Sep 17 00:00:00 2001 From: Igor Zhuravlov Date: Sun, 30 May 2021 19:35:35 +1000 Subject: [PATCH 4/4] code fixed: a check was added in xGBT02 to handle NaN while norm computing A check was borrowed from xLANGB. --- TESTING/LIN/cgbt02.f | 15 +++++++++------ TESTING/LIN/dgbt02.f | 15 +++++++++------ TESTING/LIN/sgbt02.f | 15 +++++++++------ TESTING/LIN/zgbt02.f | 15 +++++++++------ 4 files changed, 36 insertions(+), 24 deletions(-) diff --git a/TESTING/LIN/cgbt02.f b/TESTING/LIN/cgbt02.f index 12fc57284b..8c51e8b2a2 100644 --- a/TESTING/LIN/cgbt02.f +++ b/TESTING/LIN/cgbt02.f @@ -170,13 +170,13 @@ SUBROUTINE CGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * .. * .. Local Scalars .. INTEGER I1, I2, J, KD, N1 - REAL ANORM, BNORM, EPS, XNORM + REAL ANORM, BNORM, EPS, TEMP, XNORM COMPLEX CDUM * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, SISNAN REAL SCASUM, SLAMCH - EXTERNAL LSAME, SCASUM, SLAMCH + EXTERNAL LSAME, SCASUM, SISNAN, SLAMCH * .. * .. External Subroutines .. EXTERNAL CGBMV @@ -211,8 +211,10 @@ SUBROUTINE CGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, DO 10 J = 1, N I1 = MAX( KD+1-J, 1 ) I2 = MIN( KD+M-J, KL+KD ) - IF( I2.GE.I1 ) - $ ANORM = MAX( ANORM, SCASUM( I2-I1+1, A( I1, J ), 1 ) ) + IF( I2.GE.I1 ) THEN + TEMP = SCASUM( I2-I1+1, A( I1, J ), 1 ) + IF( ANORM.LT.TEMP .OR. SISNAN( TEMP ) ) ANORM = TEMP + END IF 10 CONTINUE ELSE * @@ -228,7 +230,8 @@ SUBROUTINE CGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, 14 CONTINUE 16 CONTINUE DO 18 I1 = 1, M - ANORM = MAX( ANORM, RWORK( I1 ) ) + TEMP = RWORK( I1 ) + IF( ANORM.LT.TEMP .OR. SISNAN( TEMP ) ) ANORM = TEMP 18 CONTINUE END IF IF( ANORM.LE.ZERO ) THEN diff --git a/TESTING/LIN/dgbt02.f b/TESTING/LIN/dgbt02.f index a5f669656f..ef7e00c3d6 100644 --- a/TESTING/LIN/dgbt02.f +++ b/TESTING/LIN/dgbt02.f @@ -168,12 +168,12 @@ SUBROUTINE DGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * .. * .. Local Scalars .. INTEGER I1, I2, J, KD, N1 - DOUBLE PRECISION ANORM, BNORM, EPS, XNORM + DOUBLE PRECISION ANORM, BNORM, EPS, TEMP, XNORM * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL DISNAN, LSAME DOUBLE PRECISION DASUM, DLAMCH - EXTERNAL LSAME, DASUM, DLAMCH + EXTERNAL DASUM, DISNAN, DLAMCH, LSAME * .. * .. External Subroutines .. EXTERNAL DGBMV @@ -202,8 +202,10 @@ SUBROUTINE DGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, DO 10 J = 1, N I1 = MAX( KD+1-J, 1 ) I2 = MIN( KD+M-J, KL+KD ) - IF( I2.GE.I1 ) - $ ANORM = MAX( ANORM, DASUM( I2-I1+1, A( I1, J ), 1 ) ) + IF( I2.GE.I1 ) THEN + TEMP = DASUM( I2-I1+1, A( I1, J ), 1 ) + IF( ANORM.LT.TEMP .OR. DISNAN( TEMP ) ) ANORM = TEMP + END IF 10 CONTINUE ELSE * @@ -219,7 +221,8 @@ SUBROUTINE DGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, 14 CONTINUE 16 CONTINUE DO 18 I1 = 1, M - ANORM = MAX( ANORM, RWORK( I1 ) ) + TEMP = RWORK( I1 ) + IF( ANORM.LT.TEMP .OR. DISNAN( TEMP ) ) ANORM = TEMP 18 CONTINUE END IF IF( ANORM.LE.ZERO ) THEN diff --git a/TESTING/LIN/sgbt02.f b/TESTING/LIN/sgbt02.f index b09b15bfd9..2c0734d150 100644 --- a/TESTING/LIN/sgbt02.f +++ b/TESTING/LIN/sgbt02.f @@ -168,12 +168,12 @@ SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * .. * .. Local Scalars .. INTEGER I1, I2, J, KD, N1 - REAL ANORM, BNORM, EPS, XNORM + REAL ANORM, BNORM, EPS, TEMP, XNORM * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, SISNAN REAL SASUM, SLAMCH - EXTERNAL LSAME, SASUM, SLAMCH + EXTERNAL LSAME, SASUM, SISNAN, SLAMCH * .. * .. External Subroutines .. EXTERNAL SGBMV @@ -202,8 +202,10 @@ SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, DO 10 J = 1, N I1 = MAX( KD+1-J, 1 ) I2 = MIN( KD+M-J, KL+KD ) - IF( I2.GE.I1 ) - $ ANORM = MAX( ANORM, SASUM( I2-I1+1, A( I1, J ), 1 ) ) + IF( I2.GE.I1 ) THEN + TEMP = SASUM( I2-I1+1, A( I1, J ), 1 ) + IF( ANORM.LT.TEMP .OR. SISNAN( TEMP ) ) ANORM = TEMP + END IF 10 CONTINUE ELSE * @@ -219,7 +221,8 @@ SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, 14 CONTINUE 16 CONTINUE DO 18 I1 = 1, M - ANORM = MAX( ANORM, RWORK( I1 ) ) + TEMP = RWORK( I1 ) + IF( ANORM.LT.TEMP .OR. SISNAN( TEMP ) ) ANORM = TEMP 18 CONTINUE END IF IF( ANORM.LE.ZERO ) THEN diff --git a/TESTING/LIN/zgbt02.f b/TESTING/LIN/zgbt02.f index 21977ea64d..3547ea0b5e 100644 --- a/TESTING/LIN/zgbt02.f +++ b/TESTING/LIN/zgbt02.f @@ -170,13 +170,13 @@ SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, * .. * .. Local Scalars .. INTEGER I1, I2, J, KD, N1 - DOUBLE PRECISION ANORM, BNORM, EPS, XNORM + DOUBLE PRECISION ANORM, BNORM, EPS, TEMP, XNORM COMPLEX*16 ZDUM * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL DISNAN, LSAME DOUBLE PRECISION DLAMCH, DZASUM - EXTERNAL LSAME, DLAMCH, DZASUM + EXTERNAL DISNAN, DLAMCH, DZASUM, LSAME * .. * .. External Subroutines .. EXTERNAL ZGBMV @@ -211,8 +211,10 @@ SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, DO 10 J = 1, N I1 = MAX( KD+1-J, 1 ) I2 = MIN( KD+M-J, KL+KD ) - IF( I2.GE.I1 ) - $ ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) ) + IF( I2.GE.I1 ) THEN + TEMP = DZASUM( I2-I1+1, A( I1, J ), 1 ) + IF( ANORM.LT.TEMP .OR. DISNAN( TEMP ) ) ANORM = TEMP + END IF 10 CONTINUE ELSE * @@ -228,7 +230,8 @@ SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, 14 CONTINUE 16 CONTINUE DO 18 I1 = 1, M - ANORM = MAX( ANORM, RWORK( I1 ) ) + TEMP = RWORK( I1 ) + IF( ANORM.LT.TEMP .OR. DISNAN( TEMP ) ) ANORM = TEMP 18 CONTINUE END IF IF( ANORM.LE.ZERO ) THEN