LAPACK: TESTING/LIN/sebchvxx.f Source File

Go to the documentation of this file.
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