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 )