LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
sebchvxx.f
Go to the documentation of this file.
1 *> \brief \b SEBCHVXX
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE SEBCHVXX( THRESH, PATH )
12 *
13 * .. Scalar Arguments ..
14 * REAL THRESH
15 * CHARACTER*3 PATH
16 * ..
17 *
18 *
19 *> \par Purpose:
20 * =============
21 *>
22 *> \verbatim
23 *>
24 *> SEBCHVXX will run S**SVXX on a series of Hilbert matrices and then
25 *> compare the error bounds returned by SGESVXX to see if the returned
26 *> answer indeed falls within those bounds.
27 *>
28 *> Eight test ratios will be computed. The tests will pass if they are .LT.
29 *> THRESH. There are two cases that are determined by 1 / (SQRT( N ) * EPS).
30 *> If that value is .LE. to the component wise reciprocal condition number,
31 *> it uses the guaranteed case, other wise it uses the unguaranteed case.
32 *>
33 *> Test ratios:
34 *> Let Xc be X_computed and Xt be X_truth.
35 *> The norm used is the infinity norm.
36 *>
37 *> Let A be the guaranteed case and B be the unguaranteed case.
38 *>
39 *> 1. Normwise guaranteed forward error bound.
40 *> A: norm ( abs( Xc - Xt ) / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ) and
41 *> ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N),10) * EPS.
42 *> If these conditions are met, the test ratio is set to be
43 *> ERRBND( *, nwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS.
44 *> B: For this case, SGESVXX should just return 1. If it is less than
45 *> one, treat it the same as in 1A. Otherwise it fails. (Set test
46 *> ratio to ERRBND( *, nwise_i, bnd_i ) * THRESH?)
47 *>
48 *> 2. Componentwise guaranteed forward error bound.
49 *> A: norm ( abs( Xc(j) - Xt(j) ) ) / norm (Xt(j)) .LE. ERRBND( *, cwise_i, bnd_i )
50 *> for all j .AND. ERRBND( *, cwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS.
51 *> If these conditions are met, the test ratio is set to be
52 *> ERRBND( *, cwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS.
53 *> B: Same as normwise test ratio.
54 *>
55 *> 3. Backwards error.
56 *> A: The test ratio is set to BERR/EPS.
57 *> B: Same test ratio.
58 *>
59 *> 4. Reciprocal condition number.
60 *> A: A condition number is computed with Xt and compared with the one
61 *> returned from SGESVXX. Let RCONDc be the RCOND returned by SGESVXX
62 *> and RCONDt be the RCOND from the truth value. Test ratio is set to
63 *> MAX(RCONDc/RCONDt, RCONDt/RCONDc).
64 *> B: Test ratio is set to 1 / (EPS * RCONDc).
65 *>
66 *> 5. Reciprocal normwise condition number.
67 *> A: The test ratio is set to
68 *> MAX(ERRBND( *, nwise_i, cond_i ) / NCOND, NCOND / ERRBND( *, nwise_i, cond_i )).
69 *> B: Test ratio is set to 1 / (EPS * ERRBND( *, nwise_i, cond_i )).
70 *>
71 *> 7. Reciprocal componentwise condition number.
72 *> A: Test ratio is set to
73 *> MAX(ERRBND( *, cwise_i, cond_i ) / CCOND, CCOND / ERRBND( *, cwise_i, cond_i )).
74 *> B: Test ratio is set to 1 / (EPS * ERRBND( *, cwise_i, cond_i )).
75 *>
76 *> .. Parameters ..
77 *> NMAX is determined by the largest number in the inverse of the Hilbert
78 *> matrix. Precision is exhausted when the largest entry in it is greater
79 *> than 2 to the power of the number of bits in the fraction of the data
80 *> type used plus one, which is 24 for single precision.
81 *> NMAX should be 6 for single and 11 for double.
82 *> \endverbatim
83 *
84 * Authors:
85 * ========
86 *
87 *> \author Univ. of Tennessee
88 *> \author Univ. of California Berkeley
89 *> \author Univ. of Colorado Denver
90 *> \author NAG Ltd.
91 *
92 *> \date November 2011
93 *
94 *> \ingroup single_lin
95 *
96 * =====================================================================
97  SUBROUTINE sebchvxx( THRESH, PATH )
98  IMPLICIT NONE
99 * .. Scalar Arguments ..
100  REAL thresh
101  CHARACTER*3 path
102 
103  INTEGER nmax, nparams, nerrbnd, ntests, kl, ku
104  parameter(nmax = 6, nparams = 2, nerrbnd = 3,
105  $ ntests = 6)
106 
107 * .. Local Scalars ..
108  INTEGER n, nrhs, info, i ,j, k, nfail, lda, ldab,
109  $ ldafb, n_aux_tests
110  CHARACTER fact, trans, uplo, equed
111  CHARACTER*2 c2
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
120 
121 * .. Local Arrays ..
122  REAL tstrat(ntests), rinv(nmax), params(nparams),
123  $ a(nmax, nmax), acopy(nmax, nmax),
124  $ invhilb(nmax, nmax), r(nmax), c(nmax), s(nmax),
125  $ work(nmax * 5), b(nmax, nmax), x(nmax, nmax),
126  $ diff(nmax, nmax), af(nmax, nmax),
127  $ ab( (nmax-1)+(nmax-1)+1, nmax ),
128  $ abcopy( (nmax-1)+(nmax-1)+1, nmax ),
129  $ afb( 2*(nmax-1)+(nmax-1)+1, nmax ),
130  $ errbnd_n(nmax*3), errbnd_c(nmax*3)
131  INTEGER iwork(nmax), ipiv(nmax)
132 
133 * .. External Functions ..
134  REAL slamch
135 
136 * .. External Subroutines ..
137  EXTERNAL slahilb, sgesvxx, ssysvxx, sposvxx, sgbsvxx,
138  $ slacpy, lsamen
139  LOGICAL lsamen
140 
141 * .. Intrinsic Functions ..
142  INTRINSIC sqrt, max, abs
143 
144 * .. Parameters ..
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)
149 
150 * Create the loop to test out the Hilbert matrices
151 
152  fact = 'E'
153  uplo = 'U'
154  trans = 'N'
155  equed = 'N'
156  eps = slamch('Epsilon')
157  nfail = 0
158  n_aux_tests = 0
159  lda = nmax
160  ldab = (nmax-1)+(nmax-1)+1
161  ldafb = 2*(nmax-1)+(nmax-1)+1
162  c2 = path( 2: 3 )
163 
164 * Main loop to test the different Hilbert Matrices.
165 
166  printed_guide = .false.
167 
168  DO n = 1 , nmax
169  params(1) = -1
170  params(2) = -1
171 
172  kl = n-1
173  ku = n-1
174  nrhs = n
175  m = max(sqrt(REAL(n)), 10.0)
176 
177 * Generate the Hilbert matrix, its inverse, and the
178 * right hand side, all scaled by the LCM(1,..,2N-1).
179  CALL slahilb(n, n, a, lda, invhilb, lda, b, lda, work, info)
180 
181 * Copy A into ACOPY.
182  CALL slacpy('ALL', n, n, a, nmax, acopy, nmax)
183 
184 * Store A in band format for GB tests
185  DO j = 1, n
186  DO i = 1, kl+ku+1
187  ab( i, j ) = 0.0e+0
188  END DO
189  END DO
190  DO j = 1, n
191  DO i = max( 1, j-ku ), min( n, j+kl )
192  ab( ku+1+i-j, j ) = a( i, j )
193  END DO
194  END DO
195 
196 * Copy AB into ABCOPY.
197  DO j = 1, n
198  DO i = 1, kl+ku+1
199  abcopy( i, j ) = 0.0e+0
200  END DO
201  END DO
202  CALL slacpy('ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)
203 
204 * Call S**SVXX with default PARAMS and N_ERR_BND = 3.
205  IF ( lsamen( 2, c2, 'SY' ) ) THEN
206  CALL ssysvxx(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 sposvxx(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 sgbsvxx(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,
220  $ info)
221  ELSE
222  CALL sgesvxx(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)
226  END IF
227 
228  n_aux_tests = n_aux_tests + 1
229  IF (orcond .LT. eps) THEN
230 ! Either factorization failed or the matrix is flagged, and 1 <=
231 ! INFO <= N+1. We don't decide based on rcond anymore.
232 ! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN
233 ! NFAIL = NFAIL + 1
234 ! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND
235 ! END IF
236  ELSE
237 ! Either everything succeeded (INFO == 0) or some solution failed
238 ! to converge (INFO > N+1).
239  IF (info .GT. 0 .AND. info .LE. n+1) THEN
240  nfail = nfail + 1
241  WRITE (*, fmt=8000) c2, n, info, orcond, rcond
242  END IF
243  END IF
244 
245 * Calculating the difference between S**SVXX's X and the true X.
246  DO i = 1, n
247  DO j = 1, nrhs
248  diff( i, j ) = x( i, j ) - invhilb( i, j )
249  END DO
250  END DO
251 
252 * Calculating the RCOND
253  rnorm = 0
254  rinorm = 0
255  IF ( lsamen( 2, c2, 'PO' ) .OR. lsamen( 2, c2, 'SY' ) ) THEN
256  DO i = 1, n
257  sumr = 0
258  sumri = 0
259  DO j = 1, n
260  sumr = sumr + abs(s(i) * a(i,j) * s(j))
261  sumri = sumri + abs(invhilb(i, j) / s(j) / s(i))
262  END DO
263  rnorm = max(rnorm,sumr)
264  rinorm = max(rinorm,sumri)
265  END DO
266  ELSE IF ( lsamen( 2, c2, 'GE' ) .OR. lsamen( 2, c2, 'GB' ) )
267  $ THEN
268  DO i = 1, n
269  sumr = 0
270  sumri = 0
271  DO j = 1, n
272  sumr = sumr + abs(r(i) * a(i,j) * c(j))
273  sumri = sumri + abs(invhilb(i, j) / r(j) / c(i))
274  END DO
275  rnorm = max(rnorm,sumr)
276  rinorm = max(rinorm,sumri)
277  END DO
278  END IF
279 
280  rnorm = rnorm / a(1, 1)
281  rcond = 1.0/(rnorm * rinorm)
282 
283 * Calculating the R for normwise rcond.
284  DO i = 1, n
285  rinv(i) = 0.0
286  END DO
287  DO j = 1, n
288  DO i = 1, n
289  rinv(i) = rinv(i) + abs(a(i,j))
290  END DO
291  END DO
292 
293 * Calculating the Normwise rcond.
294  rinorm = 0.0
295  DO i = 1, n
296  sumri = 0.0
297  DO j = 1, n
298  sumri = sumri + abs(invhilb(i,j) * rinv(j))
299  END DO
300  rinorm = max(rinorm, sumri)
301  END DO
302 
303 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
304 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
305  ncond = a(1,1) / rinorm
306 
307  condthresh = m * eps
308  errthresh = m * eps
309 
310  DO k = 1, nrhs
311  normt = 0.0
312  normdif = 0.0
313  cwise_err = 0.0
314  DO i = 1, n
315  normt = max(abs(invhilb(i, k)), normt)
316  normdif = max(abs(x(i,k) - invhilb(i,k)), normdif)
317  IF (invhilb(i,k) .NE. 0.0) THEN
318  cwise_err = max(abs(x(i,k) - invhilb(i,k))
319  $ /abs(invhilb(i,k)), cwise_err)
320  ELSE IF (x(i, k) .NE. 0.0) THEN
321  cwise_err = slamch('OVERFLOW')
322  END IF
323  END DO
324  IF (normt .NE. 0.0) THEN
325  nwise_err = normdif / normt
326  ELSE IF (normdif .NE. 0.0) THEN
327  nwise_err = slamch('OVERFLOW')
328  ELSE
329  nwise_err = 0.0
330  ENDIF
331 
332  DO i = 1, n
333  rinv(i) = 0.0
334  END DO
335  DO j = 1, n
336  DO i = 1, n
337  rinv(i) = rinv(i) + abs(a(i, j) * invhilb(j, k))
338  END DO
339  END DO
340  rinorm = 0.0
341  DO i = 1, n
342  sumri = 0.0
343  DO j = 1, n
344  sumri = sumri
345  $ + abs(invhilb(i, j) * rinv(j) / invhilb(i, k))
346  END DO
347  rinorm = max(rinorm, sumri)
348  END DO
349 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
350 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
351  ccond = a(1,1)/rinorm
352 
353 ! Forward error bound tests
354  nwise_bnd = errbnd_n(k + (bnd_i-1)*nrhs)
355  cwise_bnd = errbnd_c(k + (bnd_i-1)*nrhs)
356  nwise_rcond = errbnd_n(k + (cond_i-1)*nrhs)
357  cwise_rcond = errbnd_c(k + (cond_i-1)*nrhs)
358 ! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
359 ! $ condthresh, ncond.ge.condthresh
360 ! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
361 
362  IF (ncond .GE. condthresh) THEN
363  nguar = 'YES'
364  IF (nwise_bnd .GT. errthresh) THEN
365  tstrat(1) = 1/(2.0*eps)
366  ELSE
367 
368  IF (nwise_bnd .NE. 0.0) THEN
369  tstrat(1) = nwise_err / nwise_bnd
370  ELSE IF (nwise_err .NE. 0.0) THEN
371  tstrat(1) = 1/(16.0*eps)
372  ELSE
373  tstrat(1) = 0.0
374  END IF
375  IF (tstrat(1) .GT. 1.0) THEN
376  tstrat(1) = 1/(4.0*eps)
377  END IF
378  END IF
379  ELSE
380  nguar = 'NO'
381  IF (nwise_bnd .LT. 1.0) THEN
382  tstrat(1) = 1/(8.0*eps)
383  ELSE
384  tstrat(1) = 1.0
385  END IF
386  END IF
387 ! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
388 ! $ condthresh, ccond.ge.condthresh
389 ! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
390  IF (ccond .GE. condthresh) THEN
391  cguar = 'YES'
392 
393  IF (cwise_bnd .GT. errthresh) THEN
394  tstrat(2) = 1/(2.0*eps)
395  ELSE
396  IF (cwise_bnd .NE. 0.0) THEN
397  tstrat(2) = cwise_err / cwise_bnd
398  ELSE IF (cwise_err .NE. 0.0) THEN
399  tstrat(2) = 1/(16.0*eps)
400  ELSE
401  tstrat(2) = 0.0
402  END IF
403  IF (tstrat(2) .GT. 1.0) tstrat(2) = 1/(4.0*eps)
404  END IF
405  ELSE
406  cguar = 'NO'
407  IF (cwise_bnd .LT. 1.0) THEN
408  tstrat(2) = 1/(8.0*eps)
409  ELSE
410  tstrat(2) = 1.0
411  END IF
412  END IF
413 
414 ! Backwards error test
415  tstrat(3) = berr(k)/eps
416 
417 ! Condition number tests
418  tstrat(4) = rcond / orcond
419  IF (rcond .GE. condthresh .AND. tstrat(4) .LT. 1.0)
420  $ tstrat(4) = 1.0 / tstrat(4)
421 
422  tstrat(5) = ncond / nwise_rcond
423  IF (ncond .GE. condthresh .AND. tstrat(5) .LT. 1.0)
424  $ tstrat(5) = 1.0 / tstrat(5)
425 
426  tstrat(6) = ccond / nwise_rcond
427  IF (ccond .GE. condthresh .AND. tstrat(6) .LT. 1.0)
428  $ tstrat(6) = 1.0 / tstrat(6)
429 
430  DO i = 1, ntests
431  IF (tstrat(i) .GT. thresh) THEN
432  IF (.NOT.printed_guide) THEN
433  WRITE(*,*)
434  WRITE( *, 9996) 1
435  WRITE( *, 9995) 2
436  WRITE( *, 9994) 3
437  WRITE( *, 9993) 4
438  WRITE( *, 9992) 5
439  WRITE( *, 9991) 6
440  WRITE( *, 9990) 7
441  WRITE( *, 9989) 8
442  WRITE(*,*)
443  printed_guide = .true.
444  END IF
445  WRITE( *, 9999) c2, n, k, nguar, cguar, i, tstrat(i)
446  nfail = nfail + 1
447  END IF
448  END DO
449  END DO
450 
451 c$$$ WRITE(*,*)
452 c$$$ WRITE(*,*) 'Normwise Error Bounds'
453 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i)
454 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i)
455 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i)
456 c$$$ WRITE(*,*)
457 c$$$ WRITE(*,*) 'Componentwise Error Bounds'
458 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i)
459 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i)
460 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i)
461 c$$$ print *, 'Info: ', info
462 c$$$ WRITE(*,*)
463 * WRITE(*,*) 'TSTRAT: ',TSTRAT
464 
465  END DO
466 
467  WRITE(*,*)
468  IF( nfail .GT. 0 ) THEN
469  WRITE(*,9998) c2, nfail, ntests*n+n_aux_tests
470  ELSE
471  WRITE(*,9997) c2
472  END IF
473  9999 format( ' S', a2, 'SVXX: N =', i2, ', RHS = ', i2,
474  $ ', NWISE GUAR. = ', a, ', CWISE GUAR. = ', a,
475  $ ' test(',i1,') =', g12.5 )
476  9998 format( ' S', a2, 'SVXX: ', i6, ' out of ', i6,
477  $ ' tests failed to pass the threshold' )
478  9997 format( ' S', a2, 'SVXX passed the tests of error bounds' )
479 * Test ratios.
480  9996 format( 3x, i2, ': Normwise guaranteed forward error', / 5x,
481  $ 'Guaranteed case: if norm ( abs( Xc - Xt )',
482  $ .LE.' / norm ( Xt ) ERRBND( *, nwise_i, bnd_i ), then',
483  $ / 5x,
484  $ .LE.'ERRBND( *, nwise_i, bnd_i ) MAX(SQRT(N), 10) * EPS')
485  9995 format( 3x, i2, ': Componentwise guaranteed forward error' )
486  9994 format( 3x, i2, ': Backwards error' )
487  9993 format( 3x, i2, ': Reciprocal condition number' )
488  9992 format( 3x, i2, ': Reciprocal normwise condition number' )
489  9991 format( 3x, i2, ': Raw normwise error estimate' )
490  9990 format( 3x, i2, ': Reciprocal componentwise condition number' )
491  9989 format( 3x, i2, ': Raw componentwise error estimate' )
492 
493  8000 format( ' S', a2, 'SVXX: N =', i2, ', INFO = ', i3,
494  $ ', ORCOND = ', g12.5, ', real RCOND = ', g12.5 )
495 
496  END