LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
debchvxx.f
Go to the documentation of this file.
1 *> \brief \b DEBCHVXX
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 DEBCHVXX( THRESH, PATH )
12 *
13 * .. Scalar Arguments ..
14 * DOUBLE PRECISION THRESH
15 * CHARACTER*3 PATH
16 * ..
17 *
18 *
19 *> \par Purpose:
20 * =============
21 *>
22 *> \verbatim
23 *>
24 *> DEBCHVXX will run D**SVXX on a series of Hilbert matrices and then
25 *> compare the error bounds returned by D**SVXX 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, CGESVXX 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 CGESVXX. Let RCONDc be the RCOND returned by D**SVXX
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 *> 6. 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 double_lin
95 *
96 * =====================================================================
97  SUBROUTINE debchvxx( THRESH, PATH )
98  IMPLICIT NONE
99 * .. Scalar Arguments ..
100  DOUBLE PRECISION thresh
101  CHARACTER*3 path
102 
103  INTEGER nmax, nparams, nerrbnd, ntests, kl, ku
104  parameter(nmax = 10, nparams = 2, nerrbnd = 3,
105  $ ntests = 6)
106 
107 * .. Local Scalars ..
108  INTEGER n, nrhs, info, i ,j, k, nfail, lda,
109  $ n_aux_tests, ldab, ldafb
110  CHARACTER fact, trans, uplo, equed
111  CHARACTER*2 c2
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
120 
121 * .. Local Arrays ..
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),
130  $ acopy(nmax, nmax)
131  INTEGER ipiv(nmax), iwork(3*nmax)
132 
133 * .. External Functions ..
134  DOUBLE PRECISION dlamch
135 
136 * .. External Subroutines ..
137  EXTERNAL dlahilb, dgesvxx, dposvxx, dsysvxx,
138  $ dgbsvxx, dlacpy, lsamen
139  LOGICAL lsamen
140 
141 * .. Intrinsic Functions ..
142  INTRINSIC sqrt, max, abs, dble
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 = dlamch('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(dble(n)), 10.0d+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 dlahilb(n, n, a, lda, invhilb, lda, b, lda, work, info)
180 
181 * Copy A into ACOPY.
182  CALL dlacpy('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.0d+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.0d+0
200  END DO
201  END DO
202  CALL dlacpy('ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)
203 
204 * Call D**SVXX with default PARAMS and N_ERR_BND = 3.
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,
220  $ info)
221  ELSE
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)
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 D**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.0d+0
254  rinorm = 0.0d+0
255  IF ( lsamen( 2, c2, 'PO' ) .OR. lsamen( 2, c2, 'SY' ) ) THEN
256  DO i = 1, n
257  sumr = 0.0d+0
258  sumri = 0.0d+0
259  DO j = 1, n
260  sumr = sumr + s(i) * abs(a(i,j)) * s(j)
261  sumri = sumri + abs(invhilb(i, j)) / (s(j) * s(i))
262 
263  END DO
264  rnorm = max(rnorm,sumr)
265  rinorm = max(rinorm,sumri)
266  END DO
267  ELSE IF ( lsamen( 2, c2, 'GE' ) .OR. lsamen( 2, c2, 'GB' ) )
268  $ THEN
269  DO i = 1, n
270  sumr = 0.0d+0
271  sumri = 0.0d+0
272  DO j = 1, n
273  sumr = sumr + r(i) * abs(a(i,j)) * c(j)
274  sumri = sumri + abs(invhilb(i, j)) / (r(j) * c(i))
275  END DO
276  rnorm = max(rnorm,sumr)
277  rinorm = max(rinorm,sumri)
278  END DO
279  END IF
280 
281  rnorm = rnorm / abs(a(1, 1))
282  rcond = 1.0d+0/(rnorm * rinorm)
283 
284 * Calculating the R for normwise rcond.
285  DO i = 1, n
286  rinv(i) = 0.0d+0
287  END DO
288  DO j = 1, n
289  DO i = 1, n
290  rinv(i) = rinv(i) + abs(a(i,j))
291  END DO
292  END DO
293 
294 * Calculating the Normwise rcond.
295  rinorm = 0.0d+0
296  DO i = 1, n
297  sumri = 0.0d+0
298  DO j = 1, n
299  sumri = sumri + abs(invhilb(i,j) * rinv(j))
300  END DO
301  rinorm = max(rinorm, sumri)
302  END DO
303 
304 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
305 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
306  ncond = abs(a(1,1)) / rinorm
307 
308  condthresh = m * eps
309  errthresh = m * eps
310 
311  DO k = 1, nrhs
312  normt = 0.0d+0
313  normdif = 0.0d+0
314  cwise_err = 0.0d+0
315  DO i = 1, n
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')
323  END IF
324  END DO
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')
329  ELSE
330  nwise_err = 0.0d+0
331  ENDIF
332 
333  DO i = 1, n
334  rinv(i) = 0.0d+0
335  END DO
336  DO j = 1, n
337  DO i = 1, n
338  rinv(i) = rinv(i) + abs(a(i, j) * invhilb(j, k))
339  END DO
340  END DO
341  rinorm = 0.0d+0
342  DO i = 1, n
343  sumri = 0.0d+0
344  DO j = 1, n
345  sumri = sumri
346  $ + abs(invhilb(i, j) * rinv(j) / invhilb(i, k))
347  END DO
348  rinorm = max(rinorm, sumri)
349  END DO
350 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
351 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
352  ccond = abs(a(1,1))/rinorm
353 
354 ! Forward error bound tests
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)
359 ! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
360 ! $ condthresh, ncond.ge.condthresh
361 ! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
362  IF (ncond .GE. condthresh) THEN
363  nguar = 'YES'
364  IF (nwise_bnd .GT. errthresh) THEN
365  tstrat(1) = 1/(2.0d+0*eps)
366  ELSE
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)
371  ELSE
372  tstrat(1) = 0.0d+0
373  END IF
374  IF (tstrat(1) .GT. 1.0d+0) THEN
375  tstrat(1) = 1/(4.0d+0*eps)
376  END IF
377  END IF
378  ELSE
379  nguar = 'NO'
380  IF (nwise_bnd .LT. 1.0d+0) THEN
381  tstrat(1) = 1/(8.0d+0*eps)
382  ELSE
383  tstrat(1) = 1.0d+0
384  END IF
385  END IF
386 ! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
387 ! $ condthresh, ccond.ge.condthresh
388 ! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
389  IF (ccond .GE. condthresh) THEN
390  cguar = 'YES'
391  IF (cwise_bnd .GT. errthresh) THEN
392  tstrat(2) = 1/(2.0d+0*eps)
393  ELSE
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)
398  ELSE
399  tstrat(2) = 0.0d+0
400  END IF
401  IF (tstrat(2) .GT. 1.0d+0) tstrat(2) = 1/(4.0d+0*eps)
402  END IF
403  ELSE
404  cguar = 'NO'
405  IF (cwise_bnd .LT. 1.0d+0) THEN
406  tstrat(2) = 1/(8.0d+0*eps)
407  ELSE
408  tstrat(2) = 1.0d+0
409  END IF
410  END IF
411 
412 ! Backwards error test
413  tstrat(3) = berr(k)/eps
414 
415 ! Condition number tests
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)
419 
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)
423 
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)
427 
428  DO i = 1, ntests
429  IF (tstrat(i) .GT. thresh) THEN
430  IF (.NOT.printed_guide) THEN
431  WRITE(*,*)
432  WRITE( *, 9996) 1
433  WRITE( *, 9995) 2
434  WRITE( *, 9994) 3
435  WRITE( *, 9993) 4
436  WRITE( *, 9992) 5
437  WRITE( *, 9991) 6
438  WRITE( *, 9990) 7
439  WRITE( *, 9989) 8
440  WRITE(*,*)
441  printed_guide = .true.
442  END IF
443  WRITE( *, 9999) c2, n, k, nguar, cguar, i, tstrat(i)
444  nfail = nfail + 1
445  END IF
446  END DO
447  END DO
448 
449 c$$$ WRITE(*,*)
450 c$$$ WRITE(*,*) 'Normwise Error Bounds'
451 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i)
452 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i)
453 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i)
454 c$$$ WRITE(*,*)
455 c$$$ WRITE(*,*) 'Componentwise Error Bounds'
456 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i)
457 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i)
458 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i)
459 c$$$ print *, 'Info: ', info
460 c$$$ WRITE(*,*)
461 * WRITE(*,*) 'TSTRAT: ',TSTRAT
462 
463  END DO
464 
465  WRITE(*,*)
466  IF( nfail .GT. 0 ) THEN
467  WRITE(*,9998) c2, nfail, ntests*n+n_aux_tests
468  ELSE
469  WRITE(*,9997) c2
470  END IF
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' )
477 * Test ratios.
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',
481  $ / 5x,
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' )
490 
491  8000 format( ' D', a2, 'SVXX: N =', i2, ', INFO = ', i3,
492  $ ', ORCOND = ', g12.5, ', real RCOND = ', g12.5 )
493 
494  END