00001 SUBROUTINE DEBCHVXX( THRESH, PATH )
00002 IMPLICIT NONE
00003
00004 DOUBLE PRECISION THRESH
00005 CHARACTER*3 PATH
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069 INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU
00070 PARAMETER (NMAX = 10, NPARAMS = 2, NERRBND = 3,
00071 $ NTESTS = 6)
00072
00073
00074 INTEGER N, NRHS, INFO, I ,J, k, NFAIL, LDA,
00075 $ N_AUX_TESTS, LDAB, LDAFB
00076 CHARACTER FACT, TRANS, UPLO, EQUED
00077 CHARACTER*2 C2
00078 CHARACTER(3) NGUAR, CGUAR
00079 LOGICAL printed_guide
00080 DOUBLE PRECISION 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
00088 DOUBLE PRECISION TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
00089 $ S(NMAX),R(NMAX),C(NMAX), DIFF(NMAX, NMAX),
00090 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3),
00091 $ A(NMAX,NMAX),INVHILB(NMAX,NMAX),X(NMAX,NMAX),
00092 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
00093 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
00094 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX ),
00095 $ WORK(NMAX*3*5), AF(NMAX, NMAX),B(NMAX, NMAX),
00096 $ ACOPY(NMAX, NMAX)
00097 INTEGER IPIV(NMAX), IWORK(3*NMAX)
00098
00099
00100 DOUBLE PRECISION DLAMCH
00101
00102
00103 EXTERNAL DLAHILB, DGESVXX, DPOSVXX, DSYSVXX,
00104 $ DGBSVXX, DLACPY, LSAMEN
00105 LOGICAL LSAMEN
00106
00107
00108 INTRINSIC SQRT, MAX, ABS, DBLE
00109
00110
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
00117
00118 FACT = 'E'
00119 UPLO = 'U'
00120 TRANS = 'N'
00121 EQUED = 'N'
00122 EPS = DLAMCH('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
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(DBLE(N)), 10.0D+0)
00142
00143
00144
00145 CALL DLAHILB(N, N, A, LDA, INVHILB, LDA, B, LDA, WORK, INFO)
00146
00147
00148 CALL DLACPY('ALL', N, N, A, NMAX, ACOPY, NMAX)
00149
00150
00151 DO J = 1, N
00152 DO I = 1, KL+KU+1
00153 AB( I, J ) = 0.0D+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
00163 DO J = 1, N
00164 DO I = 1, KL+KU+1
00165 ABCOPY( I, J ) = 0.0D+0
00166 END DO
00167 END DO
00168 CALL DLACPY('ALL', KL+KU+1, N, AB, LDAB, ABCOPY, LDAB)
00169
00170
00171 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
00172 CALL DSYSVXX(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 DPOSVXX(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 DGBSVXX(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 DGESVXX(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
00197
00198
00199
00200
00201
00202 ELSE
00203
00204
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
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
00219 RNORM = 0.0D+0
00220 RINORM = 0.0D+0
00221 IF ( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'SY' ) ) THEN
00222 DO I = 1, N
00223 SUMR = 0.0D+0
00224 SUMRI = 0.0D+0
00225 DO J = 1, N
00226 SUMR = SUMR + S(I) * ABS(A(I,J)) * S(J)
00227 SUMRI = SUMRI + ABS(INVHILB(I, J)) / (S(J) * S(I))
00228
00229 END DO
00230 RNORM = MAX(RNORM,SUMR)
00231 RINORM = MAX(RINORM,SUMRI)
00232 END DO
00233 ELSE IF ( LSAMEN( 2, C2, 'GE' ) .OR. LSAMEN( 2, C2, 'GB' ) )
00234 $ THEN
00235 DO I = 1, N
00236 SUMR = 0.0D+0
00237 SUMRI = 0.0D+0
00238 DO J = 1, N
00239 SUMR = SUMR + R(I) * ABS(A(I,J)) * C(J)
00240 SUMRI = SUMRI + ABS(INVHILB(I, J)) / (R(J) * C(I))
00241 END DO
00242 RNORM = MAX(RNORM,SUMR)
00243 RINORM = MAX(RINORM,SUMRI)
00244 END DO
00245 END IF
00246
00247 RNORM = RNORM / ABS(A(1, 1))
00248 RCOND = 1.0D+0/(RNORM * RINORM)
00249
00250
00251 DO I = 1, N
00252 RINV(I) = 0.0D+0
00253 END DO
00254 DO J = 1, N
00255 DO I = 1, N
00256 RINV(I) = RINV(I) + ABS(A(I,J))
00257 END DO
00258 END DO
00259
00260
00261 RINORM = 0.0D+0
00262 DO I = 1, N
00263 SUMRI = 0.0D+0
00264 DO J = 1, N
00265 SUMRI = SUMRI + ABS(INVHILB(I,J) * RINV(J))
00266 END DO
00267 RINORM = MAX(RINORM, SUMRI)
00268 END DO
00269
00270
00271
00272 NCOND = ABS(A(1,1)) / RINORM
00273
00274 CONDTHRESH = M * EPS
00275 ERRTHRESH = M * EPS
00276
00277 DO K = 1, NRHS
00278 NORMT = 0.0D+0
00279 NORMDIF = 0.0D+0
00280 CWISE_ERR = 0.0D+0
00281 DO I = 1, N
00282 NORMT = MAX(ABS(INVHILB(I, K)), NORMT)
00283 NORMDIF = MAX(ABS(X(I,K) - INVHILB(I,K)), NORMDIF)
00284 IF (INVHILB(I,K) .NE. 0.0D+0) THEN
00285 CWISE_ERR = MAX(ABS(X(I,K) - INVHILB(I,K))
00286 $ /ABS(INVHILB(I,K)), CWISE_ERR)
00287 ELSE IF (X(I, K) .NE. 0.0D+0) THEN
00288 CWISE_ERR = DLAMCH('OVERFLOW')
00289 END IF
00290 END DO
00291 IF (NORMT .NE. 0.0D+0) THEN
00292 NWISE_ERR = NORMDIF / NORMT
00293 ELSE IF (NORMDIF .NE. 0.0D+0) THEN
00294 NWISE_ERR = DLAMCH('OVERFLOW')
00295 ELSE
00296 NWISE_ERR = 0.0D+0
00297 ENDIF
00298
00299 DO I = 1, N
00300 RINV(I) = 0.0D+0
00301 END DO
00302 DO J = 1, N
00303 DO I = 1, N
00304 RINV(I) = RINV(I) + ABS(A(I, J) * INVHILB(J, K))
00305 END DO
00306 END DO
00307 RINORM = 0.0D+0
00308 DO I = 1, N
00309 SUMRI = 0.0D+0
00310 DO J = 1, N
00311 SUMRI = SUMRI
00312 $ + ABS(INVHILB(I, J) * RINV(J) / INVHILB(I, K))
00313 END DO
00314 RINORM = MAX(RINORM, SUMRI)
00315 END DO
00316
00317
00318 CCOND = ABS(A(1,1))/RINORM
00319
00320
00321 NWISE_BND = ERRBND_N(K + (BND_I-1)*NRHS)
00322 CWISE_BND = ERRBND_C(K + (BND_I-1)*NRHS)
00323 NWISE_RCOND = ERRBND_N(K + (COND_I-1)*NRHS)
00324 CWISE_RCOND = ERRBND_C(K + (COND_I-1)*NRHS)
00325
00326
00327
00328 IF (NCOND .GE. CONDTHRESH) THEN
00329 NGUAR = 'YES'
00330 IF (NWISE_BND .GT. ERRTHRESH) THEN
00331 TSTRAT(1) = 1/(2.0D+0*EPS)
00332 ELSE
00333 IF (NWISE_BND .NE. 0.0D+0) THEN
00334 TSTRAT(1) = NWISE_ERR / NWISE_BND
00335 ELSE IF (NWISE_ERR .NE. 0.0D+0) THEN
00336 TSTRAT(1) = 1/(16.0*EPS)
00337 ELSE
00338 TSTRAT(1) = 0.0D+0
00339 END IF
00340 IF (TSTRAT(1) .GT. 1.0D+0) THEN
00341 TSTRAT(1) = 1/(4.0D+0*EPS)
00342 END IF
00343 END IF
00344 ELSE
00345 NGUAR = 'NO'
00346 IF (NWISE_BND .LT. 1.0D+0) THEN
00347 TSTRAT(1) = 1/(8.0D+0*EPS)
00348 ELSE
00349 TSTRAT(1) = 1.0D+0
00350 END IF
00351 END IF
00352
00353
00354
00355 IF (CCOND .GE. CONDTHRESH) THEN
00356 CGUAR = 'YES'
00357 IF (CWISE_BND .GT. ERRTHRESH) THEN
00358 TSTRAT(2) = 1/(2.0D+0*EPS)
00359 ELSE
00360 IF (CWISE_BND .NE. 0.0D+0) THEN
00361 TSTRAT(2) = CWISE_ERR / CWISE_BND
00362 ELSE IF (CWISE_ERR .NE. 0.0D+0) THEN
00363 TSTRAT(2) = 1/(16.0D+0*EPS)
00364 ELSE
00365 TSTRAT(2) = 0.0D+0
00366 END IF
00367 IF (TSTRAT(2) .GT. 1.0D+0) TSTRAT(2) = 1/(4.0D+0*EPS)
00368 END IF
00369 ELSE
00370 CGUAR = 'NO'
00371 IF (CWISE_BND .LT. 1.0D+0) THEN
00372 TSTRAT(2) = 1/(8.0D+0*EPS)
00373 ELSE
00374 TSTRAT(2) = 1.0D+0
00375 END IF
00376 END IF
00377
00378
00379 TSTRAT(3) = BERR(K)/EPS
00380
00381
00382 TSTRAT(4) = RCOND / ORCOND
00383 IF (RCOND .GE. CONDTHRESH .AND. TSTRAT(4) .LT. 1.0D+0)
00384 $ TSTRAT(4) = 1.0D+0 / TSTRAT(4)
00385
00386 TSTRAT(5) = NCOND / NWISE_RCOND
00387 IF (NCOND .GE. CONDTHRESH .AND. TSTRAT(5) .LT. 1.0D+0)
00388 $ TSTRAT(5) = 1.0D+0 / TSTRAT(5)
00389
00390 TSTRAT(6) = CCOND / NWISE_RCOND
00391 IF (CCOND .GE. CONDTHRESH .AND. TSTRAT(6) .LT. 1.0D+0)
00392 $ TSTRAT(6) = 1.0D+0 / TSTRAT(6)
00393
00394 DO I = 1, NTESTS
00395 IF (TSTRAT(I) .GT. THRESH) THEN
00396 IF (.NOT.PRINTED_GUIDE) THEN
00397 WRITE(*,*)
00398 WRITE( *, 9996) 1
00399 WRITE( *, 9995) 2
00400 WRITE( *, 9994) 3
00401 WRITE( *, 9993) 4
00402 WRITE( *, 9992) 5
00403 WRITE( *, 9991) 6
00404 WRITE( *, 9990) 7
00405 WRITE( *, 9989) 8
00406 WRITE(*,*)
00407 PRINTED_GUIDE = .TRUE.
00408 END IF
00409 WRITE( *, 9999) C2, N, K, NGUAR, CGUAR, I, TSTRAT(I)
00410 NFAIL = NFAIL + 1
00411 END IF
00412 END DO
00413 END DO
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429 END DO
00430
00431 WRITE(*,*)
00432 IF( NFAIL .GT. 0 ) THEN
00433 WRITE(*,9998) C2, NFAIL, NTESTS*N+N_AUX_TESTS
00434 ELSE
00435 WRITE(*,9997) C2
00436 END IF
00437 9999 FORMAT( ' D', A2, 'SVXX: N =', I2, ', RHS = ', I2,
00438 $ ', NWISE GUAR. = ', A, ', CWISE GUAR. = ', A,
00439 $ ' test(',I1,') =', G12.5 )
00440 9998 FORMAT( ' D', A2, 'SVXX: ', I6, ' out of ', I6,
00441 $ ' tests failed to pass the threshold' )
00442 9997 FORMAT( ' D', A2, 'SVXX passed the tests of error bounds' )
00443
00444 9996 FORMAT( 3X, I2, ': Normwise guaranteed forward error', / 5X,
00445 $ 'Guaranteed case: if norm ( abs( Xc - Xt )',
00446 $ ' / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then',
00447 $ / 5X,
00448 $ 'ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS')
00449 9995 FORMAT( 3X, I2, ': Componentwise guaranteed forward error' )
00450 9994 FORMAT( 3X, I2, ': Backwards error' )
00451 9993 FORMAT( 3X, I2, ': Reciprocal condition number' )
00452 9992 FORMAT( 3X, I2, ': Reciprocal normwise condition number' )
00453 9991 FORMAT( 3X, I2, ': Raw normwise error estimate' )
00454 9990 FORMAT( 3X, I2, ': Reciprocal componentwise condition number' )
00455 9989 FORMAT( 3X, I2, ': Raw componentwise error estimate' )
00456
00457 8000 FORMAT( ' D', A2, 'SVXX: N =', I2, ', INFO = ', I3,
00458 $ ', ORCOND = ', G12.5, ', real RCOND = ', G12.5 )
00459
00460 END