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