 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SEBCHVXX( THRESH, PATH ) 00002 IMPLICIT NONE 00003 * .. Scalar Arguments .. 00004 REAL THRESH 00005 CHARACTER*3 PATH 00006 * 00007 * Purpose 00008 * ====== 00009 * 00010 * SEBCHVXX will run S**SVXX on a series of Hilbert matrices and then 00011 * compare the error bounds returned by SGESVXX to see if the returned 00012 * answer indeed falls within those bounds. 00013 * 00014 * Eight test ratios will be computed. The tests will pass if they are .LT. 00015 * THRESH. There are two cases that are determined by 1 / (SQRT( N ) * EPS). 00016 * If that value is .LE. to the component wise reciprocal condition number, 00017 * it uses the guaranteed case, other wise it uses the unguaranteed case. 00018 * 00019 * Test ratios: 00020 * Let Xc be X_computed and Xt be X_truth. 00021 * The norm used is the infinity norm. 00022 00023 * Let A be the guaranteed case and B be the unguaranteed case. 00024 * 00025 * 1. Normwise guaranteed forward error bound. 00026 * A: norm ( abs( Xc - Xt ) / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ) and 00027 * ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N),10) * EPS. 00028 * If these conditions are met, the test ratio is set to be 00029 * ERRBND( *, nwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS. 00030 * B: For this case, SGESVXX should just return 1. If it is less than 00031 * one, treat it the same as in 1A. Otherwise it fails. (Set test 00032 * ratio to ERRBND( *, nwise_i, bnd_i ) * THRESH?) 00033 * 00034 * 2. Componentwise guaranteed forward error bound. 00035 * A: norm ( abs( Xc(j) - Xt(j) ) ) / norm (Xt(j)) .LE. ERRBND( *, cwise_i, bnd_i ) 00036 * for all j .AND. ERRBND( *, cwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS. 00037 * If these conditions are met, the test ratio is set to be 00038 * ERRBND( *, cwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS. 00039 * B: Same as normwise test ratio. 00040 * 00041 * 3. Backwards error. 00042 * A: The test ratio is set to BERR/EPS. 00043 * B: Same test ratio. 00044 * 00045 * 4. Reciprocal condition number. 00046 * A: A condition number is computed with Xt and compared with the one 00047 * returned from SGESVXX. Let RCONDc be the RCOND returned by SGESVXX 00048 * and RCONDt be the RCOND from the truth value. Test ratio is set to 00049 * MAX(RCONDc/RCONDt, RCONDt/RCONDc). 00050 * B: Test ratio is set to 1 / (EPS * RCONDc). 00051 * 00052 * 5. Reciprocal normwise condition number. 00053 * A: The test ratio is set to 00054 * MAX(ERRBND( *, nwise_i, cond_i ) / NCOND, NCOND / ERRBND( *, nwise_i, cond_i )). 00055 * B: Test ratio is set to 1 / (EPS * ERRBND( *, nwise_i, cond_i )). 00056 * 00057 * 7. Reciprocal componentwise condition number. 00058 * A: Test ratio is set to 00059 * MAX(ERRBND( *, cwise_i, cond_i ) / CCOND, CCOND / ERRBND( *, cwise_i, cond_i )). 00060 * B: Test ratio is set to 1 / (EPS * ERRBND( *, cwise_i, cond_i )). 00061 * 00062 * .. Parameters .. 00063 * NMAX is determined by the largest number in the inverse of the Hilbert 00064 * matrix. Precision is exhausted when the largest entry in it is greater 00065 * than 2 to the power of the number of bits in the fraction of the data 00066 * type used plus one, which is 24 for single precision. 00067 * NMAX should be 6 for single and 11 for double. 00068 00069 INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU 00070 PARAMETER (NMAX = 6, NPARAMS = 2, NERRBND = 3, 00071 $ NTESTS = 6) 00072 00073 * .. Local Scalars .. 00074 INTEGER N, NRHS, INFO, I ,J, k, NFAIL, LDA, LDAB, 00075 $ LDAFB, N_AUX_TESTS 00076 CHARACTER FACT, TRANS, UPLO, EQUED 00077 CHARACTER*2 C2 00078 CHARACTER(3) NGUAR, CGUAR 00079 LOGICAL printed_guide 00080 REAL NCOND, CCOND, M, NORMDIF, NORMT, RCOND, 00081 $ RNORM, RINORM, SUMR, SUMRI, EPS, 00082 $ BERR(NMAX), RPVGRW, ORCOND, 00083 $ CWISE_ERR, NWISE_ERR, CWISE_BND, NWISE_BND, 00084 $ CWISE_RCOND, NWISE_RCOND, 00085 $ CONDTHRESH, ERRTHRESH 00086 00087 * .. Local Arrays .. 00088 REAL TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS), 00089 $ A(NMAX, NMAX), ACOPY(NMAX, NMAX), 00090 $ INVHILB(NMAX, NMAX), R(NMAX), C(NMAX), S(NMAX), 00091 $ WORK(NMAX * 5), B(NMAX, NMAX), X(NMAX, NMAX), 00092 $ DIFF(NMAX, NMAX), AF(NMAX, NMAX), 00093 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ), 00094 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ), 00095 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX ), 00096 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3) 00097 INTEGER IWORK(NMAX), IPIV(NMAX) 00098 00099 * .. External Functions .. 00100 REAL SLAMCH 00101 00102 * .. External Subroutines .. 00103 EXTERNAL SLAHILB, SGESVXX, SSYSVXX, SPOSVXX, SGBSVXX, 00104 $ SLACPY, LSAMEN 00105 LOGICAL LSAMEN 00106 00107 * .. Intrinsic Functions .. 00108 INTRINSIC SQRT, MAX, ABS 00109 00110 * .. Parameters .. 00111 INTEGER NWISE_I, CWISE_I 00112 PARAMETER (NWISE_I = 1, CWISE_I = 1) 00113 INTEGER BND_I, COND_I 00114 PARAMETER (BND_I = 2, COND_I = 3) 00115 00116 * Create the loop to test out the Hilbert matrices 00117 00118 FACT = 'E' 00119 UPLO = 'U' 00120 TRANS = 'N' 00121 EQUED = 'N' 00122 EPS = SLAMCH('Epsilon') 00123 NFAIL = 0 00124 N_AUX_TESTS = 0 00125 LDA = NMAX 00126 LDAB = (NMAX-1)+(NMAX-1)+1 00127 LDAFB = 2*(NMAX-1)+(NMAX-1)+1 00128 C2 = PATH( 2: 3 ) 00129 00130 * Main loop to test the different Hilbert Matrices. 00131 00132 printed_guide = .false. 00133 00134 DO N = 1 , NMAX 00135 PARAMS(1) = -1 00136 PARAMS(2) = -1 00137 00138 KL = N-1 00139 KU = N-1 00140 NRHS = n 00141 M = MAX(SQRT(REAL(N)), 10.0) 00142 00143 * Generate the Hilbert matrix, its inverse, and the 00144 * right hand side, all scaled by the LCM(1,..,2N-1). 00145 CALL SLAHILB(N, N, A, LDA, INVHILB, LDA, B, LDA, WORK, INFO) 00146 00147 * Copy A into ACOPY. 00148 CALL SLACPY('ALL', N, N, A, NMAX, ACOPY, NMAX) 00149 00150 * Store A in band format for GB tests 00151 DO J = 1, N 00152 DO I = 1, KL+KU+1 00153 AB( I, J ) = 0.0E+0 00154 END DO 00155 END DO 00156 DO J = 1, N 00157 DO I = MAX( 1, J-KU ), MIN( N, J+KL ) 00158 AB( KU+1+I-J, J ) = A( I, J ) 00159 END DO 00160 END DO 00161 00162 * Copy AB into ABCOPY. 00163 DO J = 1, N 00164 DO I = 1, KL+KU+1 00165 ABCOPY( I, J ) = 0.0E+0 00166 END DO 00167 END DO 00168 CALL SLACPY('ALL', KL+KU+1, N, AB, LDAB, ABCOPY, LDAB) 00169 00170 * Call S**SVXX with default PARAMS and N_ERR_BND = 3. 00171 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN 00172 CALL SSYSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA, 00173 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND, 00174 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS, 00175 $ PARAMS, WORK, IWORK, INFO) 00176 ELSE IF ( LSAMEN( 2, C2, 'PO' ) ) THEN 00177 CALL SPOSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA, 00178 $ EQUED, S, B, LDA, X, LDA, ORCOND, 00179 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS, 00180 $ PARAMS, WORK, IWORK, INFO) 00181 ELSE IF ( LSAMEN( 2, C2, 'GB' ) ) THEN 00182 CALL SGBSVXX(FACT, TRANS, N, KL, KU, NRHS, ABCOPY, 00183 $ LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B, 00184 $ LDA, X, LDA, ORCOND, RPVGRW, BERR, NERRBND, 00185 $ ERRBND_N, ERRBND_C, NPARAMS, PARAMS, WORK, IWORK, 00186 $ INFO) 00187 ELSE 00188 CALL SGESVXX(FACT, TRANS, N, NRHS, ACOPY, LDA, AF, LDA, 00189 $ IPIV, EQUED, R, C, B, LDA, X, LDA, ORCOND, 00190 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS, 00191 $ PARAMS, WORK, IWORK, INFO) 00192 END IF 00193 00194 N_AUX_TESTS = N_AUX_TESTS + 1 00195 IF (ORCOND .LT. EPS) THEN 00196 ! Either factorization failed or the matrix is flagged, and 1 <= 00197 ! INFO <= N+1. We don't decide based on rcond anymore. 00198 ! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN 00199 ! NFAIL = NFAIL + 1 00200 ! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND 00201 ! END IF 00202 ELSE 00203 ! Either everything succeeded (INFO == 0) or some solution failed 00204 ! to converge (INFO > N+1). 00205 IF (INFO .GT. 0 .AND. INFO .LE. N+1) THEN 00206 NFAIL = NFAIL + 1 00207 WRITE (*, FMT=8000) C2, N, INFO, ORCOND, RCOND 00208 END IF 00209 END IF 00210 00211 * Calculating the difference between S**SVXX's X and the true X. 00212 DO I = 1, N 00213 DO J = 1, NRHS 00214 DIFF( I, J ) = X( I, J ) - INVHILB( I, J ) 00215 END DO 00216 END DO 00217 00218 * Calculating the RCOND 00219 RNORM = 0 00220 RINORM = 0 00221 IF ( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'SY' ) ) THEN 00222 DO I = 1, N 00223 SUMR = 0 00224 SUMRI = 0 00225 DO J = 1, N 00226 SUMR = SUMR + ABS(S(I) * A(I,J) * S(J)) 00227 SUMRI = SUMRI + ABS(INVHILB(I, J) / S(J) / S(I)) 00228 END DO 00229 RNORM = MAX(RNORM,SUMR) 00230 RINORM = MAX(RINORM,SUMRI) 00231 END DO 00232 ELSE IF ( LSAMEN( 2, C2, 'GE' ) .OR. LSAMEN( 2, C2, 'GB' ) ) 00233 $ THEN 00234 DO I = 1, N 00235 SUMR = 0 00236 SUMRI = 0 00237 DO J = 1, N 00238 SUMR = SUMR + ABS(R(I) * A(I,J) * C(J)) 00239 SUMRI = SUMRI + ABS(INVHILB(I, J) / R(J) / C(I)) 00240 END DO 00241 RNORM = MAX(RNORM,SUMR) 00242 RINORM = MAX(RINORM,SUMRI) 00243 END DO 00244 END IF 00245 00246 RNORM = RNORM / A(1, 1) 00247 RCOND = 1.0/(RNORM * RINORM) 00248 00249 * Calculating the R for normwise rcond. 00250 DO I = 1, N 00251 RINV(I) = 0.0 00252 END DO 00253 DO J = 1, N 00254 DO I = 1, N 00255 RINV(I) = RINV(I) + ABS(A(I,J)) 00256 END DO 00257 END DO 00258 00259 * Calculating the Normwise rcond. 00260 RINORM = 0.0 00261 DO I = 1, N 00262 SUMRI = 0.0 00263 DO J = 1, N 00264 SUMRI = SUMRI + ABS(INVHILB(I,J) * RINV(J)) 00265 END DO 00266 RINORM = MAX(RINORM, SUMRI) 00267 END DO 00268 00269 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm 00270 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix) 00271 NCOND = A(1,1) / RINORM 00272 00273 CONDTHRESH = M * EPS 00274 ERRTHRESH = M * EPS 00275 00276 DO K = 1, NRHS 00277 NORMT = 0.0 00278 NORMDIF = 0.0 00279 CWISE_ERR = 0.0 00280 DO I = 1, N 00281 NORMT = MAX(ABS(INVHILB(I, K)), NORMT) 00282 NORMDIF = MAX(ABS(X(I,K) - INVHILB(I,K)), NORMDIF) 00283 IF (INVHILB(I,K) .NE. 0.0) THEN 00284 CWISE_ERR = MAX(ABS(X(I,K) - INVHILB(I,K)) 00285 $ /ABS(INVHILB(I,K)), CWISE_ERR) 00286 ELSE IF (X(I, K) .NE. 0.0) THEN 00287 CWISE_ERR = SLAMCH('OVERFLOW') 00288 END IF 00289 END DO 00290 IF (NORMT .NE. 0.0) THEN 00291 NWISE_ERR = NORMDIF / NORMT 00292 ELSE IF (NORMDIF .NE. 0.0) THEN 00293 NWISE_ERR = SLAMCH('OVERFLOW') 00294 ELSE 00295 NWISE_ERR = 0.0 00296 ENDIF 00297 00298 DO I = 1, N 00299 RINV(I) = 0.0 00300 END DO 00301 DO J = 1, N 00302 DO I = 1, N 00303 RINV(I) = RINV(I) + ABS(A(I, J) * INVHILB(J, K)) 00304 END DO 00305 END DO 00306 RINORM = 0.0 00307 DO I = 1, N 00308 SUMRI = 0.0 00309 DO J = 1, N 00310 SUMRI = SUMRI 00311 $ + ABS(INVHILB(I, J) * RINV(J) / INVHILB(I, K)) 00312 END DO 00313 RINORM = MAX(RINORM, SUMRI) 00314 END DO 00315 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm 00316 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix) 00317 CCOND = A(1,1)/RINORM 00318 00319 ! Forward error bound tests 00320 NWISE_BND = ERRBND_N(K + (BND_I-1)*NRHS) 00321 CWISE_BND = ERRBND_C(K + (BND_I-1)*NRHS) 00322 NWISE_RCOND = ERRBND_N(K + (COND_I-1)*NRHS) 00323 CWISE_RCOND = ERRBND_C(K + (COND_I-1)*NRHS) 00324 ! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond, 00325 ! $ condthresh, ncond.ge.condthresh 00326 ! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh 00327 00328 IF (NCOND .GE. CONDTHRESH) THEN 00329 NGUAR = 'YES' 00330 IF (NWISE_BND .GT. ERRTHRESH) THEN 00331 TSTRAT(1) = 1/(2.0*EPS) 00332 ELSE 00333 00334 IF (NWISE_BND .NE. 0.0) THEN 00335 TSTRAT(1) = NWISE_ERR / NWISE_BND 00336 ELSE IF (NWISE_ERR .NE. 0.0) THEN 00337 TSTRAT(1) = 1/(16.0*EPS) 00338 ELSE 00339 TSTRAT(1) = 0.0 00340 END IF 00341 IF (TSTRAT(1) .GT. 1.0) THEN 00342 TSTRAT(1) = 1/(4.0*EPS) 00343 END IF 00344 END IF 00345 ELSE 00346 NGUAR = 'NO' 00347 IF (NWISE_BND .LT. 1.0) THEN 00348 TSTRAT(1) = 1/(8.0*EPS) 00349 ELSE 00350 TSTRAT(1) = 1.0 00351 END IF 00352 END IF 00353 ! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond, 00354 ! $ condthresh, ccond.ge.condthresh 00355 ! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh 00356 IF (CCOND .GE. CONDTHRESH) THEN 00357 CGUAR = 'YES' 00358 00359 IF (CWISE_BND .GT. ERRTHRESH) THEN 00360 TSTRAT(2) = 1/(2.0*EPS) 00361 ELSE 00362 IF (CWISE_BND .NE. 0.0) THEN 00363 TSTRAT(2) = CWISE_ERR / CWISE_BND 00364 ELSE IF (CWISE_ERR .NE. 0.0) THEN 00365 TSTRAT(2) = 1/(16.0*EPS) 00366 ELSE 00367 TSTRAT(2) = 0.0 00368 END IF 00369 IF (TSTRAT(2) .GT. 1.0) TSTRAT(2) = 1/(4.0*EPS) 00370 END IF 00371 ELSE 00372 CGUAR = 'NO' 00373 IF (CWISE_BND .LT. 1.0) THEN 00374 TSTRAT(2) = 1/(8.0*EPS) 00375 ELSE 00376 TSTRAT(2) = 1.0 00377 END IF 00378 END IF 00379 00380 ! Backwards error test 00381 TSTRAT(3) = BERR(K)/EPS 00382 00383 ! Condition number tests 00384 TSTRAT(4) = RCOND / ORCOND 00385 IF (RCOND .GE. CONDTHRESH .AND. TSTRAT(4) .LT. 1.0) 00386 $ TSTRAT(4) = 1.0 / TSTRAT(4) 00387 00388 TSTRAT(5) = NCOND / NWISE_RCOND 00389 IF (NCOND .GE. CONDTHRESH .AND. TSTRAT(5) .LT. 1.0) 00390 $ TSTRAT(5) = 1.0 / TSTRAT(5) 00391 00392 TSTRAT(6) = CCOND / NWISE_RCOND 00393 IF (CCOND .GE. CONDTHRESH .AND. TSTRAT(6) .LT. 1.0) 00394 $ TSTRAT(6) = 1.0 / TSTRAT(6) 00395 00396 DO I = 1, NTESTS 00397 IF (TSTRAT(I) .GT. THRESH) THEN 00398 IF (.NOT.PRINTED_GUIDE) THEN 00399 WRITE(*,*) 00400 WRITE( *, 9996) 1 00401 WRITE( *, 9995) 2 00402 WRITE( *, 9994) 3 00403 WRITE( *, 9993) 4 00404 WRITE( *, 9992) 5 00405 WRITE( *, 9991) 6 00406 WRITE( *, 9990) 7 00407 WRITE( *, 9989) 8 00408 WRITE(*,*) 00409 PRINTED_GUIDE = .TRUE. 00410 END IF 00411 WRITE( *, 9999) C2, N, K, NGUAR, CGUAR, I, TSTRAT(I) 00412 NFAIL = NFAIL + 1 00413 END IF 00414 END DO 00415 END DO 00416 00417 c$$$ WRITE(*,*) 00418 c$$$ WRITE(*,*) 'Normwise Error Bounds' 00419 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i) 00420 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i) 00421 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i) 00422 c$$$ WRITE(*,*) 00423 c$$$ WRITE(*,*) 'Componentwise Error Bounds' 00424 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i) 00425 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i) 00426 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i) 00427 c$$$ print *, 'Info: ', info 00428 c$$$ WRITE(*,*) 00429 * WRITE(*,*) 'TSTRAT: ',TSTRAT 00430 00431 END DO 00432 00433 WRITE(*,*) 00434 IF( NFAIL .GT. 0 ) THEN 00435 WRITE(*,9998) C2, NFAIL, NTESTS*N+N_AUX_TESTS 00436 ELSE 00437 WRITE(*,9997) C2 00438 END IF 00439 9999 FORMAT( ' S', A2, 'SVXX: N =', I2, ', RHS = ', I2, 00440 $ ', NWISE GUAR. = ', A, ', CWISE GUAR. = ', A, 00441 $ ' test(',I1,') =', G12.5 ) 00442 9998 FORMAT( ' S', A2, 'SVXX: ', I6, ' out of ', I6, 00443 $ ' tests failed to pass the threshold' ) 00444 9997 FORMAT( ' S', A2, 'SVXX passed the tests of error bounds' ) 00445 * Test ratios. 00446 9996 FORMAT( 3X, I2, ': Normwise guaranteed forward error', / 5X, 00447 $ 'Guaranteed case: if norm ( abs( Xc - Xt )', 00448 $ ' / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then', 00449 $ / 5X, 00450 $ 'ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS') 00451 9995 FORMAT( 3X, I2, ': Componentwise guaranteed forward error' ) 00452 9994 FORMAT( 3X, I2, ': Backwards error' ) 00453 9993 FORMAT( 3X, I2, ': Reciprocal condition number' ) 00454 9992 FORMAT( 3X, I2, ': Reciprocal normwise condition number' ) 00455 9991 FORMAT( 3X, I2, ': Raw normwise error estimate' ) 00456 9990 FORMAT( 3X, I2, ': Reciprocal componentwise condition number' ) 00457 9989 FORMAT( 3X, I2, ': Raw componentwise error estimate' ) 00458 00459 8000 FORMAT( ' S', A2, 'SVXX: N =', I2, ', INFO = ', I3, 00460 $ ', ORCOND = ', G12.5, ', real RCOND = ', G12.5 ) 00461 00462 END
1.7.3 
