1*> \brief \b DTBT02 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 DTBT02( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, X, 12* LDX, B, LDB, WORK, RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER DIAG, TRANS, UPLO 16* INTEGER KD, LDAB, LDB, LDX, N, NRHS 17* DOUBLE PRECISION RESID 18* .. 19* .. Array Arguments .. 20* DOUBLE PRECISION AB( LDAB, * ), B( LDB, * ), WORK( * ), 21* $ X( LDX, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> DTBT02 computes the residual for the computed solution to a 31*> triangular system of linear equations A*x = b or A' *x = b when 32*> A is a triangular band matrix. Here A' is the transpose of A and 33*> x and b are N by NRHS matrices. The test ratio is the maximum over 34*> the number of right hand sides of 35*> norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ), 36*> where op(A) denotes A or A' and EPS is the machine epsilon. 37*> \endverbatim 38* 39* Arguments: 40* ========== 41* 42*> \param[in] UPLO 43*> \verbatim 44*> UPLO is CHARACTER*1 45*> Specifies whether the matrix A is upper or lower triangular. 46*> = 'U': Upper triangular 47*> = 'L': Lower triangular 48*> \endverbatim 49*> 50*> \param[in] TRANS 51*> \verbatim 52*> TRANS is CHARACTER*1 53*> Specifies the operation applied to A. 54*> = 'N': A *x = b (No transpose) 55*> = 'T': A'*x = b (Transpose) 56*> = 'C': A'*x = b (Conjugate transpose = Transpose) 57*> \endverbatim 58*> 59*> \param[in] DIAG 60*> \verbatim 61*> DIAG is CHARACTER*1 62*> Specifies whether or not the matrix A is unit triangular. 63*> = 'N': Non-unit triangular 64*> = 'U': Unit triangular 65*> \endverbatim 66*> 67*> \param[in] N 68*> \verbatim 69*> N is INTEGER 70*> The order of the matrix A. N >= 0. 71*> \endverbatim 72*> 73*> \param[in] KD 74*> \verbatim 75*> KD is INTEGER 76*> The number of superdiagonals or subdiagonals of the 77*> triangular band matrix A. KD >= 0. 78*> \endverbatim 79*> 80*> \param[in] NRHS 81*> \verbatim 82*> NRHS is INTEGER 83*> The number of right hand sides, i.e., the number of columns 84*> of the matrices X and B. NRHS >= 0. 85*> \endverbatim 86*> 87*> \param[in] AB 88*> \verbatim 89*> AB is DOUBLE PRECISION array, dimension (LDAB,N) 90*> The upper or lower triangular band matrix A, stored in the 91*> first kd+1 rows of the array. The j-th column of A is stored 92*> in the j-th column of the array AB as follows: 93*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 94*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 95*> \endverbatim 96*> 97*> \param[in] LDAB 98*> \verbatim 99*> LDAB is INTEGER 100*> The leading dimension of the array AB. LDAB >= KD+1. 101*> \endverbatim 102*> 103*> \param[in] X 104*> \verbatim 105*> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 106*> The computed solution vectors for the system of linear 107*> equations. 108*> \endverbatim 109*> 110*> \param[in] LDX 111*> \verbatim 112*> LDX is INTEGER 113*> The leading dimension of the array X. LDX >= max(1,N). 114*> \endverbatim 115*> 116*> \param[in] B 117*> \verbatim 118*> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 119*> The right hand side vectors for the system of linear 120*> equations. 121*> \endverbatim 122*> 123*> \param[in] LDB 124*> \verbatim 125*> LDB is INTEGER 126*> The leading dimension of the array B. LDB >= max(1,N). 127*> \endverbatim 128*> 129*> \param[out] WORK 130*> \verbatim 131*> WORK is DOUBLE PRECISION array, dimension (N) 132*> \endverbatim 133*> 134*> \param[out] RESID 135*> \verbatim 136*> RESID is DOUBLE PRECISION 137*> The maximum over the number of right hand sides of 138*> norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ). 139*> \endverbatim 140* 141* Authors: 142* ======== 143* 144*> \author Univ. of Tennessee 145*> \author Univ. of California Berkeley 146*> \author Univ. of Colorado Denver 147*> \author NAG Ltd. 148* 149*> \date November 2011 150* 151*> \ingroup double_lin 152* 153* ===================================================================== 154 SUBROUTINE DTBT02( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, X, 155 $ LDX, B, LDB, WORK, RESID ) 156* 157* -- LAPACK test routine (version 3.4.0) -- 158* -- LAPACK is a software package provided by Univ. of Tennessee, -- 159* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 160* November 2011 161* 162* .. Scalar Arguments .. 163 CHARACTER DIAG, TRANS, UPLO 164 INTEGER KD, LDAB, LDB, LDX, N, NRHS 165 DOUBLE PRECISION RESID 166* .. 167* .. Array Arguments .. 168 DOUBLE PRECISION AB( LDAB, * ), B( LDB, * ), WORK( * ), 169 $ X( LDX, * ) 170* .. 171* 172* ===================================================================== 173* 174* .. Parameters .. 175 DOUBLE PRECISION ZERO, ONE 176 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 177* .. 178* .. Local Scalars .. 179 INTEGER J 180 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM 181* .. 182* .. External Functions .. 183 LOGICAL LSAME 184 DOUBLE PRECISION DASUM, DLAMCH, DLANTB 185 EXTERNAL LSAME, DASUM, DLAMCH, DLANTB 186* .. 187* .. External Subroutines .. 188 EXTERNAL DAXPY, DCOPY, DTBMV 189* .. 190* .. Intrinsic Functions .. 191 INTRINSIC MAX 192* .. 193* .. Executable Statements .. 194* 195* Quick exit if N = 0 or NRHS = 0 196* 197 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN 198 RESID = ZERO 199 RETURN 200 END IF 201* 202* Compute the 1-norm of A or A'. 203* 204 IF( LSAME( TRANS, 'N' ) ) THEN 205 ANORM = DLANTB( '1', UPLO, DIAG, N, KD, AB, LDAB, WORK ) 206 ELSE 207 ANORM = DLANTB( 'I', UPLO, DIAG, N, KD, AB, LDAB, WORK ) 208 END IF 209* 210* Exit with RESID = 1/EPS if ANORM = 0. 211* 212 EPS = DLAMCH( 'Epsilon' ) 213 IF( ANORM.LE.ZERO ) THEN 214 RESID = ONE / EPS 215 RETURN 216 END IF 217* 218* Compute the maximum over the number of right hand sides of 219* norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ). 220* 221 RESID = ZERO 222 DO 10 J = 1, NRHS 223 CALL DCOPY( N, X( 1, J ), 1, WORK, 1 ) 224 CALL DTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 ) 225 CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 ) 226 BNORM = DASUM( N, WORK, 1 ) 227 XNORM = DASUM( N, X( 1, J ), 1 ) 228 IF( XNORM.LE.ZERO ) THEN 229 RESID = ONE / EPS 230 ELSE 231 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) 232 END IF 233 10 CONTINUE 234* 235 RETURN 236* 237* End of DTBT02 238* 239 END 240