00001 SUBROUTINE CEBCHVXX( THRESH, PATH )
00002 IMPLICIT NONE
00003
00004 REAL 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 = 6, 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 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 COMPLEX ZDUM
00087
00088
00089 REAL TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
00090 $ S(NMAX), R(NMAX),C(NMAX),RWORK(3*NMAX),
00091 $ DIFF(NMAX, NMAX),
00092 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3)
00093 INTEGER IPIV(NMAX)
00094 COMPLEX A(NMAX,NMAX),INVHILB(NMAX,NMAX),X(NMAX,NMAX),
00095 $ WORK(NMAX*3*5), AF(NMAX, NMAX),B(NMAX, NMAX),
00096 $ ACOPY(NMAX, NMAX),
00097 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
00098 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
00099 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX )
00100
00101
00102 REAL SLAMCH
00103
00104
00105 EXTERNAL CLAHILB, CGESVXX, CSYSVXX, CPOSVXX,
00106 $ CGBSVXX, CLACPY, LSAMEN
00107 LOGICAL LSAMEN
00108
00109
00110 INTRINSIC SQRT, MAX, ABS, REAL, AIMAG
00111
00112
00113 REAL CABS1
00114
00115
00116 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00117
00118
00119 INTEGER NWISE_I, CWISE_I
00120 PARAMETER (NWISE_I = 1, CWISE_I = 1)
00121 INTEGER BND_I, COND_I
00122 PARAMETER (BND_I = 2, COND_I = 3)
00123
00124
00125
00126 FACT = 'E'
00127 UPLO = 'U'
00128 TRANS = 'N'
00129 EQUED = 'N'
00130 EPS = SLAMCH('Epsilon')
00131 NFAIL = 0
00132 N_AUX_TESTS = 0
00133 LDA = NMAX
00134 LDAB = (NMAX-1)+(NMAX-1)+1
00135 LDAFB = 2*(NMAX-1)+(NMAX-1)+1
00136 C2 = PATH( 2: 3 )
00137
00138
00139
00140 printed_guide = .false.
00141
00142 DO N = 1 , NMAX
00143 PARAMS(1) = -1
00144 PARAMS(2) = -1
00145
00146 KL = N-1
00147 KU = N-1
00148 NRHS = n
00149 M = MAX(SQRT(REAL(N)), 10.0)
00150
00151
00152
00153 CALL CLAHILB(N, N, A, LDA, INVHILB, LDA, B,
00154 $ LDA, WORK, INFO, PATH)
00155
00156
00157 CALL CLACPY('ALL', N, N, A, NMAX, ACOPY, NMAX)
00158
00159
00160 DO J = 1, N
00161 DO I = 1, KL+KU+1
00162 AB( I, J ) = (0.0E+0,0.0E+0)
00163 END DO
00164 END DO
00165 DO J = 1, N
00166 DO I = MAX( 1, J-KU ), MIN( N, J+KL )
00167 AB( KU+1+I-J, J ) = A( I, J )
00168 END DO
00169 END DO
00170
00171
00172 DO J = 1, N
00173 DO I = 1, KL+KU+1
00174 ABCOPY( I, J ) = (0.0E+0,0.0E+0)
00175 END DO
00176 END DO
00177 CALL CLACPY('ALL', KL+KU+1, N, AB, LDAB, ABCOPY, LDAB)
00178
00179
00180 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
00181 CALL CSYSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
00182 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
00183 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
00184 $ PARAMS, WORK, RWORK, INFO)
00185 ELSE IF ( LSAMEN( 2, C2, 'PO' ) ) THEN
00186 CALL CPOSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
00187 $ EQUED, S, B, LDA, X, LDA, ORCOND,
00188 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
00189 $ PARAMS, WORK, RWORK, INFO)
00190 ELSE IF ( LSAMEN( 2, C2, 'HE' ) ) THEN
00191 CALL CHESVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
00192 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
00193 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
00194 $ PARAMS, WORK, RWORK, INFO)
00195 ELSE IF ( LSAMEN( 2, C2, 'GB' ) ) THEN
00196 CALL CGBSVXX(FACT, TRANS, N, KL, KU, NRHS, ABCOPY,
00197 $ LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B,
00198 $ LDA, X, LDA, ORCOND, RPVGRW, BERR, NERRBND,
00199 $ ERRBND_N, ERRBND_C, NPARAMS, PARAMS, WORK, RWORK,
00200 $ INFO)
00201 ELSE
00202 CALL CGESVXX(FACT, TRANS, N, NRHS, ACOPY, LDA, AF, LDA,
00203 $ IPIV, EQUED, R, C, B, LDA, X, LDA, ORCOND,
00204 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
00205 $ PARAMS, WORK, RWORK, INFO)
00206 END IF
00207
00208 N_AUX_TESTS = N_AUX_TESTS + 1
00209 IF (ORCOND .LT. EPS) THEN
00210
00211
00212
00213
00214
00215
00216 ELSE
00217
00218
00219 IF (INFO .GT. 0 .AND. INFO .LE. N+1) THEN
00220 NFAIL = NFAIL + 1
00221 WRITE (*, FMT=8000) C2, N, INFO, ORCOND, RCOND
00222 END IF
00223 END IF
00224
00225
00226 DO I = 1,N
00227 DO J =1,NRHS
00228 DIFF(I,J) = X(I,J) - INVHILB(I,J)
00229 END DO
00230 END DO
00231
00232
00233 RNORM = 0
00234 RINORM = 0
00235 IF ( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'SY' ) .OR.
00236 $ LSAMEN( 2, C2, 'HE' ) ) THEN
00237 DO I = 1, N
00238 SUMR = 0
00239 SUMRI = 0
00240 DO J = 1, N
00241 SUMR = SUMR + S(I) * CABS1(A(I,J)) * S(J)
00242 SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (S(J) * S(I))
00243 END DO
00244 RNORM = MAX(RNORM,SUMR)
00245 RINORM = MAX(RINORM,SUMRI)
00246 END DO
00247 ELSE IF ( LSAMEN( 2, C2, 'GE' ) .OR. LSAMEN( 2, C2, 'GB' ) )
00248 $ THEN
00249 DO I = 1, N
00250 SUMR = 0
00251 SUMRI = 0
00252 DO J = 1, N
00253 SUMR = SUMR + R(I) * CABS1(A(I,J)) * C(J)
00254 SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (R(J) * C(I))
00255 END DO
00256 RNORM = MAX(RNORM,SUMR)
00257 RINORM = MAX(RINORM,SUMRI)
00258 END DO
00259 END IF
00260
00261 RNORM = RNORM / CABS1(A(1, 1))
00262 RCOND = 1.0/(RNORM * RINORM)
00263
00264
00265 DO I = 1, N
00266 RINV(I) = 0.0
00267 END DO
00268 DO J = 1, N
00269 DO I = 1, N
00270 RINV(I) = RINV(I) + CABS1(A(I,J))
00271 END DO
00272 END DO
00273
00274
00275 RINORM = 0.0
00276 DO I = 1, N
00277 SUMRI = 0.0
00278 DO J = 1, N
00279 SUMRI = SUMRI + CABS1(INVHILB(I,J) * RINV(J))
00280 END DO
00281 RINORM = MAX(RINORM, SUMRI)
00282 END DO
00283
00284
00285
00286 NCOND = CABS1(A(1,1)) / RINORM
00287
00288 CONDTHRESH = M * EPS
00289 ERRTHRESH = M * EPS
00290
00291 DO K = 1, NRHS
00292 NORMT = 0.0
00293 NORMDIF = 0.0
00294 CWISE_ERR = 0.0
00295 DO I = 1, N
00296 NORMT = MAX(CABS1(INVHILB(I, K)), NORMT)
00297 NORMDIF = MAX(CABS1(X(I,K) - INVHILB(I,K)), NORMDIF)
00298 IF (INVHILB(I,K) .NE. 0.0) THEN
00299 CWISE_ERR = MAX(CABS1(X(I,K) - INVHILB(I,K))
00300 $ /CABS1(INVHILB(I,K)), CWISE_ERR)
00301 ELSE IF (X(I, K) .NE. 0.0) THEN
00302 CWISE_ERR = SLAMCH('OVERFLOW')
00303 END IF
00304 END DO
00305 IF (NORMT .NE. 0.0) THEN
00306 NWISE_ERR = NORMDIF / NORMT
00307 ELSE IF (NORMDIF .NE. 0.0) THEN
00308 NWISE_ERR = SLAMCH('OVERFLOW')
00309 ELSE
00310 NWISE_ERR = 0.0
00311 ENDIF
00312
00313 DO I = 1, N
00314 RINV(I) = 0.0
00315 END DO
00316 DO J = 1, N
00317 DO I = 1, N
00318 RINV(I) = RINV(I) + CABS1(A(I, J) * INVHILB(J, K))
00319 END DO
00320 END DO
00321 RINORM = 0.0
00322 DO I = 1, N
00323 SUMRI = 0.0
00324 DO J = 1, N
00325 SUMRI = SUMRI
00326 $ + CABS1(INVHILB(I, J) * RINV(J) / INVHILB(I, K))
00327 END DO
00328 RINORM = MAX(RINORM, SUMRI)
00329 END DO
00330
00331
00332 CCOND = CABS1(A(1,1))/RINORM
00333
00334
00335 NWISE_BND = ERRBND_N(K + (BND_I-1)*NRHS)
00336 CWISE_BND = ERRBND_C(K + (BND_I-1)*NRHS)
00337 NWISE_RCOND = ERRBND_N(K + (COND_I-1)*NRHS)
00338 CWISE_RCOND = ERRBND_C(K + (COND_I-1)*NRHS)
00339
00340
00341
00342 IF (NCOND .GE. CONDTHRESH) THEN
00343 NGUAR = 'YES'
00344 IF (NWISE_BND .GT. ERRTHRESH) THEN
00345 TSTRAT(1) = 1/(2.0*EPS)
00346 ELSE
00347 IF (NWISE_BND .NE. 0.0) THEN
00348 TSTRAT(1) = NWISE_ERR / NWISE_BND
00349 ELSE IF (NWISE_ERR .NE. 0.0) THEN
00350 TSTRAT(1) = 1/(16.0*EPS)
00351 ELSE
00352 TSTRAT(1) = 0.0
00353 END IF
00354 IF (TSTRAT(1) .GT. 1.0) THEN
00355 TSTRAT(1) = 1/(4.0*EPS)
00356 END IF
00357 END IF
00358 ELSE
00359 NGUAR = 'NO'
00360 IF (NWISE_BND .LT. 1.0) THEN
00361 TSTRAT(1) = 1/(8.0*EPS)
00362 ELSE
00363 TSTRAT(1) = 1.0
00364 END IF
00365 END IF
00366
00367
00368
00369 IF (CCOND .GE. CONDTHRESH) THEN
00370 CGUAR = 'YES'
00371 IF (CWISE_BND .GT. ERRTHRESH) THEN
00372 TSTRAT(2) = 1/(2.0*EPS)
00373 ELSE
00374 IF (CWISE_BND .NE. 0.0) THEN
00375 TSTRAT(2) = CWISE_ERR / CWISE_BND
00376 ELSE IF (CWISE_ERR .NE. 0.0) THEN
00377 TSTRAT(2) = 1/(16.0*EPS)
00378 ELSE
00379 TSTRAT(2) = 0.0
00380 END IF
00381 IF (TSTRAT(2) .GT. 1.0) TSTRAT(2) = 1/(4.0*EPS)
00382 END IF
00383 ELSE
00384 CGUAR = 'NO'
00385 IF (CWISE_BND .LT. 1.0) THEN
00386 TSTRAT(2) = 1/(8.0*EPS)
00387 ELSE
00388 TSTRAT(2) = 1.0
00389 END IF
00390 END IF
00391
00392
00393 TSTRAT(3) = BERR(K)/EPS
00394
00395
00396 TSTRAT(4) = RCOND / ORCOND
00397 IF (RCOND .GE. CONDTHRESH .AND. TSTRAT(4) .LT. 1.0)
00398 $ TSTRAT(4) = 1.0 / TSTRAT(4)
00399
00400 TSTRAT(5) = NCOND / NWISE_RCOND
00401 IF (NCOND .GE. CONDTHRESH .AND. TSTRAT(5) .LT. 1.0)
00402 $ TSTRAT(5) = 1.0 / TSTRAT(5)
00403
00404 TSTRAT(6) = CCOND / NWISE_RCOND
00405 IF (CCOND .GE. CONDTHRESH .AND. TSTRAT(6) .LT. 1.0)
00406 $ TSTRAT(6) = 1.0 / TSTRAT(6)
00407
00408 DO I = 1, NTESTS
00409 IF (TSTRAT(I) .GT. THRESH) THEN
00410 IF (.NOT.PRINTED_GUIDE) THEN
00411 WRITE(*,*)
00412 WRITE( *, 9996) 1
00413 WRITE( *, 9995) 2
00414 WRITE( *, 9994) 3
00415 WRITE( *, 9993) 4
00416 WRITE( *, 9992) 5
00417 WRITE( *, 9991) 6
00418 WRITE( *, 9990) 7
00419 WRITE( *, 9989) 8
00420 WRITE(*,*)
00421 PRINTED_GUIDE = .TRUE.
00422 END IF
00423 WRITE( *, 9999) C2, N, K, NGUAR, CGUAR, I, TSTRAT(I)
00424 NFAIL = NFAIL + 1
00425 END IF
00426 END DO
00427 END DO
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443 END DO
00444
00445 WRITE(*,*)
00446 IF( NFAIL .GT. 0 ) THEN
00447 WRITE(*,9998) C2, NFAIL, NTESTS*N+N_AUX_TESTS
00448 ELSE
00449 WRITE(*,9997) C2
00450 END IF
00451 9999 FORMAT( ' C', A2, 'SVXX: N =', I2, ', RHS = ', I2,
00452 $ ', NWISE GUAR. = ', A, ', CWISE GUAR. = ', A,
00453 $ ' test(',I1,') =', G12.5 )
00454 9998 FORMAT( ' C', A2, 'SVXX: ', I6, ' out of ', I6,
00455 $ ' tests failed to pass the threshold' )
00456 9997 FORMAT( ' C', A2, 'SVXX passed the tests of error bounds' )
00457
00458 9996 FORMAT( 3X, I2, ': Normwise guaranteed forward error', / 5X,
00459 $ 'Guaranteed case: if norm ( abs( Xc - Xt )',
00460 $ ' / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then',
00461 $ / 5X,
00462 $ 'ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS')
00463 9995 FORMAT( 3X, I2, ': Componentwise guaranteed forward error' )
00464 9994 FORMAT( 3X, I2, ': Backwards error' )
00465 9993 FORMAT( 3X, I2, ': Reciprocal condition number' )
00466 9992 FORMAT( 3X, I2, ': Reciprocal normwise condition number' )
00467 9991 FORMAT( 3X, I2, ': Raw normwise error estimate' )
00468 9990 FORMAT( 3X, I2, ': Reciprocal componentwise condition number' )
00469 9989 FORMAT( 3X, I2, ': Raw componentwise error estimate' )
00470
00471 8000 FORMAT( ' C', A2, 'SVXX: N =', I2, ', INFO = ', I3,
00472 $ ', ORCOND = ', G12.5, ', real RCOND = ', G12.5 )
00473
00474 END