1       SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B,
  2      $                   LDB, RESID )
  3 *
  4 *  -- LAPACK test routine (version 3.1) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     November 2006
  7 *
  8 *     .. Scalar Arguments ..
  9       CHARACTER          TRANS
 10       INTEGER            KL, KU, LDA, LDB, LDX, M, N, NRHS
 11       REAL               RESID
 12 *     ..
 13 *     .. Array Arguments ..
 14       REAL               A( LDA, * ), B( LDB, * ), X( LDX, * )
 15 *     ..
 16 *
 17 *  Purpose
 18 *  =======
 19 *
 20 *  SGBT02 computes the residual for a solution of a banded system of
 21 *  equations  A*x = b  or  A'*x = b:
 22 *     RESID = norm( B - A*X ) / ( norm(A) * norm(X) * EPS).
 23 *  where EPS is the machine precision.
 24 *
 25 *  Arguments
 26 *  =========
 27 *
 28 *  TRANS   (input) CHARACTER*1
 29 *          Specifies the form of the system of equations:
 30 *          = 'N':  A *x = b
 31 *          = 'T':  A'*x = b, where A' is the transpose of A
 32 *          = 'C':  A'*x = b, where A' is the transpose of A
 33 *
 34 *  M       (input) INTEGER
 35 *          The number of rows of the matrix A.  M >= 0.
 36 *
 37 *  N       (input) INTEGER
 38 *          The number of columns of the matrix A.  N >= 0.
 39 *
 40 *  KL      (input) INTEGER
 41 *          The number of subdiagonals within the band of A.  KL >= 0.
 42 *
 43 *  KU      (input) INTEGER
 44 *          The number of superdiagonals within the band of A.  KU >= 0.
 45 *
 46 *  NRHS    (input) INTEGER
 47 *          The number of columns of B.  NRHS >= 0.
 48 *
 49 *  A       (input) REAL array, dimension (LDA,N)
 50 *          The original matrix A in band storage, stored in rows 1 to
 51 *          KL+KU+1.
 52 *
 53 *  LDA     (input) INTEGER
 54 *          The leading dimension of the array A.  LDA >= max(1,KL+KU+1).
 55 *
 56 *  X       (input) REAL array, dimension (LDX,NRHS)
 57 *          The computed solution vectors for the system of linear
 58 *          equations.
 59 *
 60 *  LDX     (input) INTEGER
 61 *          The leading dimension of the array X.  If TRANS = 'N',
 62 *          LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M).
 63 *
 64 *  B       (input/output) REAL array, dimension (LDB,NRHS)
 65 *          On entry, the right hand side vectors for the system of
 66 *          linear equations.
 67 *          On exit, B is overwritten with the difference B - A*X.
 68 *
 69 *  LDB     (input) INTEGER
 70 *          The leading dimension of the array B.  IF TRANS = 'N',
 71 *          LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N).
 72 *
 73 *  RESID   (output) REAL
 74 *          The maximum over the number of right hand sides of
 75 *          norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
 76 *
 77 *  =====================================================================
 78 *
 79 *     .. Parameters ..
 80       REAL               ZERO, ONE
 81       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 82 *     ..
 83 *     .. Local Scalars ..
 84       INTEGER            I1, I2, J, KD, N1
 85       REAL               ANORM, BNORM, EPS, XNORM
 86 *     ..
 87 *     .. External Functions ..
 88       LOGICAL            LSAME
 89       REAL               SASUM, SLAMCH
 90       EXTERNAL           LSAME, SASUM, SLAMCH
 91 *     ..
 92 *     .. External Subroutines ..
 93       EXTERNAL           SGBMV
 94 *     ..
 95 *     .. Intrinsic Functions ..
 96       INTRINSIC          MAXMIN
 97 *     ..
 98 *     .. Executable Statements ..
 99 *
100 *     Quick return if N = 0 pr NRHS = 0
101 *
102       IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.LE.0 ) THEN
103          RESID = ZERO
104          RETURN
105       END IF
106 *
107 *     Exit with RESID = 1/EPS if ANORM = 0.
108 *
109       EPS = SLAMCH( 'Epsilon' )
110       KD = KU + 1
111       ANORM = ZERO
112       DO 10 J = 1, N
113          I1 = MAX( KD+1-J, 1 )
114          I2 = MIN( KD+M-J, KL+KD )
115          ANORM = MAX( ANORM, SASUM( I2-I1+1, A( I1, J ), 1 ) )
116    10 CONTINUE
117       IF( ANORM.LE.ZERO ) THEN
118          RESID = ONE / EPS
119          RETURN
120       END IF
121 *
122       IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
123          N1 = N
124       ELSE
125          N1 = M
126       END IF
127 *
128 *     Compute  B - A*X (or  B - A'*X )
129 *
130       DO 20 J = 1, NRHS
131          CALL SGBMV( TRANS, M, N, KL, KU, -ONE, A, LDA, X( 1, J ), 1,
132      $               ONE, B( 1, J ), 1 )
133    20 CONTINUE
134 *
135 *     Compute the maximum over the number of right hand sides of
136 *        norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
137 *
138       RESID = ZERO
139       DO 30 J = 1, NRHS
140          BNORM = SASUM( N1, B( 1, J ), 1 )
141          XNORM = SASUM( N1, X( 1, J ), 1 )
142          IF( XNORM.LE.ZERO ) THEN
143             RESID = ONE / EPS
144          ELSE
145             RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
146          END IF
147    30 CONTINUE
148 *
149       RETURN
150 *
151 *     End of SGBT02
152 *
153       END