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