1*> \brief \b STPT03 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 STPT03( UPLO, TRANS, DIAG, N, NRHS, AP, SCALE, CNORM, 12* TSCAL, X, LDX, B, LDB, WORK, RESID ) 13* 14* .. Scalar Arguments .. 15* CHARACTER DIAG, TRANS, UPLO 16* INTEGER LDB, LDX, N, NRHS 17* REAL RESID, SCALE, TSCAL 18* .. 19* .. Array Arguments .. 20* REAL AP( * ), B( LDB, * ), CNORM( * ), WORK( * ), 21* $ X( LDX, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> STPT03 computes the residual for the solution to a scaled triangular 31*> system of equations A*x = s*b or A'*x = s*b when the triangular 32*> matrix A is stored in packed format. Here A' is the transpose of A, 33*> s is a scalar, and x and b are N by NRHS matrices. The test ratio is 34*> the maximum over the number of right hand sides of 35*> norm(s*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 = s*b (No transpose) 55*> = 'T': A'*x = s*b (Transpose) 56*> = 'C': A'*x = s*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] AP 81*> \verbatim 82*> AP is REAL array, dimension (N*(N+1)/2) 83*> The upper or lower triangular matrix A, packed columnwise in 84*> a linear array. The j-th column of A is stored in the array 85*> AP as follows: 86*> if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j; 87*> if UPLO = 'L', 88*> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n. 89*> \endverbatim 90*> 91*> \param[in] SCALE 92*> \verbatim 93*> SCALE is REAL 94*> The scaling factor s used in solving the triangular system. 95*> \endverbatim 96*> 97*> \param[in] CNORM 98*> \verbatim 99*> CNORM is REAL array, dimension (N) 100*> The 1-norms of the columns of A, not counting the diagonal. 101*> \endverbatim 102*> 103*> \param[in] TSCAL 104*> \verbatim 105*> TSCAL is REAL 106*> The scaling factor used in computing the 1-norms in CNORM. 107*> CNORM actually contains the column norms of TSCAL*A. 108*> \endverbatim 109*> 110*> \param[in] X 111*> \verbatim 112*> X is REAL array, dimension (LDX,NRHS) 113*> The computed solution vectors for the system of linear 114*> equations. 115*> \endverbatim 116*> 117*> \param[in] LDX 118*> \verbatim 119*> LDX is INTEGER 120*> The leading dimension of the array X. LDX >= max(1,N). 121*> \endverbatim 122*> 123*> \param[in] B 124*> \verbatim 125*> B is REAL array, dimension (LDB,NRHS) 126*> The right hand side vectors for the system of linear 127*> equations. 128*> \endverbatim 129*> 130*> \param[in] LDB 131*> \verbatim 132*> LDB is INTEGER 133*> The leading dimension of the array B. LDB >= max(1,N). 134*> \endverbatim 135*> 136*> \param[out] WORK 137*> \verbatim 138*> WORK is REAL array, dimension (N) 139*> \endverbatim 140*> 141*> \param[out] RESID 142*> \verbatim 143*> RESID is REAL 144*> The maximum over the number of right hand sides of 145*> norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ). 146*> \endverbatim 147* 148* Authors: 149* ======== 150* 151*> \author Univ. of Tennessee 152*> \author Univ. of California Berkeley 153*> \author Univ. of Colorado Denver 154*> \author NAG Ltd. 155* 156*> \date December 2016 157* 158*> \ingroup single_lin 159* 160* ===================================================================== 161 SUBROUTINE STPT03( UPLO, TRANS, DIAG, N, NRHS, AP, SCALE, CNORM, 162 $ TSCAL, X, LDX, B, LDB, WORK, RESID ) 163* 164* -- LAPACK test routine (version 3.7.0) -- 165* -- LAPACK is a software package provided by Univ. of Tennessee, -- 166* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 167* December 2016 168* 169* .. Scalar Arguments .. 170 CHARACTER DIAG, TRANS, UPLO 171 INTEGER LDB, LDX, N, NRHS 172 REAL RESID, SCALE, TSCAL 173* .. 174* .. Array Arguments .. 175 REAL AP( * ), B( LDB, * ), CNORM( * ), WORK( * ), 176 $ X( LDX, * ) 177* .. 178* 179* ===================================================================== 180* 181* .. Parameters .. 182 REAL ONE, ZERO 183 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 184* .. 185* .. Local Scalars .. 186 INTEGER IX, J, JJ 187 REAL BIGNUM, EPS, ERR, SMLNUM, TNORM, XNORM, XSCAL 188* .. 189* .. External Functions .. 190 LOGICAL LSAME 191 INTEGER ISAMAX 192 REAL SLAMCH 193 EXTERNAL LSAME, ISAMAX, SLAMCH 194* .. 195* .. External Subroutines .. 196 EXTERNAL SAXPY, SCOPY, SLABAD, SSCAL, STPMV 197* .. 198* .. Intrinsic Functions .. 199 INTRINSIC ABS, MAX, REAL 200* .. 201* .. Executable Statements .. 202* 203* Quick exit if N = 0. 204* 205 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN 206 RESID = ZERO 207 RETURN 208 END IF 209 EPS = SLAMCH( 'Epsilon' ) 210 SMLNUM = SLAMCH( 'Safe minimum' ) 211 BIGNUM = ONE / SMLNUM 212 CALL SLABAD( SMLNUM, BIGNUM ) 213* 214* Compute the norm of the triangular matrix A using the column 215* norms already computed by SLATPS. 216* 217 TNORM = ZERO 218 IF( LSAME( DIAG, 'N' ) ) THEN 219 IF( LSAME( UPLO, 'U' ) ) THEN 220 JJ = 1 221 DO 10 J = 1, N 222 TNORM = MAX( TNORM, TSCAL*ABS( AP( JJ ) )+CNORM( J ) ) 223 JJ = JJ + J + 1 224 10 CONTINUE 225 ELSE 226 JJ = 1 227 DO 20 J = 1, N 228 TNORM = MAX( TNORM, TSCAL*ABS( AP( JJ ) )+CNORM( J ) ) 229 JJ = JJ + N - J + 1 230 20 CONTINUE 231 END IF 232 ELSE 233 DO 30 J = 1, N 234 TNORM = MAX( TNORM, TSCAL+CNORM( J ) ) 235 30 CONTINUE 236 END IF 237* 238* Compute the maximum over the number of right hand sides of 239* norm(op(A)*x - s*b) / ( norm(op(A)) * norm(x) * EPS ). 240* 241 RESID = ZERO 242 DO 40 J = 1, NRHS 243 CALL SCOPY( N, X( 1, J ), 1, WORK, 1 ) 244 IX = ISAMAX( N, WORK, 1 ) 245 XNORM = MAX( ONE, ABS( X( IX, J ) ) ) 246 XSCAL = ( ONE / XNORM ) / REAL( N ) 247 CALL SSCAL( N, XSCAL, WORK, 1 ) 248 CALL STPMV( UPLO, TRANS, DIAG, N, AP, WORK, 1 ) 249 CALL SAXPY( N, -SCALE*XSCAL, B( 1, J ), 1, WORK, 1 ) 250 IX = ISAMAX( N, WORK, 1 ) 251 ERR = TSCAL*ABS( WORK( IX ) ) 252 IX = ISAMAX( N, X( 1, J ), 1 ) 253 XNORM = ABS( X( IX, J ) ) 254 IF( ERR*SMLNUM.LE.XNORM ) THEN 255 IF( XNORM.GT.ZERO ) 256 $ ERR = ERR / XNORM 257 ELSE 258 IF( ERR.GT.ZERO ) 259 $ ERR = ONE / EPS 260 END IF 261 IF( ERR*SMLNUM.LE.TNORM ) THEN 262 IF( TNORM.GT.ZERO ) 263 $ ERR = ERR / TNORM 264 ELSE 265 IF( ERR.GT.ZERO ) 266 $ ERR = ONE / EPS 267 END IF 268 RESID = MAX( RESID, ERR ) 269 40 CONTINUE 270* 271 RETURN 272* 273* End of STPT03 274* 275 END 276