1*> \brief \b SEBCHVXX 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 SEBCHVXX( THRESH, PATH ) 12* 13* .. Scalar Arguments .. 14* REAL THRESH 15* CHARACTER*3 PATH 16* .. 17* 18* 19*> \par Purpose: 20* ============= 21*> 22*> \verbatim 23*> 24*> SEBCHVXX will run S**SVXX on a series of Hilbert matrices and then 25*> compare the error bounds returned by SGESVXX 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, SGESVXX 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 SGESVXX. Let RCONDc be the RCOND returned by SGESVXX 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*> 7. 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 single_lin 93* 94* ===================================================================== 95 SUBROUTINE SEBCHVXX( 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, LDAB, 107 $ LDAFB, N_AUX_TESTS 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 119* .. Local Arrays .. 120 REAL TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS), 121 $ A(NMAX, NMAX), ACOPY(NMAX, NMAX), 122 $ INVHILB(NMAX, NMAX), R(NMAX), C(NMAX), S(NMAX), 123 $ WORK(NMAX * 5), B(NMAX, NMAX), X(NMAX, NMAX), 124 $ DIFF(NMAX, NMAX), AF(NMAX, NMAX), 125 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ), 126 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ), 127 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX ), 128 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3) 129 INTEGER IWORK(NMAX), IPIV(NMAX) 130 131* .. External Functions .. 132 REAL SLAMCH 133 134* .. External Subroutines .. 135 EXTERNAL SLAHILB, SGESVXX, SSYSVXX, SPOSVXX, SGBSVXX, 136 $ SLACPY, LSAMEN 137 LOGICAL LSAMEN 138 139* .. Intrinsic Functions .. 140 INTRINSIC SQRT, MAX, ABS 141 142* .. Parameters .. 143 INTEGER NWISE_I, CWISE_I 144 PARAMETER (NWISE_I = 1, CWISE_I = 1) 145 INTEGER BND_I, COND_I 146 PARAMETER (BND_I = 2, COND_I = 3) 147 148* Create the loop to test out the Hilbert matrices 149 150 FACT = 'E' 151 UPLO = 'U' 152 TRANS = 'N' 153 EQUED = 'N' 154 EPS = SLAMCH('Epsilon') 155 NFAIL = 0 156 N_AUX_TESTS = 0 157 LDA = NMAX 158 LDAB = (NMAX-1)+(NMAX-1)+1 159 LDAFB = 2*(NMAX-1)+(NMAX-1)+1 160 C2 = PATH( 2: 3 ) 161 162* Main loop to test the different Hilbert Matrices. 163 164 printed_guide = .false. 165 166 DO N = 1 , NMAX 167 PARAMS(1) = -1 168 PARAMS(2) = -1 169 170 KL = N-1 171 KU = N-1 172 NRHS = n 173 M = MAX(SQRT(REAL(N)), 10.0) 174 175* Generate the Hilbert matrix, its inverse, and the 176* right hand side, all scaled by the LCM(1,..,2N-1). 177 CALL SLAHILB(N, N, A, LDA, INVHILB, LDA, B, LDA, WORK, INFO) 178 179* Copy A into ACOPY. 180 CALL SLACPY('ALL', N, N, A, NMAX, ACOPY, NMAX) 181 182* Store A in band format for GB tests 183 DO J = 1, N 184 DO I = 1, KL+KU+1 185 AB( I, J ) = 0.0E+0 186 END DO 187 END DO 188 DO J = 1, N 189 DO I = MAX( 1, J-KU ), MIN( N, J+KL ) 190 AB( KU+1+I-J, J ) = A( I, J ) 191 END DO 192 END DO 193 194* Copy AB into ABCOPY. 195 DO J = 1, N 196 DO I = 1, KL+KU+1 197 ABCOPY( I, J ) = 0.0E+0 198 END DO 199 END DO 200 CALL SLACPY('ALL', KL+KU+1, N, AB, LDAB, ABCOPY, LDAB) 201 202* Call S**SVXX with default PARAMS and N_ERR_BND = 3. 203 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 204 CALL SSYSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA, 205 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND, 206 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS, 207 $ PARAMS, WORK, IWORK, INFO) 208 ELSE IF ( LSAMEN( 2, C2, 'PO' ) ) THEN 209 CALL SPOSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA, 210 $ EQUED, S, B, LDA, X, LDA, ORCOND, 211 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS, 212 $ PARAMS, WORK, IWORK, INFO) 213 ELSE IF ( LSAMEN( 2, C2, 'GB' ) ) THEN 214 CALL SGBSVXX(FACT, TRANS, N, KL, KU, NRHS, ABCOPY, 215 $ LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B, 216 $ LDA, X, LDA, ORCOND, RPVGRW, BERR, NERRBND, 217 $ ERRBND_N, ERRBND_C, NPARAMS, PARAMS, WORK, IWORK, 218 $ INFO) 219 ELSE 220 CALL SGESVXX(FACT, TRANS, N, NRHS, ACOPY, LDA, AF, LDA, 221 $ IPIV, EQUED, R, C, B, LDA, X, LDA, ORCOND, 222 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS, 223 $ PARAMS, WORK, IWORK, INFO) 224 END IF 225 226 N_AUX_TESTS = N_AUX_TESTS + 1 227 IF (ORCOND .LT. EPS) THEN 228! Either factorization failed or the matrix is flagged, and 1 <= 229! INFO <= N+1. We don't decide based on rcond anymore. 230! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN 231! NFAIL = NFAIL + 1 232! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND 233! END IF 234 ELSE 235! Either everything succeeded (INFO == 0) or some solution failed 236! to converge (INFO > N+1). 237 IF (INFO .GT. 0 .AND. INFO .LE. N+1) THEN 238 NFAIL = NFAIL + 1 239 WRITE (*, FMT=8000) C2, N, INFO, ORCOND, RCOND 240 END IF 241 END IF 242 243* Calculating the difference between S**SVXX's X and the true X. 244 DO I = 1, N 245 DO J = 1, NRHS 246 DIFF( I, J ) = X( I, J ) - INVHILB( I, J ) 247 END DO 248 END DO 249 250* Calculating the RCOND 251 RNORM = 0 252 RINORM = 0 253 IF ( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'SY' ) ) THEN 254 DO I = 1, N 255 SUMR = 0 256 SUMRI = 0 257 DO J = 1, N 258 SUMR = SUMR + ABS(S(I) * A(I,J) * S(J)) 259 SUMRI = SUMRI + ABS(INVHILB(I, J) / S(J) / S(I)) 260 END DO 261 RNORM = MAX(RNORM,SUMR) 262 RINORM = MAX(RINORM,SUMRI) 263 END DO 264 ELSE IF ( LSAMEN( 2, C2, 'GE' ) .OR. LSAMEN( 2, C2, 'GB' ) ) 265 $ THEN 266 DO I = 1, N 267 SUMR = 0 268 SUMRI = 0 269 DO J = 1, N 270 SUMR = SUMR + ABS(R(I) * A(I,J) * C(J)) 271 SUMRI = SUMRI + ABS(INVHILB(I, J) / R(J) / C(I)) 272 END DO 273 RNORM = MAX(RNORM,SUMR) 274 RINORM = MAX(RINORM,SUMRI) 275 END DO 276 END IF 277 278 RNORM = RNORM / A(1, 1) 279 RCOND = 1.0/(RNORM * RINORM) 280 281* Calculating the R for normwise rcond. 282 DO I = 1, N 283 RINV(I) = 0.0 284 END DO 285 DO J = 1, N 286 DO I = 1, N 287 RINV(I) = RINV(I) + ABS(A(I,J)) 288 END DO 289 END DO 290 291* Calculating the Normwise rcond. 292 RINORM = 0.0 293 DO I = 1, N 294 SUMRI = 0.0 295 DO J = 1, N 296 SUMRI = SUMRI + ABS(INVHILB(I,J) * RINV(J)) 297 END DO 298 RINORM = MAX(RINORM, SUMRI) 299 END DO 300 301! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm 302! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix) 303 NCOND = A(1,1) / RINORM 304 305 CONDTHRESH = M * EPS 306 ERRTHRESH = M * EPS 307 308 DO K = 1, NRHS 309 NORMT = 0.0 310 NORMDIF = 0.0 311 CWISE_ERR = 0.0 312 DO I = 1, N 313 NORMT = MAX(ABS(INVHILB(I, K)), NORMT) 314 NORMDIF = MAX(ABS(X(I,K) - INVHILB(I,K)), NORMDIF) 315 IF (INVHILB(I,K) .NE. 0.0) THEN 316 CWISE_ERR = MAX(ABS(X(I,K) - INVHILB(I,K)) 317 $ /ABS(INVHILB(I,K)), CWISE_ERR) 318 ELSE IF (X(I, K) .NE. 0.0) THEN 319 CWISE_ERR = SLAMCH('OVERFLOW') 320 END IF 321 END DO 322 IF (NORMT .NE. 0.0) THEN 323 NWISE_ERR = NORMDIF / NORMT 324 ELSE IF (NORMDIF .NE. 0.0) THEN 325 NWISE_ERR = SLAMCH('OVERFLOW') 326 ELSE 327 NWISE_ERR = 0.0 328 ENDIF 329 330 DO I = 1, N 331 RINV(I) = 0.0 332 END DO 333 DO J = 1, N 334 DO I = 1, N 335 RINV(I) = RINV(I) + ABS(A(I, J) * INVHILB(J, K)) 336 END DO 337 END DO 338 RINORM = 0.0 339 DO I = 1, N 340 SUMRI = 0.0 341 DO J = 1, N 342 SUMRI = SUMRI 343 $ + ABS(INVHILB(I, J) * RINV(J) / INVHILB(I, K)) 344 END DO 345 RINORM = MAX(RINORM, SUMRI) 346 END DO 347! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm 348! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix) 349 CCOND = A(1,1)/RINORM 350 351! Forward error bound tests 352 NWISE_BND = ERRBND_N(K + (BND_I-1)*NRHS) 353 CWISE_BND = ERRBND_C(K + (BND_I-1)*NRHS) 354 NWISE_RCOND = ERRBND_N(K + (COND_I-1)*NRHS) 355 CWISE_RCOND = ERRBND_C(K + (COND_I-1)*NRHS) 356! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond, 357! $ condthresh, ncond.ge.condthresh 358! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh 359 360 IF (NCOND .GE. CONDTHRESH) THEN 361 NGUAR = 'YES' 362 IF (NWISE_BND .GT. ERRTHRESH) THEN 363 TSTRAT(1) = 1/(2.0*EPS) 364 ELSE 365 366 IF (NWISE_BND .NE. 0.0) THEN 367 TSTRAT(1) = NWISE_ERR / NWISE_BND 368 ELSE IF (NWISE_ERR .NE. 0.0) THEN 369 TSTRAT(1) = 1/(16.0*EPS) 370 ELSE 371 TSTRAT(1) = 0.0 372 END IF 373 IF (TSTRAT(1) .GT. 1.0) THEN 374 TSTRAT(1) = 1/(4.0*EPS) 375 END IF 376 END IF 377 ELSE 378 NGUAR = 'NO' 379 IF (NWISE_BND .LT. 1.0) THEN 380 TSTRAT(1) = 1/(8.0*EPS) 381 ELSE 382 TSTRAT(1) = 1.0 383 END IF 384 END IF 385! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond, 386! $ condthresh, ccond.ge.condthresh 387! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh 388 IF (CCOND .GE. CONDTHRESH) THEN 389 CGUAR = 'YES' 390 391 IF (CWISE_BND .GT. ERRTHRESH) THEN 392 TSTRAT(2) = 1/(2.0*EPS) 393 ELSE 394 IF (CWISE_BND .NE. 0.0) THEN 395 TSTRAT(2) = CWISE_ERR / CWISE_BND 396 ELSE IF (CWISE_ERR .NE. 0.0) THEN 397 TSTRAT(2) = 1/(16.0*EPS) 398 ELSE 399 TSTRAT(2) = 0.0 400 END IF 401 IF (TSTRAT(2) .GT. 1.0) TSTRAT(2) = 1/(4.0*EPS) 402 END IF 403 ELSE 404 CGUAR = 'NO' 405 IF (CWISE_BND .LT. 1.0) THEN 406 TSTRAT(2) = 1/(8.0*EPS) 407 ELSE 408 TSTRAT(2) = 1.0 409 END IF 410 END IF 411 412! Backwards error test 413 TSTRAT(3) = BERR(K)/EPS 414 415! Condition number tests 416 TSTRAT(4) = RCOND / ORCOND 417 IF (RCOND .GE. CONDTHRESH .AND. TSTRAT(4) .LT. 1.0) 418 $ TSTRAT(4) = 1.0 / TSTRAT(4) 419 420 TSTRAT(5) = NCOND / NWISE_RCOND 421 IF (NCOND .GE. CONDTHRESH .AND. TSTRAT(5) .LT. 1.0) 422 $ TSTRAT(5) = 1.0 / TSTRAT(5) 423 424 TSTRAT(6) = CCOND / NWISE_RCOND 425 IF (CCOND .GE. CONDTHRESH .AND. TSTRAT(6) .LT. 1.0) 426 $ TSTRAT(6) = 1.0 / TSTRAT(6) 427 428 DO I = 1, NTESTS 429 IF (TSTRAT(I) .GT. THRESH) THEN 430 IF (.NOT.PRINTED_GUIDE) THEN 431 WRITE(*,*) 432 WRITE( *, 9996) 1 433 WRITE( *, 9995) 2 434 WRITE( *, 9994) 3 435 WRITE( *, 9993) 4 436 WRITE( *, 9992) 5 437 WRITE( *, 9991) 6 438 WRITE( *, 9990) 7 439 WRITE( *, 9989) 8 440 WRITE(*,*) 441 PRINTED_GUIDE = .TRUE. 442 END IF 443 WRITE( *, 9999) C2, N, K, NGUAR, CGUAR, I, TSTRAT(I) 444 NFAIL = NFAIL + 1 445 END IF 446 END DO 447 END DO 448 449c$$$ WRITE(*,*) 450c$$$ WRITE(*,*) 'Normwise Error Bounds' 451c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i) 452c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i) 453c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i) 454c$$$ WRITE(*,*) 455c$$$ WRITE(*,*) 'Componentwise Error Bounds' 456c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i) 457c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i) 458c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i) 459c$$$ print *, 'Info: ', info 460c$$$ WRITE(*,*) 461* WRITE(*,*) 'TSTRAT: ',TSTRAT 462 463 END DO 464 465 WRITE(*,*) 466 IF( NFAIL .GT. 0 ) THEN 467 WRITE(*,9998) C2, NFAIL, NTESTS*N+N_AUX_TESTS 468 ELSE 469 WRITE(*,9997) C2 470 END IF 471 9999 FORMAT( ' S', A2, 'SVXX: N =', I2, ', RHS = ', I2, 472 $ ', NWISE GUAR. = ', A, ', CWISE GUAR. = ', A, 473 $ ' test(',I1,') =', G12.5 ) 474 9998 FORMAT( ' S', A2, 'SVXX: ', I6, ' out of ', I6, 475 $ ' tests failed to pass the threshold' ) 476 9997 FORMAT( ' S', A2, 'SVXX passed the tests of error bounds' ) 477* Test ratios. 478 9996 FORMAT( 3X, I2, ': Normwise guaranteed forward error', / 5X, 479 $ 'Guaranteed case: if norm ( abs( Xc - Xt )', 480 $ ' / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then', 481 $ / 5X, 482 $ 'ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS') 483 9995 FORMAT( 3X, I2, ': Componentwise guaranteed forward error' ) 484 9994 FORMAT( 3X, I2, ': Backwards error' ) 485 9993 FORMAT( 3X, I2, ': Reciprocal condition number' ) 486 9992 FORMAT( 3X, I2, ': Reciprocal normwise condition number' ) 487 9991 FORMAT( 3X, I2, ': Raw normwise error estimate' ) 488 9990 FORMAT( 3X, I2, ': Reciprocal componentwise condition number' ) 489 9989 FORMAT( 3X, I2, ': Raw componentwise error estimate' ) 490 491 8000 FORMAT( ' S', A2, 'SVXX: N =', I2, ', INFO = ', I3, 492 $ ', ORCOND = ', G12.5, ', real RCOND = ', G12.5 ) 493* 494* End of SEBCHVXX 495* 496 END 497