103 INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU
104 parameter (nmax = 6, 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 REAL 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 REAL 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 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 )
144 INTRINSIC sqrt, max, abs,
REAL, AIMAG
150 cabs1( zdum ) = abs(
REAL( ZDUM ) ) + abs( AIMAG( 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 = slamch(
'Epsilon')
168 ldab = (nmax-1)+(nmax-1)+1
169 ldafb = 2*(nmax-1)+(nmax-1)+1
174 printed_guide = .false.
183 m = max(sqrt(
REAL(n)), 10.0)
187 CALL clahilb(n, n, a, lda, invhilb, lda, b,
188 $ lda, work, info, path)
191 CALL clacpy(
'ALL', n, n, a, nmax, acopy, nmax)
196 ab( i, j ) = (0.0e+0,0.0e+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.0e+0,0.0e+0)
211 CALL clacpy(
'ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)
214 IF ( lsamen( 2, c2,
'SY' ) )
THEN
215 CALL csysvxx(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 cposvxx(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 chesvxx(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 cgbsvxx(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 cgesvxx(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.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.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.0)
THEN
336 cwise_err = slamch(
'OVERFLOW')
339 IF (normt .NE. 0.0)
THEN
340 nwise_err = normdif / normt
341 ELSE IF (normdif .NE. 0.0)
THEN
342 nwise_err = slamch(
'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.0*eps)
381 IF (nwise_bnd .NE. 0.0)
THEN
382 tstrat(1) = nwise_err / nwise_bnd
383 ELSE IF (nwise_err .NE. 0.0)
THEN
384 tstrat(1) = 1/(16.0*eps)
388 IF (tstrat(1) .GT. 1.0)
THEN
389 tstrat(1) = 1/(4.0*eps)
394 IF (nwise_bnd .LT. 1.0)
THEN
395 tstrat(1) = 1/(8.0*eps)
403 IF (ccond .GE. condthresh)
THEN
405 IF (cwise_bnd .GT. errthresh)
THEN
406 tstrat(2) = 1/(2.0*eps)
408 IF (cwise_bnd .NE. 0.0)
THEN
409 tstrat(2) = cwise_err / cwise_bnd
410 ELSE IF (cwise_err .NE. 0.0)
THEN
411 tstrat(2) = 1/(16.0*eps)
415 IF (tstrat(2) .GT. 1.0) tstrat(2) = 1/(4.0*eps)
419 IF (cwise_bnd .LT. 1.0)
THEN
420 tstrat(2) = 1/(8.0*eps)
427 tstrat(3) = berr(k)/eps
430 tstrat(4) = rcond / orcond
431 IF (rcond .GE. condthresh .AND. tstrat(4) .LT. 1.0)
432 $ tstrat(4) = 1.0 / tstrat(4)
434 tstrat(5) = ncond / nwise_rcond
435 IF (ncond .GE. condthresh .AND. tstrat(5) .LT. 1.0)
436 $ tstrat(5) = 1.0 / tstrat(5)
438 tstrat(6) = ccond / nwise_rcond
439 IF (ccond .GE. condthresh .AND. tstrat(6) .LT. 1.0)
440 $ tstrat(6) = 1.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(
' C', a2,
'SVXX: N =', i2,
', RHS = ', i2,
486 $
', NWISE GUAR. = ', a,
', CWISE GUAR. = ', a,
487 $
' test(',i1,
') =', g12.5 )
488 9998
FORMAT(
' C', a2,
'SVXX: ', i6,
' out of ', i6,
489 $
' tests failed to pass the threshold' )
490 9997
FORMAT(
' C', 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(
' C', a2,
'SVXX: N =', i2,
', INFO = ', i3,
506 $
', ORCOND = ', g12.5,
', real RCOND = ', g12.5 )
subroutine cposvxx(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)
CPOSVXX computes the solution to system of linear equations A * X = B for PO matrices ...
subroutine chesvxx(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)
CHESVXX computes the solution to system of linear equations A * X = B for HE matrices ...
subroutine cebchvxx(THRESH, PATH)
CEBCHVXX
subroutine csysvxx(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)
CSYSVXX computes the solution to system of linear equations A * X = B for SY matrices ...
logical function lsamen(N, CA, CB)
LSAMEN
subroutine cgbsvxx(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)
CGBSVXX computes the solution to system of linear equations A * X = B for GB matrices ...
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine clahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO, PATH)
CLAHILB
subroutine cgesvxx(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)
CGESVXX computes the solution to system of linear equations A * X = B for GE matrices ...