101 INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU
102 parameter(nmax = 6, nparams = 2, nerrbnd = 3,
106 INTEGER N, NRHS, INFO, I ,J, k, NFAIL, LDA, LDAB,
108 CHARACTER FACT, TRANS, UPLO, EQUED
110 CHARACTER(3) NGUAR, CGUAR
111 LOGICAL printed_guide
112 REAL 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 REAL TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
121 $ A(NMAX, NMAX), ACOPY(NMAX, NMAX),
122 $ INVHILB(NMAX, NMAX), R(NMAX), C(NMAX), S(NMAX),
123 $ WORK(NMAX * 5), B(NMAX, NMAX), X(NMAX, NMAX),
124 $ DIFF(NMAX, NMAX), AF(NMAX, NMAX),
125 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
126 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
127 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX ),
128 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3)
129 INTEGER IWORK(NMAX), IPIV(NMAX)
140 INTRINSIC sqrt, max, abs
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 = slamch(
'Epsilon')
158 ldab = (nmax-1)+(nmax-1)+1
159 ldafb = 2*(nmax-1)+(nmax-1)+1
164 printed_guide = .false.
173 m = max(sqrt(real(n)), 10.0)
177 CALL slahilb(n, n, a, lda, invhilb, lda, b, lda, work, info)
180 CALL slacpy(
'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.0e+0
200 CALL slacpy(
'ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)
203 IF ( lsamen( 2, c2,
'SY' ) )
THEN
204 CALL ssysvxx(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 sposvxx(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 sgbsvxx(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 sgesvxx(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 + abs(s(i) * a(i,j) * s(j))
259 sumri = sumri + abs(invhilb(i, j) / s(j) / s(i))
261 rnorm = max(rnorm,sumr)
262 rinorm = max(rinorm,sumri)
264 ELSE IF ( lsamen( 2, c2,
'GE' ) .OR. lsamen( 2, c2,
'GB' ) )
270 sumr = sumr + abs(r(i) * a(i,j) * c(j))
271 sumri = sumri + abs(invhilb(i, j) / r(j) / c(i))
273 rnorm = max(rnorm,sumr)
274 rinorm = max(rinorm,sumri)
278 rnorm = rnorm / a(1, 1)
279 rcond = 1.0/(rnorm * rinorm)
287 rinv(i) = rinv(i) + abs(a(i,j))
296 sumri = sumri + abs(invhilb(i,j) * rinv(j))
298 rinorm = max(rinorm, sumri)
303 ncond = a(1,1) / rinorm
313 normt = max(abs(invhilb(i, k)), normt)
314 normdif = max(abs(x(i,k) - invhilb(i,k)), normdif)
315 IF (invhilb(i,k) .NE. 0.0)
THEN
316 cwise_err = max(abs(x(i,k) - invhilb(i,k))
317 $ /abs(invhilb(i,k)), cwise_err)
318 ELSE IF (x(i, k) .NE. 0.0)
THEN
319 cwise_err = slamch(
'OVERFLOW')
322 IF (normt .NE. 0.0)
THEN
323 nwise_err = normdif / normt
324 ELSE IF (normdif .NE. 0.0)
THEN
325 nwise_err = slamch(
'OVERFLOW')
335 rinv(i) = rinv(i) + abs(a(i, j) * invhilb(j, k))
343 $ + abs(invhilb(i, j) * rinv(j) / invhilb(i, k))
345 rinorm = max(rinorm, sumri)
349 ccond = a(1,1)/rinorm
352 nwise_bnd = errbnd_n(k + (bnd_i-1)*nrhs)
353 cwise_bnd = errbnd_c(k + (bnd_i-1)*nrhs)
354 nwise_rcond = errbnd_n(k + (cond_i-1)*nrhs)
355 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.0*eps)
366 IF (nwise_bnd .NE. 0.0)
THEN
367 tstrat(1) = nwise_err / nwise_bnd
368 ELSE IF (nwise_err .NE. 0.0)
THEN
369 tstrat(1) = 1/(16.0*eps)
373 IF (tstrat(1) .GT. 1.0)
THEN
374 tstrat(1) = 1/(4.0*eps)
379 IF (nwise_bnd .LT. 1.0)
THEN
380 tstrat(1) = 1/(8.0*eps)
388 IF (ccond .GE. condthresh)
THEN
391 IF (cwise_bnd .GT. errthresh)
THEN
392 tstrat(2) = 1/(2.0*eps)
394 IF (cwise_bnd .NE. 0.0)
THEN
395 tstrat(2) = cwise_err / cwise_bnd
396 ELSE IF (cwise_err .NE. 0.0)
THEN
397 tstrat(2) = 1/(16.0*eps)
401 IF (tstrat(2) .GT. 1.0) tstrat(2) = 1/(4.0*eps)
405 IF (cwise_bnd .LT. 1.0)
THEN
406 tstrat(2) = 1/(8.0*eps)
413 tstrat(3) = berr(k)/eps
416 tstrat(4) = rcond / orcond
417 IF (rcond .GE. condthresh .AND. tstrat(4) .LT. 1.0)
418 $ tstrat(4) = 1.0 / tstrat(4)
420 tstrat(5) = ncond / nwise_rcond
421 IF (ncond .GE. condthresh .AND. tstrat(5) .LT. 1.0)
422 $ tstrat(5) = 1.0 / tstrat(5)
424 tstrat(6) = ccond / nwise_rcond
425 IF (ccond .GE. condthresh .AND. tstrat(6) .LT. 1.0)
426 $ tstrat(6) = 1.0 / tstrat(6)
429 IF (tstrat(i) .GT. thresh)
THEN
430 IF (.NOT.printed_guide)
THEN
441 printed_guide = .true.
443 WRITE( *, 9999) c2, n, k, nguar, cguar, i, tstrat(i)
466 IF( nfail .GT. 0 )
THEN
467 WRITE(*,9998) c2, nfail, ntests*n+n_aux_tests
471 9999
FORMAT(
' S', a2,
'SVXX: N =', i2,
', RHS = ', i2,
472 $
', NWISE GUAR. = ', a,
', CWISE GUAR. = ', a,
473 $
' test(',i1,
') =', g12.5 )
474 9998
FORMAT(
' S', a2,
'SVXX: ', i6,
' out of ', i6,
475 $
' tests failed to pass the threshold' )
476 9997
FORMAT(
' S', a2,
'SVXX passed the tests of error bounds' )
478 9996
FORMAT( 3x, i2,
': Normwise guaranteed forward error', / 5x,
479 $
'Guaranteed case: if norm ( abs( Xc - Xt )',
480 $ .LE.
' / norm ( Xt ) ERRBND( *, nwise_i, bnd_i ), then',
482 $ .LE.
'ERRBND( *, nwise_i, bnd_i ) MAX(SQRT(N), 10) * EPS')
483 9995
FORMAT( 3x, i2,
': Componentwise guaranteed forward error' )
484 9994
FORMAT( 3x, i2,
': Backwards error' )
485 9993
FORMAT( 3x, i2,
': Reciprocal condition number' )
486 9992
FORMAT( 3x, i2,
': Reciprocal normwise condition number' )
487 9991
FORMAT( 3x, i2,
': Raw normwise error estimate' )
488 9990
FORMAT( 3x, i2,
': Reciprocal componentwise condition number' )
489 9989
FORMAT( 3x, i2,
': Raw componentwise error estimate' )
491 8000
FORMAT(
' S', a2,
'SVXX: N =', i2,
', INFO = ', i3,
492 $
', ORCOND = ', g12.5,
', real RCOND = ', g12.5 )