Merge pull request #571 from jip/fix

modify xGBT02 and xGET02 to use 1-norm of op(A) not A (no.2 from PR #562)

Author: langou

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

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

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

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

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

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

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.
 *