1*> \brief \b DTRT02 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 DTRT02( UPLO, TRANS, DIAG, N, NRHS, A, LDA, X, LDX, B, 12* LDB, WORK, RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER DIAG, TRANS, UPLO 16* INTEGER LDA, LDB, LDX, N, NRHS 17* DOUBLE PRECISION RESID 18* .. 19* .. Array Arguments .. 20* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ), 21* $ X( LDX, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> DTRT02 computes the residual for the computed solution to a 31*> triangular system of linear equations A*x = b or A'*x = b. 32*> Here A is a triangular matrix, A' is the transpose of A, and x and b 33*> are N by NRHS matrices. The test ratio is the maximum over the 34*> 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] NRHS 74*> \verbatim 75*> NRHS is INTEGER 76*> The number of right hand sides, i.e., the number of columns 77*> of the matrices X and B. NRHS >= 0. 78*> \endverbatim 79*> 80*> \param[in] A 81*> \verbatim 82*> A is DOUBLE PRECISION array, dimension (LDA,N) 83*> The triangular matrix A. If UPLO = 'U', the leading n by n 84*> upper triangular part of the array A contains the upper 85*> triangular matrix, and the strictly lower triangular part of 86*> A is not referenced. If UPLO = 'L', the leading n by n lower 87*> triangular part of the array A contains the lower triangular 88*> matrix, and the strictly upper triangular part of A is not 89*> referenced. If DIAG = 'U', the diagonal elements of A are 90*> also not referenced and are assumed to be 1. 91*> \endverbatim 92*> 93*> \param[in] LDA 94*> \verbatim 95*> LDA is INTEGER 96*> The leading dimension of the array A. LDA >= max(1,N). 97*> \endverbatim 98*> 99*> \param[in] X 100*> \verbatim 101*> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 102*> The computed solution vectors for the system of linear 103*> equations. 104*> \endverbatim 105*> 106*> \param[in] LDX 107*> \verbatim 108*> LDX is INTEGER 109*> The leading dimension of the array X. LDX >= max(1,N). 110*> \endverbatim 111*> 112*> \param[in] B 113*> \verbatim 114*> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 115*> The right hand side vectors for the system of linear 116*> equations. 117*> \endverbatim 118*> 119*> \param[in] LDB 120*> \verbatim 121*> LDB is INTEGER 122*> The leading dimension of the array B. LDB >= max(1,N). 123*> \endverbatim 124*> 125*> \param[out] WORK 126*> \verbatim 127*> WORK is DOUBLE PRECISION array, dimension (N) 128*> \endverbatim 129*> 130*> \param[out] RESID 131*> \verbatim 132*> RESID is DOUBLE PRECISION 133*> The maximum over the number of right hand sides of 134*> norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ). 135*> \endverbatim 136* 137* Authors: 138* ======== 139* 140*> \author Univ. of Tennessee 141*> \author Univ. of California Berkeley 142*> \author Univ. of Colorado Denver 143*> \author NAG Ltd. 144* 145*> \date November 2011 146* 147*> \ingroup double_lin 148* 149* ===================================================================== 150 SUBROUTINE DTRT02( UPLO, TRANS, DIAG, N, NRHS, A, LDA, X, LDX, B, 151 $ LDB, WORK, RESID ) 152* 153* -- LAPACK test routine (version 3.4.0) -- 154* -- LAPACK is a software package provided by Univ. of Tennessee, -- 155* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 156* November 2011 157* 158* .. Scalar Arguments .. 159 CHARACTER DIAG, TRANS, UPLO 160 INTEGER LDA, LDB, LDX, N, NRHS 161 DOUBLE PRECISION RESID 162* .. 163* .. Array Arguments .. 164 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ), 165 $ X( LDX, * ) 166* .. 167* 168* ===================================================================== 169* 170* .. Parameters .. 171 DOUBLE PRECISION ZERO, ONE 172 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 173* .. 174* .. Local Scalars .. 175 INTEGER J 176 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM 177* .. 178* .. External Functions .. 179 LOGICAL LSAME 180 DOUBLE PRECISION DASUM, DLAMCH, DLANTR 181 EXTERNAL LSAME, DASUM, DLAMCH, DLANTR 182* .. 183* .. External Subroutines .. 184 EXTERNAL DAXPY, DCOPY, DTRMV 185* .. 186* .. Intrinsic Functions .. 187 INTRINSIC MAX 188* .. 189* .. Executable Statements .. 190* 191* Quick exit if N = 0 or NRHS = 0 192* 193 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN 194 RESID = ZERO 195 RETURN 196 END IF 197* 198* Compute the 1-norm of A or A'. 199* 200 IF( LSAME( TRANS, 'N' ) ) THEN 201 ANORM = DLANTR( '1', UPLO, DIAG, N, N, A, LDA, WORK ) 202 ELSE 203 ANORM = DLANTR( 'I', UPLO, DIAG, N, N, A, LDA, WORK ) 204 END IF 205* 206* Exit with RESID = 1/EPS if ANORM = 0. 207* 208 EPS = DLAMCH( 'Epsilon' ) 209 IF( ANORM.LE.ZERO ) THEN 210 RESID = ONE / EPS 211 RETURN 212 END IF 213* 214* Compute the maximum over the number of right hand sides of 215* norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ) 216* 217 RESID = ZERO 218 DO 10 J = 1, NRHS 219 CALL DCOPY( N, X( 1, J ), 1, WORK, 1 ) 220 CALL DTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 ) 221 CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 ) 222 BNORM = DASUM( N, WORK, 1 ) 223 XNORM = DASUM( N, X( 1, J ), 1 ) 224 IF( XNORM.LE.ZERO ) THEN 225 RESID = ONE / EPS 226 ELSE 227 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) 228 END IF 229 10 CONTINUE 230* 231 RETURN 232* 233* End of DTRT02 234* 235 END 236