1*> \brief \b ZPPT05 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 ZPPT05( 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* DOUBLE PRECISION BERR( * ), FERR( * ), RESLTS( * ) 20* COMPLEX*16 AP( * ), B( LDB, * ), X( LDX, * ), 21* $ XACT( LDXACT, * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> ZPPT05 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*16 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*16 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*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*> \date December 2016 153* 154*> \ingroup complex16_lin 155* 156* ===================================================================== 157 SUBROUTINE ZPPT05( UPLO, N, NRHS, AP, B, LDB, X, LDX, XACT, 158 $ LDXACT, FERR, BERR, RESLTS ) 159* 160* -- LAPACK test routine (version 3.7.0) -- 161* -- LAPACK is a software package provided by Univ. of Tennessee, -- 162* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 163* December 2016 164* 165* .. Scalar Arguments .. 166 CHARACTER UPLO 167 INTEGER LDB, LDX, LDXACT, N, NRHS 168* .. 169* .. Array Arguments .. 170 DOUBLE PRECISION BERR( * ), FERR( * ), RESLTS( * ) 171 COMPLEX*16 AP( * ), B( LDB, * ), X( LDX, * ), 172 $ XACT( LDXACT, * ) 173* .. 174* 175* ===================================================================== 176* 177* .. Parameters .. 178 DOUBLE PRECISION ZERO, ONE 179 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 180* .. 181* .. Local Scalars .. 182 LOGICAL UPPER 183 INTEGER I, IMAX, J, JC, K 184 DOUBLE PRECISION AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM 185 COMPLEX*16 ZDUM 186* .. 187* .. External Functions .. 188 LOGICAL LSAME 189 INTEGER IZAMAX 190 DOUBLE PRECISION DLAMCH 191 EXTERNAL LSAME, IZAMAX, DLAMCH 192* .. 193* .. Intrinsic Functions .. 194 INTRINSIC ABS, DBLE, DIMAG, MAX, MIN 195* .. 196* .. Statement Functions .. 197 DOUBLE PRECISION CABS1 198* .. 199* .. Statement Function definitions .. 200 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 201* .. 202* .. Executable Statements .. 203* 204* Quick exit if N = 0 or NRHS = 0. 205* 206 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN 207 RESLTS( 1 ) = ZERO 208 RESLTS( 2 ) = ZERO 209 RETURN 210 END IF 211* 212 EPS = DLAMCH( 'Epsilon' ) 213 UNFL = DLAMCH( 'Safe minimum' ) 214 OVFL = ONE / UNFL 215 UPPER = LSAME( UPLO, 'U' ) 216* 217* Test 1: Compute the maximum of 218* norm(X - XACT) / ( norm(X) * FERR ) 219* over all the vectors X and XACT using the infinity-norm. 220* 221 ERRBND = ZERO 222 DO 30 J = 1, NRHS 223 IMAX = IZAMAX( N, X( 1, J ), 1 ) 224 XNORM = MAX( CABS1( X( IMAX, J ) ), UNFL ) 225 DIFF = ZERO 226 DO 10 I = 1, N 227 DIFF = MAX( DIFF, CABS1( X( I, J )-XACT( I, J ) ) ) 228 10 CONTINUE 229* 230 IF( XNORM.GT.ONE ) THEN 231 GO TO 20 232 ELSE IF( DIFF.LE.OVFL*XNORM ) THEN 233 GO TO 20 234 ELSE 235 ERRBND = ONE / EPS 236 GO TO 30 237 END IF 238* 239 20 CONTINUE 240 IF( DIFF / XNORM.LE.FERR( J ) ) THEN 241 ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) ) 242 ELSE 243 ERRBND = ONE / EPS 244 END IF 245 30 CONTINUE 246 RESLTS( 1 ) = ERRBND 247* 248* Test 2: Compute the maximum of BERR / ( (n+1)*EPS + (*) ), where 249* (*) = (n+1)*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i ) 250* 251 DO 90 K = 1, NRHS 252 DO 80 I = 1, N 253 TMP = CABS1( B( I, K ) ) 254 IF( UPPER ) THEN 255 JC = ( ( I-1 )*I ) / 2 256 DO 40 J = 1, I - 1 257 TMP = TMP + CABS1( AP( JC+J ) )*CABS1( X( J, K ) ) 258 40 CONTINUE 259 TMP = TMP + ABS( DBLE( AP( JC+I ) ) )*CABS1( X( I, K ) ) 260 JC = JC + I + I 261 DO 50 J = I + 1, N 262 TMP = TMP + CABS1( AP( JC ) )*CABS1( X( J, K ) ) 263 JC = JC + J 264 50 CONTINUE 265 ELSE 266 JC = I 267 DO 60 J = 1, I - 1 268 TMP = TMP + CABS1( AP( JC ) )*CABS1( X( J, K ) ) 269 JC = JC + N - J 270 60 CONTINUE 271 TMP = TMP + ABS( DBLE( AP( JC ) ) )*CABS1( X( I, K ) ) 272 DO 70 J = I + 1, N 273 TMP = TMP + CABS1( AP( JC+J-I ) )*CABS1( X( J, K ) ) 274 70 CONTINUE 275 END IF 276 IF( I.EQ.1 ) THEN 277 AXBI = TMP 278 ELSE 279 AXBI = MIN( AXBI, TMP ) 280 END IF 281 80 CONTINUE 282 TMP = BERR( K ) / ( ( N+1 )*EPS+( N+1 )*UNFL / 283 $ MAX( AXBI, ( N+1 )*UNFL ) ) 284 IF( K.EQ.1 ) THEN 285 RESLTS( 2 ) = TMP 286 ELSE 287 RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP ) 288 END IF 289 90 CONTINUE 290* 291 RETURN 292* 293* End of ZPPT05 294* 295 END 296