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
123 DOUBLE PRECISION TSTRAT(ntests), RINV(nmax), PARAMS(nparams),
124 $ s(nmax),r(nmax),c(nmax),rwork(3*nmax),
126 $ errbnd_n(nmax*3), errbnd_c(nmax*3)
128 COMPLEX*16 A(nmax,nmax),INVHILB(nmax,nmax),X(nmax,nmax),
129 $ work(nmax*3*5), af(nmax, nmax),b(nmax, nmax),
131 $ ab( (nmax-1)+(nmax-1)+1, nmax ),
132 $ abcopy( (nmax-1)+(nmax-1)+1, nmax ),
133 $ afb( 2*(nmax-1)+(nmax-1)+1, nmax )
136 DOUBLE PRECISION DLAMCH
144 INTRINSIC sqrt, max, abs, dble, dimag
147 DOUBLE PRECISION CABS1
150 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
153 INTEGER NWISE_I, CWISE_I
154 parameter (nwise_i = 1, cwise_i = 1)
155 INTEGER BND_I, COND_I
156 parameter (bnd_i = 2, cond_i = 3)
164 eps = dlamch(
'Epsilon')
168 ldab = (nmax-1)+(nmax-1)+1
169 ldafb = 2*(nmax-1)+(nmax-1)+1
174 printed_guide = .false.
183 m = max(sqrt(dble(n)), 10.0d+0)
187 CALL zlahilb(n, n, a, lda, invhilb, lda, b,
188 $ lda, work, info, path)
191 CALL zlacpy(
'ALL', n, n, a, nmax, acopy, nmax)
196 ab( i, j ) = (0.0d+0,0.0d+0)
200 DO i = max( 1, j-ku ), min( n, j+kl )
201 ab( ku+1+i-j, j ) = a( i, j )
208 abcopy( i, j ) = (0.0d+0,0.0d+0)
211 CALL zlacpy(
'ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)
214 IF ( lsamen( 2, c2,
'SY' ) )
THEN
215 CALL zsysvxx(fact, uplo, n, nrhs, acopy, lda, af, lda,
216 $ ipiv, equed, s, b, lda, x, lda, orcond,
217 $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
218 $ params, work, rwork, info)
219 ELSE IF ( lsamen( 2, c2,
'PO' ) )
THEN
220 CALL zposvxx(fact, uplo, n, nrhs, acopy, lda, af, lda,
221 $ equed, s, b, lda, x, lda, orcond,
222 $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
223 $ params, work, rwork, info)
224 ELSE IF ( lsamen( 2, c2,
'HE' ) )
THEN
225 CALL zhesvxx(fact, uplo, n, nrhs, acopy, lda, af, lda,
226 $ ipiv, equed, s, b, lda, x, lda, orcond,
227 $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
228 $ params, work, rwork, info)
229 ELSE IF ( lsamen( 2, c2,
'GB' ) )
THEN
230 CALL zgbsvxx(fact, trans, n, kl, ku, nrhs, abcopy,
231 $ ldab, afb, ldafb, ipiv, equed, r, c, b,
232 $ lda, x, lda, orcond, rpvgrw, berr, nerrbnd,
233 $ errbnd_n, errbnd_c, nparams, params, work, rwork,
236 CALL zgesvxx(fact, trans, n, nrhs, acopy, lda, af, lda,
237 $ ipiv, equed, r, c, b, lda, x, lda, orcond,
238 $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
239 $ params, work, rwork, info)
242 n_aux_tests = n_aux_tests + 1
243 IF (orcond .LT. eps)
THEN
253 IF (info .GT. 0 .AND. info .LE. n+1)
THEN
255 WRITE (*, fmt=8000) c2, n, info, orcond, rcond
262 diff(i,j) = x(i,j) - invhilb(i,j)
269 IF ( lsamen( 2, c2,
'PO' ) .OR. lsamen( 2, c2,
'SY' ) .OR.
270 $ lsamen( 2, c2,
'HE' ) )
THEN
275 sumr = sumr + s(i) * cabs1(a(i,j)) * s(j)
276 sumri = sumri + cabs1(invhilb(i, j)) / (s(j) * s(i))
278 rnorm = max(rnorm,sumr)
279 rinorm = max(rinorm,sumri)
281 ELSE IF ( lsamen( 2, c2,
'GE' ) .OR. lsamen( 2, c2,
'GB' ) )
287 sumr = sumr + r(i) * cabs1(a(i,j)) * c(j)
288 sumri = sumri + cabs1(invhilb(i, j)) / (r(j) * c(i))
290 rnorm = max(rnorm,sumr)
291 rinorm = max(rinorm,sumri)
295 rnorm = rnorm / cabs1(a(1, 1))
296 rcond = 1.0d+0/(rnorm * rinorm)
304 rinv(i) = rinv(i) + cabs1(a(i,j))
313 sumri = sumri + cabs1(invhilb(i,j) * rinv(j))
315 rinorm = max(rinorm, sumri)
320 ncond = cabs1(a(1,1)) / rinorm
330 normt = max(cabs1(invhilb(i, k)), normt)
331 normdif = max(cabs1(x(i,k) - invhilb(i,k)), normdif)
332 IF (invhilb(i,k) .NE. 0.0d+0)
THEN
333 cwise_err = max(cabs1(x(i,k) - invhilb(i,k))
334 $ /cabs1(invhilb(i,k)), cwise_err)
335 ELSE IF (x(i, k) .NE. 0.0d+0)
THEN
336 cwise_err = dlamch(
'OVERFLOW')
339 IF (normt .NE. 0.0d+0)
THEN
340 nwise_err = normdif / normt
341 ELSE IF (normdif .NE. 0.0d+0)
THEN
342 nwise_err = dlamch(
'OVERFLOW')
352 rinv(i) = rinv(i) + cabs1(a(i, j) * invhilb(j, k))
360 $ + cabs1(invhilb(i, j) * rinv(j) / invhilb(i, k))
362 rinorm = max(rinorm, sumri)
366 ccond = cabs1(a(1,1))/rinorm
369 nwise_bnd = errbnd_n(k + (bnd_i-1)*nrhs)
370 cwise_bnd = errbnd_c(k + (bnd_i-1)*nrhs)
371 nwise_rcond = errbnd_n(k + (cond_i-1)*nrhs)
372 cwise_rcond = errbnd_c(k + (cond_i-1)*nrhs)
376 IF (ncond .GE. condthresh)
THEN
378 IF (nwise_bnd .GT. errthresh)
THEN
379 tstrat(1) = 1/(2.0d+0*eps)
381 IF (nwise_bnd .NE. 0.0d+0)
THEN
382 tstrat(1) = nwise_err / nwise_bnd
383 ELSE IF (nwise_err .NE. 0.0d+0)
THEN
384 tstrat(1) = 1/(16.0*eps)
388 IF (tstrat(1) .GT. 1.0d+0)
THEN
389 tstrat(1) = 1/(4.0d+0*eps)
394 IF (nwise_bnd .LT. 1.0d+0)
THEN
395 tstrat(1) = 1/(8.0d+0*eps)
403 IF (ccond .GE. condthresh)
THEN
405 IF (cwise_bnd .GT. errthresh)
THEN
406 tstrat(2) = 1/(2.0d+0*eps)
408 IF (cwise_bnd .NE. 0.0d+0)
THEN
409 tstrat(2) = cwise_err / cwise_bnd
410 ELSE IF (cwise_err .NE. 0.0d+0)
THEN
411 tstrat(2) = 1/(16.0d+0*eps)
415 IF (tstrat(2) .GT. 1.0d+0) tstrat(2) = 1/(4.0d+0*eps)
419 IF (cwise_bnd .LT. 1.0d+0)
THEN
420 tstrat(2) = 1/(8.0d+0*eps)
427 tstrat(3) = berr(k)/eps
430 tstrat(4) = rcond / orcond
431 IF (rcond .GE. condthresh .AND. tstrat(4) .LT. 1.0d+0)
432 $ tstrat(4) = 1.0d+0 / tstrat(4)
434 tstrat(5) = ncond / nwise_rcond
435 IF (ncond .GE. condthresh .AND. tstrat(5) .LT. 1.0d+0)
436 $ tstrat(5) = 1.0d+0 / tstrat(5)
438 tstrat(6) = ccond / nwise_rcond
439 IF (ccond .GE. condthresh .AND. tstrat(6) .LT. 1.0d+0)
440 $ tstrat(6) = 1.0d+0 / tstrat(6)
443 IF (tstrat(i) .GT. thresh)
THEN
444 IF (.NOT.printed_guide)
THEN
455 printed_guide = .true.
457 WRITE( *, 9999) c2, n, k, nguar, cguar, i, tstrat(i)
480 IF( nfail .GT. 0 )
THEN
481 WRITE(*,9998) c2, nfail, ntests*n+n_aux_tests
485 9999
FORMAT(
' Z', a2,
'SVXX: N =', i2,
', RHS = ', i2,
486 $
', NWISE GUAR. = ', a,
', CWISE GUAR. = ', a,
487 $
' test(',i1,
') =', g12.5 )
488 9998
FORMAT(
' Z', a2,
'SVXX: ', i6,
' out of ', i6,
489 $
' tests failed to pass the threshold' )
490 9997
FORMAT(
' Z', a2,
'SVXX passed the tests of error bounds' )
492 9996
FORMAT( 3x, i2,
': Normwise guaranteed forward error', / 5x,
493 $
'Guaranteed case: if norm ( abs( Xc - Xt )',
494 $ .LE.
' / norm ( Xt ) ERRBND( *, nwise_i, bnd_i ), then',
496 $ .LE.
'ERRBND( *, nwise_i, bnd_i ) MAX(SQRT(N), 10) * EPS')
497 9995
FORMAT( 3x, i2,
': Componentwise guaranteed forward error' )
498 9994
FORMAT( 3x, i2,
': Backwards error' )
499 9993
FORMAT( 3x, i2,
': Reciprocal condition number' )
500 9992
FORMAT( 3x, i2,
': Reciprocal normwise condition number' )
501 9991
FORMAT( 3x, i2,
': Raw normwise error estimate' )
502 9990
FORMAT( 3x, i2,
': Reciprocal componentwise condition number' )
503 9989
FORMAT( 3x, i2,
': Raw componentwise error estimate' )
505 8000
FORMAT(
' Z', a2,
'SVXX: N =', i2,
', INFO = ', i3,
506 $
', ORCOND = ', g12.5,
', real RCOND = ', g12.5 )
subroutine zgesvxx(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, RWORK, INFO)
ZGESVXX computes the solution to system of linear equations A * X = B for GE matrices ...
subroutine zgbsvxx(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, RWORK, INFO)
ZGBSVXX computes the solution to system of linear equations A * X = B for GB matrices ...
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zsysvxx(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, RWORK, INFO)
ZSYSVXX computes the solution to system of linear equations A * X = B for SY matrices ...
logical function lsamen(N, CA, CB)
LSAMEN
subroutine zlahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO, PATH)
ZLAHILB
subroutine zebchvxx(THRESH, PATH)
ZEBCHVXX
subroutine zhesvxx(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, RWORK, INFO)
ZHESVXX computes the solution to system of linear equations A * X = B for HE matrices ...
subroutine zposvxx(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, RWORK, INFO)
ZPOSVXX computes the solution to system of linear equations A * X = B for PO matrices ...