1*> \brief \b CTRT05 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 CTRT05( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X, 12* LDX, XACT, LDXACT, FERR, BERR, RESLTS ) 13* 14* .. Scalar Arguments .. 15* CHARACTER DIAG, TRANS, UPLO 16* INTEGER LDA, LDB, LDX, LDXACT, N, NRHS 17* .. 18* .. Array Arguments .. 19* REAL BERR( * ), FERR( * ), RESLTS( * ) 20* COMPLEX A( LDA, * ), B( LDB, * ), X( LDX, * ), 21* $ XACT( LDXACT, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> CTRT05 tests the error bounds from iterative refinement for the 31*> computed solution to a system of equations A*X = B, where A is a 32*> triangular n by n matrix. 33*> 34*> RESLTS(1) = test of the error bound 35*> = norm(X - XACT) / ( norm(X) * FERR ) 36*> 37*> A large value is returned if this ratio is not less than one. 38*> 39*> RESLTS(2) = residual from the iterative refinement routine 40*> = the maximum of BERR / ( (n+1)*EPS + (*) ), where 41*> (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i ) 42*> \endverbatim 43* 44* Arguments: 45* ========== 46* 47*> \param[in] UPLO 48*> \verbatim 49*> UPLO is CHARACTER*1 50*> Specifies whether the matrix A is upper or lower triangular. 51*> = 'U': Upper triangular 52*> = 'L': Lower triangular 53*> \endverbatim 54*> 55*> \param[in] TRANS 56*> \verbatim 57*> TRANS is CHARACTER*1 58*> Specifies the form of the system of equations. 59*> = 'N': A * X = B (No transpose) 60*> = 'T': A'* X = B (Transpose) 61*> = 'C': A'* X = B (Conjugate transpose = Transpose) 62*> \endverbatim 63*> 64*> \param[in] DIAG 65*> \verbatim 66*> DIAG is CHARACTER*1 67*> Specifies whether or not the matrix A is unit triangular. 68*> = 'N': Non-unit triangular 69*> = 'U': Unit triangular 70*> \endverbatim 71*> 72*> \param[in] N 73*> \verbatim 74*> N is INTEGER 75*> The number of rows of the matrices X, B, and XACT, and the 76*> order of the matrix A. N >= 0. 77*> \endverbatim 78*> 79*> \param[in] NRHS 80*> \verbatim 81*> NRHS is INTEGER 82*> The number of columns of the matrices X, B, and XACT. 83*> NRHS >= 0. 84*> \endverbatim 85*> 86*> \param[in] A 87*> \verbatim 88*> A is COMPLEX array, dimension (LDA,N) 89*> The triangular matrix A. If UPLO = 'U', the leading n by n 90*> upper triangular part of the array A contains the upper 91*> triangular matrix, and the strictly lower triangular part of 92*> A is not referenced. If UPLO = 'L', the leading n by n lower 93*> triangular part of the array A contains the lower triangular 94*> matrix, and the strictly upper triangular part of A is not 95*> referenced. If DIAG = 'U', the diagonal elements of A are 96*> also not referenced and are assumed to be 1. 97*> \endverbatim 98*> 99*> \param[in] LDA 100*> \verbatim 101*> LDA is INTEGER 102*> The leading dimension of the array A. LDA >= max(1,N). 103*> \endverbatim 104*> 105*> \param[in] B 106*> \verbatim 107*> B is COMPLEX array, dimension (LDB,NRHS) 108*> The right hand side vectors for the system of linear 109*> equations. 110*> \endverbatim 111*> 112*> \param[in] LDB 113*> \verbatim 114*> LDB is INTEGER 115*> The leading dimension of the array B. LDB >= max(1,N). 116*> \endverbatim 117*> 118*> \param[in] X 119*> \verbatim 120*> X is COMPLEX array, dimension (LDX,NRHS) 121*> The computed solution vectors. Each vector is stored as a 122*> column of the matrix X. 123*> \endverbatim 124*> 125*> \param[in] LDX 126*> \verbatim 127*> LDX is INTEGER 128*> The leading dimension of the array X. LDX >= max(1,N). 129*> \endverbatim 130*> 131*> \param[in] XACT 132*> \verbatim 133*> XACT is COMPLEX array, dimension (LDX,NRHS) 134*> The exact solution vectors. Each vector is stored as a 135*> column of the matrix XACT. 136*> \endverbatim 137*> 138*> \param[in] LDXACT 139*> \verbatim 140*> LDXACT is INTEGER 141*> The leading dimension of the array XACT. LDXACT >= max(1,N). 142*> \endverbatim 143*> 144*> \param[in] FERR 145*> \verbatim 146*> FERR is REAL array, dimension (NRHS) 147*> The estimated forward error bounds for each solution vector 148*> X. If XTRUE is the true solution, FERR bounds the magnitude 149*> of the largest entry in (X - XTRUE) divided by the magnitude 150*> of the largest entry in X. 151*> \endverbatim 152*> 153*> \param[in] BERR 154*> \verbatim 155*> BERR is REAL array, dimension (NRHS) 156*> The componentwise relative backward error of each solution 157*> vector (i.e., the smallest relative change in any entry of A 158*> or B that makes X an exact solution). 159*> \endverbatim 160*> 161*> \param[out] RESLTS 162*> \verbatim 163*> RESLTS is REAL array, dimension (2) 164*> The maximum over the NRHS solution vectors of the ratios: 165*> RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR ) 166*> RESLTS(2) = BERR / ( (n+1)*EPS + (*) ) 167*> \endverbatim 168* 169* Authors: 170* ======== 171* 172*> \author Univ. of Tennessee 173*> \author Univ. of California Berkeley 174*> \author Univ. of Colorado Denver 175*> \author NAG Ltd. 176* 177*> \ingroup complex_lin 178* 179* ===================================================================== 180 SUBROUTINE CTRT05( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X, 181 $ LDX, XACT, LDXACT, FERR, BERR, RESLTS ) 182* 183* -- LAPACK test routine -- 184* -- LAPACK is a software package provided by Univ. of Tennessee, -- 185* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 186* 187* .. Scalar Arguments .. 188 CHARACTER DIAG, TRANS, UPLO 189 INTEGER LDA, LDB, LDX, LDXACT, N, NRHS 190* .. 191* .. Array Arguments .. 192 REAL BERR( * ), FERR( * ), RESLTS( * ) 193 COMPLEX A( LDA, * ), B( LDB, * ), X( LDX, * ), 194 $ XACT( LDXACT, * ) 195* .. 196* 197* ===================================================================== 198* 199* .. Parameters .. 200 REAL ZERO, ONE 201 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 202* .. 203* .. Local Scalars .. 204 LOGICAL NOTRAN, UNIT, UPPER 205 INTEGER I, IFU, IMAX, J, K 206 REAL AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM 207 COMPLEX ZDUM 208* .. 209* .. External Functions .. 210 LOGICAL LSAME 211 INTEGER ICAMAX 212 REAL SLAMCH 213 EXTERNAL LSAME, ICAMAX, SLAMCH 214* .. 215* .. Intrinsic Functions .. 216 INTRINSIC ABS, AIMAG, MAX, MIN, REAL 217* .. 218* .. Statement Functions .. 219 REAL CABS1 220* .. 221* .. Statement Function definitions .. 222 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 223* .. 224* .. Executable Statements .. 225* 226* Quick exit if N = 0 or NRHS = 0. 227* 228 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN 229 RESLTS( 1 ) = ZERO 230 RESLTS( 2 ) = ZERO 231 RETURN 232 END IF 233* 234 EPS = SLAMCH( 'Epsilon' ) 235 UNFL = SLAMCH( 'Safe minimum' ) 236 OVFL = ONE / UNFL 237 UPPER = LSAME( UPLO, 'U' ) 238 NOTRAN = LSAME( TRANS, 'N' ) 239 UNIT = LSAME( DIAG, 'U' ) 240* 241* Test 1: Compute the maximum of 242* norm(X - XACT) / ( norm(X) * FERR ) 243* over all the vectors X and XACT using the infinity-norm. 244* 245 ERRBND = ZERO 246 DO 30 J = 1, NRHS 247 IMAX = ICAMAX( N, X( 1, J ), 1 ) 248 XNORM = MAX( CABS1( X( IMAX, J ) ), UNFL ) 249 DIFF = ZERO 250 DO 10 I = 1, N 251 DIFF = MAX( DIFF, CABS1( X( I, J )-XACT( I, J ) ) ) 252 10 CONTINUE 253* 254 IF( XNORM.GT.ONE ) THEN 255 GO TO 20 256 ELSE IF( DIFF.LE.OVFL*XNORM ) THEN 257 GO TO 20 258 ELSE 259 ERRBND = ONE / EPS 260 GO TO 30 261 END IF 262* 263 20 CONTINUE 264 IF( DIFF / XNORM.LE.FERR( J ) ) THEN 265 ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) ) 266 ELSE 267 ERRBND = ONE / EPS 268 END IF 269 30 CONTINUE 270 RESLTS( 1 ) = ERRBND 271* 272* Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where 273* (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i ) 274* 275 IFU = 0 276 IF( UNIT ) 277 $ IFU = 1 278 DO 90 K = 1, NRHS 279 DO 80 I = 1, N 280 TMP = CABS1( B( I, K ) ) 281 IF( UPPER ) THEN 282 IF( .NOT.NOTRAN ) THEN 283 DO 40 J = 1, I - IFU 284 TMP = TMP + CABS1( A( J, I ) )*CABS1( X( J, K ) ) 285 40 CONTINUE 286 IF( UNIT ) 287 $ TMP = TMP + CABS1( X( I, K ) ) 288 ELSE 289 IF( UNIT ) 290 $ TMP = TMP + CABS1( X( I, K ) ) 291 DO 50 J = I + IFU, N 292 TMP = TMP + CABS1( A( I, J ) )*CABS1( X( J, K ) ) 293 50 CONTINUE 294 END IF 295 ELSE 296 IF( NOTRAN ) THEN 297 DO 60 J = 1, I - IFU 298 TMP = TMP + CABS1( A( I, J ) )*CABS1( X( J, K ) ) 299 60 CONTINUE 300 IF( UNIT ) 301 $ TMP = TMP + CABS1( X( I, K ) ) 302 ELSE 303 IF( UNIT ) 304 $ TMP = TMP + CABS1( X( I, K ) ) 305 DO 70 J = I + IFU, N 306 TMP = TMP + CABS1( A( J, I ) )*CABS1( X( J, K ) ) 307 70 CONTINUE 308 END IF 309 END IF 310 IF( I.EQ.1 ) THEN 311 AXBI = TMP 312 ELSE 313 AXBI = MIN( AXBI, TMP ) 314 END IF 315 80 CONTINUE 316 TMP = BERR( K ) / ( ( N+1 )*EPS+( N+1 )*UNFL / 317 $ MAX( AXBI, ( N+1 )*UNFL ) ) 318 IF( K.EQ.1 ) THEN 319 RESLTS( 2 ) = TMP 320 ELSE 321 RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP ) 322 END IF 323 90 CONTINUE 324* 325 RETURN 326* 327* End of CTRT05 328* 329 END 330