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 433b45f4ec..9bb6f29b79 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 @@ -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 b0a8be002f..389725eaed 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 @@ -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 747cb40adb..68ff66e07e 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 @@ -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 9576c47ce4..8c51e8b2a2 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, * ) * .. * @@ -160,18 +170,25 @@ 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 * .. +* .. 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,13 +202,38 @@ 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 ) - 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 ) THEN + TEMP = SCASUM( I2-I1+1, A( I1, J ), 1 ) + IF( ANORM.LT.TEMP .OR. SISNAN( TEMP ) ) ANORM = TEMP + END IF + 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 + TEMP = RWORK( I1 ) + IF( ANORM.LT.TEMP .OR. SISNAN( TEMP ) ) ANORM = TEMP + 18 CONTINUE + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN @@ -203,7 +245,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, @@ -211,7 +253,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 c8a7efeb76..fa8dce5cff 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 @@ -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 ed0e18a850..ad323d5ec0 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 @@ -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 b8236703ad..8a9fdb645a 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 @@ -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 92bf550d23..ef7e00c3d6 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( * ) * .. * * ===================================================================== @@ -158,18 +168,18 @@ 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 * .. * .. Intrinsic Functions .. - INTRINSIC MAX, MIN + INTRINSIC ABS, MAX, MIN * .. * .. Executable Statements .. * @@ -183,13 +193,38 @@ 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 ) - 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 ) THEN + TEMP = DASUM( I2-I1+1, A( I1, J ), 1 ) + IF( ANORM.LT.TEMP .OR. DISNAN( TEMP ) ) ANORM = TEMP + END IF + 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 + TEMP = RWORK( I1 ) + IF( ANORM.LT.TEMP .OR. DISNAN( TEMP ) ) ANORM = TEMP + 18 CONTINUE + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN @@ -201,7 +236,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, @@ -209,7 +244,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 7eee800aa4..9ab8c0e018 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 @@ -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 8c65b6ae9b..d7c293ab6a 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 @@ -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 a38c3486a4..342e25ac25 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 @@ -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 98725ca67b..2c0734d150 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( * ) * .. * * ===================================================================== @@ -158,18 +168,18 @@ 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 * .. * .. Intrinsic Functions .. - INTRINSIC MAX, MIN + INTRINSIC ABS, MAX, MIN * .. * .. Executable Statements .. * @@ -183,13 +193,38 @@ 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 ) - 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 ) THEN + TEMP = SASUM( I2-I1+1, A( I1, J ), 1 ) + IF( ANORM.LT.TEMP .OR. SISNAN( TEMP ) ) ANORM = TEMP + END IF + 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 + TEMP = RWORK( I1 ) + IF( ANORM.LT.TEMP .OR. SISNAN( TEMP ) ) ANORM = TEMP + 18 CONTINUE + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN @@ -201,7 +236,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, @@ -209,7 +244,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 39de25197c..d4d261d93a 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 @@ -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 3431b4b31d..aa5712c3b5 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 @@ -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 cd0e077549..9bf3560f20 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 @@ -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 250efea3d2..3547ea0b5e 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, * ) * .. * @@ -160,18 +170,25 @@ 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 * .. +* .. 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,13 +202,38 @@ 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 ) - 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 ) THEN + TEMP = DZASUM( I2-I1+1, A( I1, J ), 1 ) + IF( ANORM.LT.TEMP .OR. DISNAN( TEMP ) ) ANORM = TEMP + END IF + 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 + TEMP = RWORK( I1 ) + IF( ANORM.LT.TEMP .OR. DISNAN( TEMP ) ) ANORM = TEMP + 18 CONTINUE + END IF IF( ANORM.LE.ZERO ) THEN RESID = ONE / EPS RETURN @@ -203,7 +245,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, @@ -211,7 +253,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