1*> \brief \b CEBCHVXX 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 CEBCHVXX( THRESH, PATH ) 12* 13* .. Scalar Arguments .. 14* REAL THRESH 15* CHARACTER*3 PATH 16* .. 17* 18* Purpose 19* ====== 20* 21*> \details \b Purpose: 22*> \verbatim 23*> 24*> CEBCHVXX will run CGESVXX on a series of Hilbert matrices and then 25*> compare the error bounds returned by CGESVXX to see if the returned 26*> answer indeed falls within those bounds. 27*> 28*> Eight test ratios will be computed. The tests will pass if they are .LT. 29*> THRESH. There are two cases that are determined by 1 / (SQRT( N ) * EPS). 30*> If that value is .LE. to the component wise reciprocal condition number, 31*> it uses the guaranteed case, other wise it uses the unguaranteed case. 32*> 33*> Test ratios: 34*> Let Xc be X_computed and Xt be X_truth. 35*> The norm used is the infinity norm. 36*> 37*> Let A be the guaranteed case and B be the unguaranteed case. 38*> 39*> 1. Normwise guaranteed forward error bound. 40*> A: norm ( abs( Xc - Xt ) / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ) and 41*> ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N),10) * EPS. 42*> If these conditions are met, the test ratio is set to be 43*> ERRBND( *, nwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS. 44*> B: For this case, CGESVXX should just return 1. If it is less than 45*> one, treat it the same as in 1A. Otherwise it fails. (Set test 46*> ratio to ERRBND( *, nwise_i, bnd_i ) * THRESH?) 47*> 48*> 2. Componentwise guaranteed forward error bound. 49*> A: norm ( abs( Xc(j) - Xt(j) ) ) / norm (Xt(j)) .LE. ERRBND( *, cwise_i, bnd_i ) 50*> for all j .AND. ERRBND( *, cwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS. 51*> If these conditions are met, the test ratio is set to be 52*> ERRBND( *, cwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS. 53*> B: Same as normwise test ratio. 54*> 55*> 3. Backwards error. 56*> A: The test ratio is set to BERR/EPS. 57*> B: Same test ratio. 58*> 59*> 4. Reciprocal condition number. 60*> A: A condition number is computed with Xt and compared with the one 61*> returned from CGESVXX. Let RCONDc be the RCOND returned by CGESVXX 62*> and RCONDt be the RCOND from the truth value. Test ratio is set to 63*> MAX(RCONDc/RCONDt, RCONDt/RCONDc). 64*> B: Test ratio is set to 1 / (EPS * RCONDc). 65*> 66*> 5. Reciprocal normwise condition number. 67*> A: The test ratio is set to 68*> MAX(ERRBND( *, nwise_i, cond_i ) / NCOND, NCOND / ERRBND( *, nwise_i, cond_i )). 69*> B: Test ratio is set to 1 / (EPS * ERRBND( *, nwise_i, cond_i )). 70*> 71*> 6. Reciprocal componentwise condition number. 72*> A: Test ratio is set to 73*> MAX(ERRBND( *, cwise_i, cond_i ) / CCOND, CCOND / ERRBND( *, cwise_i, cond_i )). 74*> B: Test ratio is set to 1 / (EPS * ERRBND( *, cwise_i, cond_i )). 75*> 76*> .. Parameters .. 77*> NMAX is determined by the largest number in the inverse of the hilbert 78*> matrix. Precision is exhausted when the largest entry in it is greater 79*> than 2 to the power of the number of bits in the fraction of the data 80*> type used plus one, which is 24 for single precision. 81*> NMAX should be 6 for single and 11 for double. 82*> \endverbatim 83* 84* Authors: 85* ======== 86* 87*> \author Univ. of Tennessee 88*> \author Univ. of California Berkeley 89*> \author Univ. of Colorado Denver 90*> \author NAG Ltd. 91* 92*> \ingroup complex_lin 93* 94* ===================================================================== 95 SUBROUTINE CEBCHVXX( THRESH, PATH ) 96 IMPLICIT NONE 97* .. Scalar Arguments .. 98 REAL THRESH 99 CHARACTER*3 PATH 100 101 INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU 102 PARAMETER (NMAX = 6, NPARAMS = 2, NERRBND = 3, 103 $ NTESTS = 6) 104 105* .. Local Scalars .. 106 INTEGER N, NRHS, INFO, I ,J, k, NFAIL, LDA, 107 $ N_AUX_TESTS, LDAB, LDAFB 108 CHARACTER FACT, TRANS, UPLO, EQUED 109 CHARACTER*2 C2 110 CHARACTER(3) NGUAR, CGUAR 111 LOGICAL printed_guide 112 REAL NCOND, CCOND, M, NORMDIF, NORMT, RCOND, 113 $ RNORM, RINORM, SUMR, SUMRI, EPS, 114 $ BERR(NMAX), RPVGRW, ORCOND, 115 $ CWISE_ERR, NWISE_ERR, CWISE_BND, NWISE_BND, 116 $ CWISE_RCOND, NWISE_RCOND, 117 $ CONDTHRESH, ERRTHRESH 118 COMPLEX ZDUM 119 120* .. Local Arrays .. 121 REAL TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS), 122 $ S(NMAX), R(NMAX),C(NMAX),RWORK(3*NMAX), 123 $ DIFF(NMAX, NMAX), 124 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3) 125 INTEGER IPIV(NMAX) 126 COMPLEX A(NMAX,NMAX),INVHILB(NMAX,NMAX),X(NMAX,NMAX), 127 $ WORK(NMAX*3*5), AF(NMAX, NMAX),B(NMAX, NMAX), 128 $ ACOPY(NMAX, NMAX), 129 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ), 130 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ), 131 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX ) 132 133* .. External Functions .. 134 REAL SLAMCH 135 136* .. External Subroutines .. 137 EXTERNAL CLAHILB, CGESVXX, CSYSVXX, CPOSVXX, 138 $ CGBSVXX, CLACPY, LSAMEN 139 LOGICAL LSAMEN 140 141* .. Intrinsic Functions .. 142 INTRINSIC SQRT, MAX, ABS, REAL, AIMAG 143 144* .. Statement Functions .. 145 REAL CABS1 146* .. 147* .. Statement Function Definitions .. 148 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 149 150* .. Parameters .. 151 INTEGER NWISE_I, CWISE_I 152 PARAMETER (NWISE_I = 1, CWISE_I = 1) 153 INTEGER BND_I, COND_I 154 PARAMETER (BND_I = 2, COND_I = 3) 155 156* Create the loop to test out the Hilbert matrices 157 158 FACT = 'E' 159 UPLO = 'U' 160 TRANS = 'N' 161 EQUED = 'N' 162 EPS = SLAMCH('Epsilon') 163 NFAIL = 0 164 N_AUX_TESTS = 0 165 LDA = NMAX 166 LDAB = (NMAX-1)+(NMAX-1)+1 167 LDAFB = 2*(NMAX-1)+(NMAX-1)+1 168 C2 = PATH( 2: 3 ) 169 170* Main loop to test the different Hilbert Matrices. 171 172 printed_guide = .false. 173 174 DO N = 1 , NMAX 175 PARAMS(1) = -1 176 PARAMS(2) = -1 177 178 KL = N-1 179 KU = N-1 180 NRHS = n 181 M = MAX(SQRT(REAL(N)), 10.0) 182 183* Generate the Hilbert matrix, its inverse, and the 184* right hand side, all scaled by the LCM(1,..,2N-1). 185 CALL CLAHILB(N, N, A, LDA, INVHILB, LDA, B, 186 $ LDA, WORK, INFO, PATH) 187 188* Copy A into ACOPY. 189 CALL CLACPY('ALL', N, N, A, NMAX, ACOPY, NMAX) 190 191* Store A in band format for GB tests 192 DO J = 1, N 193 DO I = 1, KL+KU+1 194 AB( I, J ) = (0.0E+0,0.0E+0) 195 END DO 196 END DO 197 DO J = 1, N 198 DO I = MAX( 1, J-KU ), MIN( N, J+KL ) 199 AB( KU+1+I-J, J ) = A( I, J ) 200 END DO 201 END DO 202 203* Copy AB into ABCOPY. 204 DO J = 1, N 205 DO I = 1, KL+KU+1 206 ABCOPY( I, J ) = (0.0E+0,0.0E+0) 207 END DO 208 END DO 209 CALL CLACPY('ALL', KL+KU+1, N, AB, LDAB, ABCOPY, LDAB) 210 211* Call C**SVXX with default PARAMS and N_ERR_BND = 3. 212 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 213 CALL CSYSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA, 214 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND, 215 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS, 216 $ PARAMS, WORK, RWORK, INFO) 217 ELSE IF ( LSAMEN( 2, C2, 'PO' ) ) THEN 218 CALL CPOSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA, 219 $ EQUED, S, B, LDA, X, LDA, ORCOND, 220 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS, 221 $ PARAMS, WORK, RWORK, INFO) 222 ELSE IF ( LSAMEN( 2, C2, 'HE' ) ) THEN 223 CALL CHESVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA, 224 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND, 225 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS, 226 $ PARAMS, WORK, RWORK, INFO) 227 ELSE IF ( LSAMEN( 2, C2, 'GB' ) ) THEN 228 CALL CGBSVXX(FACT, TRANS, N, KL, KU, NRHS, ABCOPY, 229 $ LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B, 230 $ LDA, X, LDA, ORCOND, RPVGRW, BERR, NERRBND, 231 $ ERRBND_N, ERRBND_C, NPARAMS, PARAMS, WORK, RWORK, 232 $ INFO) 233 ELSE 234 CALL CGESVXX(FACT, TRANS, N, NRHS, ACOPY, LDA, AF, LDA, 235 $ IPIV, EQUED, R, C, B, LDA, X, LDA, ORCOND, 236 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS, 237 $ PARAMS, WORK, RWORK, INFO) 238 END IF 239 240 N_AUX_TESTS = N_AUX_TESTS + 1 241 IF (ORCOND .LT. EPS) THEN 242! Either factorization failed or the matrix is flagged, and 1 <= 243! INFO <= N+1. We don't decide based on rcond anymore. 244! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN 245! NFAIL = NFAIL + 1 246! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND 247! END IF 248 ELSE 249! Either everything succeeded (INFO == 0) or some solution failed 250! to converge (INFO > N+1). 251 IF (INFO .GT. 0 .AND. INFO .LE. N+1) THEN 252 NFAIL = NFAIL + 1 253 WRITE (*, FMT=8000) C2, N, INFO, ORCOND, RCOND 254 END IF 255 END IF 256 257* Calculating the difference between C**SVXX's X and the true X. 258 DO I = 1,N 259 DO J =1,NRHS 260 DIFF(I,J) = X(I,J) - INVHILB(I,J) 261 END DO 262 END DO 263 264* Calculating the RCOND 265 RNORM = 0 266 RINORM = 0 267 IF ( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'SY' ) .OR. 268 $ LSAMEN( 2, C2, 'HE' ) ) THEN 269 DO I = 1, N 270 SUMR = 0 271 SUMRI = 0 272 DO J = 1, N 273 SUMR = SUMR + S(I) * CABS1(A(I,J)) * S(J) 274 SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (S(J) * S(I)) 275 END DO 276 RNORM = MAX(RNORM,SUMR) 277 RINORM = MAX(RINORM,SUMRI) 278 END DO 279 ELSE IF ( LSAMEN( 2, C2, 'GE' ) .OR. LSAMEN( 2, C2, 'GB' ) ) 280 $ THEN 281 DO I = 1, N 282 SUMR = 0 283 SUMRI = 0 284 DO J = 1, N 285 SUMR = SUMR + R(I) * CABS1(A(I,J)) * C(J) 286 SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (R(J) * C(I)) 287 END DO 288 RNORM = MAX(RNORM,SUMR) 289 RINORM = MAX(RINORM,SUMRI) 290 END DO 291 END IF 292 293 RNORM = RNORM / CABS1(A(1, 1)) 294 RCOND = 1.0/(RNORM * RINORM) 295 296* Calculating the R for normwise rcond. 297 DO I = 1, N 298 RINV(I) = 0.0 299 END DO 300 DO J = 1, N 301 DO I = 1, N 302 RINV(I) = RINV(I) + CABS1(A(I,J)) 303 END DO 304 END DO 305 306* Calculating the Normwise rcond. 307 RINORM = 0.0 308 DO I = 1, N 309 SUMRI = 0.0 310 DO J = 1, N 311 SUMRI = SUMRI + CABS1(INVHILB(I,J) * RINV(J)) 312 END DO 313 RINORM = MAX(RINORM, SUMRI) 314 END DO 315 316! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm 317! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix) 318 NCOND = CABS1(A(1,1)) / RINORM 319 320 CONDTHRESH = M * EPS 321 ERRTHRESH = M * EPS 322 323 DO K = 1, NRHS 324 NORMT = 0.0 325 NORMDIF = 0.0 326 CWISE_ERR = 0.0 327 DO I = 1, N 328 NORMT = MAX(CABS1(INVHILB(I, K)), NORMT) 329 NORMDIF = MAX(CABS1(X(I,K) - INVHILB(I,K)), NORMDIF) 330 IF (INVHILB(I,K) .NE. 0.0) THEN 331 CWISE_ERR = MAX(CABS1(X(I,K) - INVHILB(I,K)) 332 $ /CABS1(INVHILB(I,K)), CWISE_ERR) 333 ELSE IF (X(I, K) .NE. 0.0) THEN 334 CWISE_ERR = SLAMCH('OVERFLOW') 335 END IF 336 END DO 337 IF (NORMT .NE. 0.0) THEN 338 NWISE_ERR = NORMDIF / NORMT 339 ELSE IF (NORMDIF .NE. 0.0) THEN 340 NWISE_ERR = SLAMCH('OVERFLOW') 341 ELSE 342 NWISE_ERR = 0.0 343 ENDIF 344 345 DO I = 1, N 346 RINV(I) = 0.0 347 END DO 348 DO J = 1, N 349 DO I = 1, N 350 RINV(I) = RINV(I) + CABS1(A(I, J) * INVHILB(J, K)) 351 END DO 352 END DO 353 RINORM = 0.0 354 DO I = 1, N 355 SUMRI = 0.0 356 DO J = 1, N 357 SUMRI = SUMRI 358 $ + CABS1(INVHILB(I, J) * RINV(J) / INVHILB(I, K)) 359 END DO 360 RINORM = MAX(RINORM, SUMRI) 361 END DO 362! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm 363! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix) 364 CCOND = CABS1(A(1,1))/RINORM 365 366! Forward error bound tests 367 NWISE_BND = ERRBND_N(K + (BND_I-1)*NRHS) 368 CWISE_BND = ERRBND_C(K + (BND_I-1)*NRHS) 369 NWISE_RCOND = ERRBND_N(K + (COND_I-1)*NRHS) 370 CWISE_RCOND = ERRBND_C(K + (COND_I-1)*NRHS) 371! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond, 372! $ condthresh, ncond.ge.condthresh 373! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh 374 IF (NCOND .GE. CONDTHRESH) THEN 375 NGUAR = 'YES' 376 IF (NWISE_BND .GT. ERRTHRESH) THEN 377 TSTRAT(1) = 1/(2.0*EPS) 378 ELSE 379 IF (NWISE_BND .NE. 0.0) THEN 380 TSTRAT(1) = NWISE_ERR / NWISE_BND 381 ELSE IF (NWISE_ERR .NE. 0.0) THEN 382 TSTRAT(1) = 1/(16.0*EPS) 383 ELSE 384 TSTRAT(1) = 0.0 385 END IF 386 IF (TSTRAT(1) .GT. 1.0) THEN 387 TSTRAT(1) = 1/(4.0*EPS) 388 END IF 389 END IF 390 ELSE 391 NGUAR = 'NO' 392 IF (NWISE_BND .LT. 1.0) THEN 393 TSTRAT(1) = 1/(8.0*EPS) 394 ELSE 395 TSTRAT(1) = 1.0 396 END IF 397 END IF 398! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond, 399! $ condthresh, ccond.ge.condthresh 400! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh 401 IF (CCOND .GE. CONDTHRESH) THEN 402 CGUAR = 'YES' 403 IF (CWISE_BND .GT. ERRTHRESH) THEN 404 TSTRAT(2) = 1/(2.0*EPS) 405 ELSE 406 IF (CWISE_BND .NE. 0.0) THEN 407 TSTRAT(2) = CWISE_ERR / CWISE_BND 408 ELSE IF (CWISE_ERR .NE. 0.0) THEN 409 TSTRAT(2) = 1/(16.0*EPS) 410 ELSE 411 TSTRAT(2) = 0.0 412 END IF 413 IF (TSTRAT(2) .GT. 1.0) TSTRAT(2) = 1/(4.0*EPS) 414 END IF 415 ELSE 416 CGUAR = 'NO' 417 IF (CWISE_BND .LT. 1.0) THEN 418 TSTRAT(2) = 1/(8.0*EPS) 419 ELSE 420 TSTRAT(2) = 1.0 421 END IF 422 END IF 423 424! Backwards error test 425 TSTRAT(3) = BERR(K)/EPS 426 427! Condition number tests 428 TSTRAT(4) = RCOND / ORCOND 429 IF (RCOND .GE. CONDTHRESH .AND. TSTRAT(4) .LT. 1.0) 430 $ TSTRAT(4) = 1.0 / TSTRAT(4) 431 432 TSTRAT(5) = NCOND / NWISE_RCOND 433 IF (NCOND .GE. CONDTHRESH .AND. TSTRAT(5) .LT. 1.0) 434 $ TSTRAT(5) = 1.0 / TSTRAT(5) 435 436 TSTRAT(6) = CCOND / NWISE_RCOND 437 IF (CCOND .GE. CONDTHRESH .AND. TSTRAT(6) .LT. 1.0) 438 $ TSTRAT(6) = 1.0 / TSTRAT(6) 439 440 DO I = 1, NTESTS 441 IF (TSTRAT(I) .GT. THRESH) THEN 442 IF (.NOT.PRINTED_GUIDE) THEN 443 WRITE(*,*) 444 WRITE( *, 9996) 1 445 WRITE( *, 9995) 2 446 WRITE( *, 9994) 3 447 WRITE( *, 9993) 4 448 WRITE( *, 9992) 5 449 WRITE( *, 9991) 6 450 WRITE( *, 9990) 7 451 WRITE( *, 9989) 8 452 WRITE(*,*) 453 PRINTED_GUIDE = .TRUE. 454 END IF 455 WRITE( *, 9999) C2, N, K, NGUAR, CGUAR, I, TSTRAT(I) 456 NFAIL = NFAIL + 1 457 END IF 458 END DO 459 END DO 460 461c$$$ WRITE(*,*) 462c$$$ WRITE(*,*) 'Normwise Error Bounds' 463c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i) 464c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i) 465c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i) 466c$$$ WRITE(*,*) 467c$$$ WRITE(*,*) 'Componentwise Error Bounds' 468c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i) 469c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i) 470c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i) 471c$$$ print *, 'Info: ', info 472c$$$ WRITE(*,*) 473* WRITE(*,*) 'TSTRAT: ',TSTRAT 474 475 END DO 476 477 WRITE(*,*) 478 IF( NFAIL .GT. 0 ) THEN 479 WRITE(*,9998) C2, NFAIL, NTESTS*N+N_AUX_TESTS 480 ELSE 481 WRITE(*,9997) C2 482 END IF 483 9999 FORMAT( ' C', A2, 'SVXX: N =', I2, ', RHS = ', I2, 484 $ ', NWISE GUAR. = ', A, ', CWISE GUAR. = ', A, 485 $ ' test(',I1,') =', G12.5 ) 486 9998 FORMAT( ' C', A2, 'SVXX: ', I6, ' out of ', I6, 487 $ ' tests failed to pass the threshold' ) 488 9997 FORMAT( ' C', A2, 'SVXX passed the tests of error bounds' ) 489* Test ratios. 490 9996 FORMAT( 3X, I2, ': Normwise guaranteed forward error', / 5X, 491 $ 'Guaranteed case: if norm ( abs( Xc - Xt )', 492 $ ' / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then', 493 $ / 5X, 494 $ 'ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS') 495 9995 FORMAT( 3X, I2, ': Componentwise guaranteed forward error' ) 496 9994 FORMAT( 3X, I2, ': Backwards error' ) 497 9993 FORMAT( 3X, I2, ': Reciprocal condition number' ) 498 9992 FORMAT( 3X, I2, ': Reciprocal normwise condition number' ) 499 9991 FORMAT( 3X, I2, ': Raw normwise error estimate' ) 500 9990 FORMAT( 3X, I2, ': Reciprocal componentwise condition number' ) 501 9989 FORMAT( 3X, I2, ': Raw componentwise error estimate' ) 502 503 8000 FORMAT( ' C', A2, 'SVXX: N =', I2, ', INFO = ', I3, 504 $ ', ORCOND = ', G12.5, ', real RCOND = ', G12.5 ) 505* 506* End of CEBCHVXX 507* 508 END 509