98 DOUBLE PRECISION THRESH
101 INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU
102 parameter(nmax = 10, nparams = 2, nerrbnd = 3,
106 INTEGER N, NRHS, INFO, I ,J, k, NFAIL, LDA,
107 $ N_AUX_TESTS, LDAB, LDAFB
108 CHARACTER FACT, TRANS, UPLO, EQUED
110 CHARACTER(3) NGUAR, CGUAR
111 LOGICAL printed_guide
112 DOUBLE PRECISION 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
120 DOUBLE PRECISION TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
121 $ S(NMAX),R(NMAX),C(NMAX), DIFF(NMAX, NMAX),
122 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3),
123 $ A(NMAX,NMAX),INVHILB(NMAX,NMAX),X(NMAX,NMAX),
124 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
125 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
126 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX ),
127 $ WORK(NMAX*3*5), AF(NMAX, NMAX),B(NMAX, NMAX),
129 INTEGER IPIV(NMAX), IWORK(3*NMAX)
132 DOUBLE PRECISION DLAMCH
140 INTRINSIC sqrt, max, abs, dble
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)
154 eps = dlamch(
'Epsilon')
158 ldab = (nmax-1)+(nmax-1)+1
159 ldafb = 2*(nmax-1)+(nmax-1)+1
164 printed_guide = .false.
173 m = max(sqrt(dble(n)), 10.0d+0)
177 CALL dlahilb(n, n, a, lda, invhilb, lda, b, lda, work, info)
180 CALL dlacpy(
'ALL', n, n, a, nmax, acopy, nmax)
189 DO i = max( 1, j-ku ), min( n, j+kl )
190 ab( ku+1+i-j, j ) = a( i, j )
197 abcopy( i, j ) = 0.0d+0
200 CALL dlacpy(
'ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)
203 IF ( lsamen( 2, c2,
'SY' ) )
THEN
204 CALL dsysvxx(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 dposvxx(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 dgbsvxx(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,
220 CALL dgesvxx(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)
226 n_aux_tests = n_aux_tests + 1
227 IF (orcond .LT. eps)
THEN
237 IF (info .GT. 0 .AND. info .LE. n+1)
THEN
239 WRITE (*, fmt=8000) c2, n, info, orcond, rcond
246 diff(i,j) = x(i,j) - invhilb(i,j)
253 IF ( lsamen( 2, c2,
'PO' ) .OR. lsamen( 2, c2,
'SY' ) )
THEN
258 sumr = sumr + s(i) * abs(a(i,j)) * s(j)
259 sumri = sumri + abs(invhilb(i, j)) / (s(j) * s(i))
262 rnorm = max(rnorm,sumr)
263 rinorm = max(rinorm,sumri)
265 ELSE IF ( lsamen( 2, c2,
'GE' ) .OR. lsamen( 2, c2,
'GB' ) )
271 sumr = sumr + r(i) * abs(a(i,j)) * c(j)
272 sumri = sumri + abs(invhilb(i, j)) / (r(j) * c(i))
274 rnorm = max(rnorm,sumr)
275 rinorm = max(rinorm,sumri)
279 rnorm = rnorm / abs(a(1, 1))
280 rcond = 1.0d+0/(rnorm * rinorm)
288 rinv(i) = rinv(i) + abs(a(i,j))
297 sumri = sumri + abs(invhilb(i,j) * rinv(j))
299 rinorm = max(rinorm, sumri)
304 ncond = abs(a(1,1)) / rinorm
314 normt = max(abs(invhilb(i, k)), normt)
315 normdif = max(abs(x(i,k) - invhilb(i,k)), normdif)
316 IF (invhilb(i,k) .NE. 0.0d+0)
THEN
317 cwise_err = max(abs(x(i,k) - invhilb(i,k))
318 $ /abs(invhilb(i,k)), cwise_err)
319 ELSE IF (x(i, k) .NE. 0.0d+0)
THEN
320 cwise_err = dlamch(
'OVERFLOW')
323 IF (normt .NE. 0.0d+0)
THEN
324 nwise_err = normdif / normt
325 ELSE IF (normdif .NE. 0.0d+0)
THEN
326 nwise_err = dlamch(
'OVERFLOW')
336 rinv(i) = rinv(i) + abs(a(i, j) * invhilb(j, k))
344 $ + abs(invhilb(i, j) * rinv(j) / invhilb(i, k))
346 rinorm = max(rinorm, sumri)
350 ccond = abs(a(1,1))/rinorm
353 nwise_bnd = errbnd_n(k + (bnd_i-1)*nrhs)
354 cwise_bnd = errbnd_c(k + (bnd_i-1)*nrhs)
355 nwise_rcond = errbnd_n(k + (cond_i-1)*nrhs)
356 cwise_rcond = errbnd_c(k + (cond_i-1)*nrhs)
360 IF (ncond .GE. condthresh)
THEN
362 IF (nwise_bnd .GT. errthresh)
THEN
363 tstrat(1) = 1/(2.0d+0*eps)
365 IF (nwise_bnd .NE. 0.0d+0)
THEN
366 tstrat(1) = nwise_err / nwise_bnd
367 ELSE IF (nwise_err .NE. 0.0d+0)
THEN
368 tstrat(1) = 1/(16.0*eps)
372 IF (tstrat(1) .GT. 1.0d+0)
THEN
373 tstrat(1) = 1/(4.0d+0*eps)
378 IF (nwise_bnd .LT. 1.0d+0)
THEN
379 tstrat(1) = 1/(8.0d+0*eps)
387 IF (ccond .GE. condthresh)
THEN
389 IF (cwise_bnd .GT. errthresh)
THEN
390 tstrat(2) = 1/(2.0d+0*eps)
392 IF (cwise_bnd .NE. 0.0d+0)
THEN
393 tstrat(2) = cwise_err / cwise_bnd
394 ELSE IF (cwise_err .NE. 0.0d+0)
THEN
395 tstrat(2) = 1/(16.0d+0*eps)
399 IF (tstrat(2) .GT. 1.0d+0) tstrat(2) = 1/(4.0d+0*eps)
403 IF (cwise_bnd .LT. 1.0d+0)
THEN
404 tstrat(2) = 1/(8.0d+0*eps)
411 tstrat(3) = berr(k)/eps
414 tstrat(4) = rcond / orcond
415 IF (rcond .GE. condthresh .AND. tstrat(4) .LT. 1.0d+0)
416 $ tstrat(4) = 1.0d+0 / tstrat(4)
418 tstrat(5) = ncond / nwise_rcond
419 IF (ncond .GE. condthresh .AND. tstrat(5) .LT. 1.0d+0)
420 $ tstrat(5) = 1.0d+0 / tstrat(5)
422 tstrat(6) = ccond / nwise_rcond
423 IF (ccond .GE. condthresh .AND. tstrat(6) .LT. 1.0d+0)
424 $ tstrat(6) = 1.0d+0 / tstrat(6)
427 IF (tstrat(i) .GT. thresh)
THEN
428 IF (.NOT.printed_guide)
THEN
439 printed_guide = .true.
441 WRITE( *, 9999) c2, n, k, nguar, cguar, i, tstrat(i)
464 IF( nfail .GT. 0 )
THEN
465 WRITE(*,9998) c2, nfail, ntests*n+n_aux_tests
469 9999
FORMAT(
' D', a2,
'SVXX: N =', i2,
', RHS = ', i2,
470 $
', NWISE GUAR. = ', a,
', CWISE GUAR. = ', a,
471 $
' test(',i1,
') =', g12.5 )
472 9998
FORMAT(
' D', a2,
'SVXX: ', i6,
' out of ', i6,
473 $
' tests failed to pass the threshold' )
474 9997
FORMAT(
' D', a2,
'SVXX passed the tests of error bounds' )
476 9996
FORMAT( 3x, i2,
': Normwise guaranteed forward error', / 5x,
477 $
'Guaranteed case: if norm ( abs( Xc - Xt )',
478 $ .LE.
' / norm ( Xt ) ERRBND( *, nwise_i, bnd_i ), then',
480 $ .LE.
'ERRBND( *, nwise_i, bnd_i ) MAX(SQRT(N), 10) * EPS')
481 9995
FORMAT( 3x, i2,
': Componentwise guaranteed forward error' )
482 9994
FORMAT( 3x, i2,
': Backwards error' )
483 9993
FORMAT( 3x, i2,
': Reciprocal condition number' )
484 9992
FORMAT( 3x, i2,
': Reciprocal normwise condition number' )
485 9991
FORMAT( 3x, i2,
': Raw normwise error estimate' )
486 9990
FORMAT( 3x, i2,
': Reciprocal componentwise condition number' )
487 9989
FORMAT( 3x, i2,
': Raw componentwise error estimate' )
489 8000
FORMAT(
' D', a2,
'SVXX: N =', i2,
', INFO = ', i3,
490 $
', ORCOND = ', g12.5,
', real RCOND = ', g12.5 )