100 DOUBLE PRECISION thresh
103 INTEGER nmax, nparams, nerrbnd, ntests, kl, ku
104 parameter (nmax = 10, nparams = 2, nerrbnd = 3,
108 INTEGER n, nrhs, info, i ,j, k, nfail, lda,
109 $ n_aux_tests, ldab, ldafb
110 CHARACTER fact, trans, uplo, equed
112 CHARACTER(3) nguar, cguar
113 LOGICAL printed_guide
114 DOUBLE PRECISION ncond, ccond, m, normdif, normt, rcond,
115 $ rnorm, rinorm, sumr, sumri, eps,
116 $ berr(nmax), rpvgrw, orcond,
117 $ cwise_err, nwise_err, cwise_bnd, nwise_bnd,
118 $ cwise_rcond, nwise_rcond,
119 $ condthresh, errthresh
122 DOUBLE PRECISION tstrat(ntests), rinv(nmax), params(nparams),
123 $ s(nmax),r(nmax),c(nmax), diff(nmax, nmax),
124 $ errbnd_n(nmax*3), errbnd_c(nmax*3),
125 $ a(nmax,nmax),invhilb(nmax,nmax),x(nmax,nmax),
126 $ ab( (nmax-1)+(nmax-1)+1, nmax ),
127 $ abcopy( (nmax-1)+(nmax-1)+1, nmax ),
128 $ afb( 2*(nmax-1)+(nmax-1)+1, nmax ),
129 $ work(nmax*3*5), af(nmax, nmax),b(nmax, nmax),
131 INTEGER ipiv(nmax), iwork(3*nmax)
142 INTRINSIC sqrt, max, abs, dble
145 INTEGER nwise_i, cwise_i
146 parameter (nwise_i = 1, cwise_i = 1)
147 INTEGER bnd_i, cond_i
148 parameter (bnd_i = 2, cond_i = 3)
160 ldab = (nmax-1)+(nmax-1)+1
161 ldafb = 2*(nmax-1)+(nmax-1)+1
166 printed_guide = .false.
175 m = max(sqrt(dble(n)), 10.0d+0)
179 CALL dlahilb(n, n, a, lda, invhilb, lda, b, lda, work, info)
182 CALL dlacpy(
'ALL', n, n, a, nmax, acopy, nmax)
191 DO i = max( 1, j-ku ), min( n, j+kl )
192 ab( ku+1+i-j, j ) = a( i, j )
199 abcopy( i, j ) = 0.0d+0
202 CALL dlacpy(
'ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)
205 IF (
lsamen( 2, c2,
'SY' ) )
THEN
206 CALL dsysvxx(fact, uplo, n, nrhs, acopy, lda, af, lda,
207 $ ipiv, equed, s, b, lda, x, lda, orcond,
208 $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
209 $ params, work, iwork, info)
210 ELSE IF (
lsamen( 2, c2,
'PO' ) )
THEN
211 CALL dposvxx(fact, uplo, n, nrhs, acopy, lda, af, lda,
212 $ equed, s, b, lda, x, lda, orcond,
213 $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
214 $ params, work, iwork, info)
215 ELSE IF (
lsamen( 2, c2,
'GB' ) )
THEN
216 CALL dgbsvxx(fact, trans, n, kl, ku, nrhs, abcopy,
217 $ ldab, afb, ldafb, ipiv, equed, r, c, b,
218 $ lda, x, lda, orcond, rpvgrw, berr, nerrbnd,
219 $ errbnd_n, errbnd_c, nparams, params, work, iwork,
222 CALL dgesvxx(fact, trans, n, nrhs, acopy, lda, af, lda,
223 $ ipiv, equed, r, c, b, lda, x, lda, orcond,
224 $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
225 $ params, work, iwork, info)
228 n_aux_tests = n_aux_tests + 1
229 IF (orcond .LT. eps)
THEN
239 IF (info .GT. 0 .AND. info .LE. n+1)
THEN
241 WRITE (*, fmt=8000) c2, n, info, orcond, rcond
248 diff(i,j) = x(i,j) - invhilb(i,j)
255 IF (
lsamen( 2, c2,
'PO' ) .OR.
lsamen( 2, c2,
'SY' ) )
THEN
260 sumr = sumr + s(i) * abs(a(i,j)) * s(j)
261 sumri = sumri + abs(invhilb(i, j)) / (s(j) * s(i))
264 rnorm = max(rnorm,sumr)
265 rinorm = max(rinorm,sumri)
267 ELSE IF (
lsamen( 2, c2,
'GE' ) .OR.
lsamen( 2, c2,
'GB' ) )
273 sumr = sumr + r(i) * abs(a(i,j)) * c(j)
274 sumri = sumri + abs(invhilb(i, j)) / (r(j) * c(i))
276 rnorm = max(rnorm,sumr)
277 rinorm = max(rinorm,sumri)
281 rnorm = rnorm / abs(a(1, 1))
282 rcond = 1.0d+0/(rnorm * rinorm)
290 rinv(i) = rinv(i) + abs(a(i,j))
299 sumri = sumri + abs(invhilb(i,j) * rinv(j))
301 rinorm = max(rinorm, sumri)
306 ncond = abs(a(1,1)) / rinorm
316 normt = max(abs(invhilb(i, k)), normt)
317 normdif = max(abs(x(i,k) - invhilb(i,k)), normdif)
318 IF (invhilb(i,k) .NE. 0.0d+0)
THEN
319 cwise_err = max(abs(x(i,k) - invhilb(i,k))
320 $ /abs(invhilb(i,k)), cwise_err)
321 ELSE IF (x(i, k) .NE. 0.0d+0)
THEN
322 cwise_err =
dlamch(
'OVERFLOW')
325 IF (normt .NE. 0.0d+0)
THEN
326 nwise_err = normdif / normt
327 ELSE IF (normdif .NE. 0.0d+0)
THEN
328 nwise_err =
dlamch(
'OVERFLOW')
338 rinv(i) = rinv(i) + abs(a(i, j) * invhilb(j, k))
346 $ + abs(invhilb(i, j) * rinv(j) / invhilb(i, k))
348 rinorm = max(rinorm, sumri)
352 ccond = abs(a(1,1))/rinorm
355 nwise_bnd = errbnd_n(k + (bnd_i-1)*nrhs)
356 cwise_bnd = errbnd_c(k + (bnd_i-1)*nrhs)
357 nwise_rcond = errbnd_n(k + (cond_i-1)*nrhs)
358 cwise_rcond = errbnd_c(k + (cond_i-1)*nrhs)
362 IF (ncond .GE. condthresh)
THEN
364 IF (nwise_bnd .GT. errthresh)
THEN
365 tstrat(1) = 1/(2.0d+0*eps)
367 IF (nwise_bnd .NE. 0.0d+0)
THEN
368 tstrat(1) = nwise_err / nwise_bnd
369 ELSE IF (nwise_err .NE. 0.0d+0)
THEN
370 tstrat(1) = 1/(16.0*eps)
374 IF (tstrat(1) .GT. 1.0d+0)
THEN
375 tstrat(1) = 1/(4.0d+0*eps)
380 IF (nwise_bnd .LT. 1.0d+0)
THEN
381 tstrat(1) = 1/(8.0d+0*eps)
389 IF (ccond .GE. condthresh)
THEN
391 IF (cwise_bnd .GT. errthresh)
THEN
392 tstrat(2) = 1/(2.0d+0*eps)
394 IF (cwise_bnd .NE. 0.0d+0)
THEN
395 tstrat(2) = cwise_err / cwise_bnd
396 ELSE IF (cwise_err .NE. 0.0d+0)
THEN
397 tstrat(2) = 1/(16.0d+0*eps)
401 IF (tstrat(2) .GT. 1.0d+0) tstrat(2) = 1/(4.0d+0*eps)
405 IF (cwise_bnd .LT. 1.0d+0)
THEN
406 tstrat(2) = 1/(8.0d+0*eps)
413 tstrat(3) = berr(k)/eps
416 tstrat(4) = rcond / orcond
417 IF (rcond .GE. condthresh .AND. tstrat(4) .LT. 1.0d+0)
418 $ tstrat(4) = 1.0d+0 / tstrat(4)
420 tstrat(5) = ncond / nwise_rcond
421 IF (ncond .GE. condthresh .AND. tstrat(5) .LT. 1.0d+0)
422 $ tstrat(5) = 1.0d+0 / tstrat(5)
424 tstrat(6) = ccond / nwise_rcond
425 IF (ccond .GE. condthresh .AND. tstrat(6) .LT. 1.0d+0)
426 $ tstrat(6) = 1.0d+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(
' D', a2,
'SVXX: N =', i2,
', RHS = ', i2,
472 $
', NWISE GUAR. = ', a,
', CWISE GUAR. = ', a,
473 $
' test(',i1,
') =', g12.5 )
474 9998
FORMAT(
' D', a2,
'SVXX: ', i6,
' out of ', i6,
475 $
' tests failed to pass the threshold' )
476 9997
FORMAT(
' D', 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(
' D', a2,
'SVXX: N =', i2,
', INFO = ', i3,
492 $
', ORCOND = ', g12.5,
', real RCOND = ', g12.5 )
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(CMACH)
DLAMCH
subroutine dlahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
DLAHILB
subroutine dgesvxx(FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO)
DGESVXX computes the solution to system of linear equations A * X = B for GE matrices ...
logical function lsamen(N, CA, CB)
LSAMEN
subroutine dgbsvxx(FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO)
DGBSVXX computes the solution to system of linear equations A * X = B for GB matrices ...
subroutine dsysvxx(FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, EQUED, S, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO)
DSYSVXX
subroutine dposvxx(FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, S, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO)
DPOSVXX computes the solution to system of linear equations A * X = B for PO matrices ...