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