1*> \brief \b SGBT02 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, 12* LDB, RWORK, RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER TRANS 16* INTEGER KL, KU, LDA, LDB, LDX, M, N, NRHS 17* REAL RESID 18* .. 19* .. Array Arguments .. 20* REAL A( LDA, * ), B( LDB, * ), X( LDX, * ), 21* RWORK( * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> SGBT02 computes the residual for a solution of a banded system of 31*> equations op(A)*X = B: 32*> RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ), 33*> where op(A) = A or A**T, depending on TRANS, and EPS is the 34*> machine epsilon. 35*> The norm used is the 1-norm. 36*> \endverbatim 37* 38* Arguments: 39* ========== 40* 41*> \param[in] TRANS 42*> \verbatim 43*> TRANS is CHARACTER*1 44*> Specifies the form of the system of equations: 45*> = 'N': A * X = B (No transpose) 46*> = 'T': A**T * X = B (Transpose) 47*> = 'C': A**H * X = B (Conjugate transpose = Transpose) 48*> \endverbatim 49*> 50*> \param[in] M 51*> \verbatim 52*> M is INTEGER 53*> The number of rows of the matrix A. M >= 0. 54*> \endverbatim 55*> 56*> \param[in] N 57*> \verbatim 58*> N is INTEGER 59*> The number of columns of the matrix A. N >= 0. 60*> \endverbatim 61*> 62*> \param[in] KL 63*> \verbatim 64*> KL is INTEGER 65*> The number of subdiagonals within the band of A. KL >= 0. 66*> \endverbatim 67*> 68*> \param[in] KU 69*> \verbatim 70*> KU is INTEGER 71*> The number of superdiagonals within the band of A. KU >= 0. 72*> \endverbatim 73*> 74*> \param[in] NRHS 75*> \verbatim 76*> NRHS is INTEGER 77*> The number of columns of B. NRHS >= 0. 78*> \endverbatim 79*> 80*> \param[in] A 81*> \verbatim 82*> A is REAL array, dimension (LDA,N) 83*> The original matrix A in band storage, stored in rows 1 to 84*> KL+KU+1. 85*> \endverbatim 86*> 87*> \param[in] LDA 88*> \verbatim 89*> LDA is INTEGER 90*> The leading dimension of the array A. LDA >= max(1,KL+KU+1). 91*> \endverbatim 92*> 93*> \param[in] X 94*> \verbatim 95*> X is REAL array, dimension (LDX,NRHS) 96*> The computed solution vectors for the system of linear 97*> equations. 98*> \endverbatim 99*> 100*> \param[in] LDX 101*> \verbatim 102*> LDX is INTEGER 103*> The leading dimension of the array X. If TRANS = 'N', 104*> LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M). 105*> \endverbatim 106*> 107*> \param[in,out] B 108*> \verbatim 109*> B is REAL array, dimension (LDB,NRHS) 110*> On entry, the right hand side vectors for the system of 111*> linear equations. 112*> On exit, B is overwritten with the difference B - A*X. 113*> \endverbatim 114*> 115*> \param[in] LDB 116*> \verbatim 117*> LDB is INTEGER 118*> The leading dimension of the array B. IF TRANS = 'N', 119*> LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N). 120*> \endverbatim 121*> 122*> \param[out] RWORK 123*> \verbatim 124*> RWORK is REAL array, dimension (MAX(1,LRWORK)), 125*> where LRWORK >= M when TRANS = 'T' or 'C'; otherwise, RWORK 126*> is not referenced. 127*> \endverbatim 128* 129*> \param[out] RESID 130*> \verbatim 131*> RESID is REAL 132*> The maximum over the number of right hand sides of 133*> norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). 134*> \endverbatim 135* 136* Authors: 137* ======== 138* 139*> \author Univ. of Tennessee 140*> \author Univ. of California Berkeley 141*> \author Univ. of Colorado Denver 142*> \author NAG Ltd. 143* 144*> \ingroup single_lin 145* 146* ===================================================================== 147 SUBROUTINE SGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, 148 $ LDB, RWORK, RESID ) 149* 150* -- LAPACK test routine -- 151* -- LAPACK is a software package provided by Univ. of Tennessee, -- 152* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 153* 154* .. Scalar Arguments .. 155 CHARACTER TRANS 156 INTEGER KL, KU, LDA, LDB, LDX, M, N, NRHS 157 REAL RESID 158* .. 159* .. Array Arguments .. 160 REAL A( LDA, * ), B( LDB, * ), X( LDX, * ), 161 $ RWORK( * ) 162* .. 163* 164* ===================================================================== 165* 166* .. Parameters .. 167 REAL ZERO, ONE 168 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 169* .. 170* .. Local Scalars .. 171 INTEGER I1, I2, J, KD, N1 172 REAL ANORM, BNORM, EPS, TEMP, XNORM 173* .. 174* .. External Functions .. 175 LOGICAL LSAME, SISNAN 176 REAL SASUM, SLAMCH 177 EXTERNAL LSAME, SASUM, SISNAN, SLAMCH 178* .. 179* .. External Subroutines .. 180 EXTERNAL SGBMV 181* .. 182* .. Intrinsic Functions .. 183 INTRINSIC ABS, MAX, MIN 184* .. 185* .. Executable Statements .. 186* 187* Quick return if N = 0 pr NRHS = 0 188* 189 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.LE.0 ) THEN 190 RESID = ZERO 191 RETURN 192 END IF 193* 194* Exit with RESID = 1/EPS if ANORM = 0. 195* 196 EPS = SLAMCH( 'Epsilon' ) 197 ANORM = ZERO 198 IF( LSAME( TRANS, 'N' ) ) THEN 199* 200* Find norm1(A). 201* 202 KD = KU + 1 203 DO 10 J = 1, N 204 I1 = MAX( KD+1-J, 1 ) 205 I2 = MIN( KD+M-J, KL+KD ) 206 IF( I2.GE.I1 ) THEN 207 TEMP = SASUM( I2-I1+1, A( I1, J ), 1 ) 208 IF( ANORM.LT.TEMP .OR. SISNAN( TEMP ) ) ANORM = TEMP 209 END IF 210 10 CONTINUE 211 ELSE 212* 213* Find normI(A). 214* 215 DO 12 I1 = 1, M 216 RWORK( I1 ) = ZERO 217 12 CONTINUE 218 DO 16 J = 1, N 219 KD = KU + 1 - J 220 DO 14 I1 = MAX( 1, J-KU ), MIN( M, J+KL ) 221 RWORK( I1 ) = RWORK( I1 ) + ABS( A( KD+I1, J ) ) 222 14 CONTINUE 223 16 CONTINUE 224 DO 18 I1 = 1, M 225 TEMP = RWORK( I1 ) 226 IF( ANORM.LT.TEMP .OR. SISNAN( TEMP ) ) ANORM = TEMP 227 18 CONTINUE 228 END IF 229 IF( ANORM.LE.ZERO ) THEN 230 RESID = ONE / EPS 231 RETURN 232 END IF 233* 234 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN 235 N1 = N 236 ELSE 237 N1 = M 238 END IF 239* 240* Compute B - op(A)*X 241* 242 DO 20 J = 1, NRHS 243 CALL SGBMV( TRANS, M, N, KL, KU, -ONE, A, LDA, X( 1, J ), 1, 244 $ ONE, B( 1, J ), 1 ) 245 20 CONTINUE 246* 247* Compute the maximum over the number of right hand sides of 248* norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ). 249* 250 RESID = ZERO 251 DO 30 J = 1, NRHS 252 BNORM = SASUM( N1, B( 1, J ), 1 ) 253 XNORM = SASUM( N1, X( 1, J ), 1 ) 254 IF( XNORM.LE.ZERO ) THEN 255 RESID = ONE / EPS 256 ELSE 257 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) 258 END IF 259 30 CONTINUE 260* 261 RETURN 262* 263* End of SGBT02 264* 265 END 266