LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zebchvxx.f
Go to the documentation of this file.
1 *> \brief \b ZEBCHVXX
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 ZEBCHVXX( THRESH, PATH )
12 *
13 * .. Scalar Arguments ..
14 * DOUBLE PRECISION THRESH
15 * CHARACTER*3 PATH
16 * ..
17 *
18 * Purpose
19 * ======
20 *
21 *> \details \b Purpose:
22 *> \verbatim
23 *>
24 *> ZEBCHVXX will run Z**SVXX on a series of Hilbert matrices and then
25 *> compare the error bounds returned by Z**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 CGESVXX
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 complex16_lin
95 *
96 * =====================================================================
97  SUBROUTINE zebchvxx( 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  COMPLEX*16 zdum
121 
122 * .. Local Arrays ..
123  DOUBLE PRECISION tstrat(ntests), rinv(nmax), params(nparams),
124  $ s(nmax),r(nmax),c(nmax),rwork(3*nmax),
125  $ diff(nmax, nmax),
126  $ errbnd_n(nmax*3), errbnd_c(nmax*3)
127  INTEGER ipiv(nmax)
128  COMPLEX*16 a(nmax,nmax),invhilb(nmax,nmax),x(nmax,nmax),
129  $ work(nmax*3*5), af(nmax, nmax),b(nmax, nmax),
130  $ acopy(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 )
134 
135 * .. External Functions ..
136  DOUBLE PRECISION dlamch
137 
138 * .. External Subroutines ..
139  EXTERNAL zlahilb, zgesvxx, zposvxx, zsysvxx,
140  $ zgbsvxx, zlacpy, lsamen
141  LOGICAL lsamen
142 
143 * .. Intrinsic Functions ..
144  INTRINSIC sqrt, max, abs, dble, dimag
145 
146 * .. Statement Functions ..
147  DOUBLE PRECISION cabs1
148 
149 * .. Statement Function Definitions ..
150  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
151 
152 * .. Parameters ..
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)
157 
158 * Create the loop to test out the Hilbert matrices
159 
160  fact = 'E'
161  uplo = 'U'
162  trans = 'N'
163  equed = 'N'
164  eps = dlamch('Epsilon')
165  nfail = 0
166  n_aux_tests = 0
167  lda = nmax
168  ldab = (nmax-1)+(nmax-1)+1
169  ldafb = 2*(nmax-1)+(nmax-1)+1
170  c2 = path( 2: 3 )
171 
172 * Main loop to test the different Hilbert Matrices.
173 
174  printed_guide = .false.
175 
176  DO n = 1 , nmax
177  params(1) = -1
178  params(2) = -1
179 
180  kl = n-1
181  ku = n-1
182  nrhs = n
183  m = max(sqrt(dble(n)), 10.0d+0)
184 
185 * Generate the Hilbert matrix, its inverse, and the
186 * right hand side, all scaled by the LCM(1,..,2N-1).
187  CALL zlahilb(n, n, a, lda, invhilb, lda, b,
188  $ lda, work, info, path)
189 
190 * Copy A into ACOPY.
191  CALL zlacpy('ALL', n, n, a, nmax, acopy, nmax)
192 
193 * Store A in band format for GB tests
194  DO j = 1, n
195  DO i = 1, kl+ku+1
196  ab( i, j ) = (0.0d+0,0.0d+0)
197  END DO
198  END DO
199  DO j = 1, n
200  DO i = max( 1, j-ku ), min( n, j+kl )
201  ab( ku+1+i-j, j ) = a( i, j )
202  END DO
203  END DO
204 
205 * Copy AB into ABCOPY.
206  DO j = 1, n
207  DO i = 1, kl+ku+1
208  abcopy( i, j ) = (0.0d+0,0.0d+0)
209  END DO
210  END DO
211  CALL zlacpy('ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)
212 
213 * Call Z**SVXX with default PARAMS and N_ERR_BND = 3.
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,
234  $ info)
235  ELSE
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)
240  END IF
241 
242  n_aux_tests = n_aux_tests + 1
243  IF (orcond .LT. eps) THEN
244 ! Either factorization failed or the matrix is flagged, and 1 <=
245 ! INFO <= N+1. We don't decide based on rcond anymore.
246 ! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN
247 ! NFAIL = NFAIL + 1
248 ! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND
249 ! END IF
250  ELSE
251 ! Either everything succeeded (INFO == 0) or some solution failed
252 ! to converge (INFO > N+1).
253  IF (info .GT. 0 .AND. info .LE. n+1) THEN
254  nfail = nfail + 1
255  WRITE (*, fmt=8000) c2, n, info, orcond, rcond
256  END IF
257  END IF
258 
259 * Calculating the difference between Z**SVXX's X and the true X.
260  DO i = 1,n
261  DO j =1,nrhs
262  diff(i,j) = x(i,j) - invhilb(i,j)
263  END DO
264  END DO
265 
266 * Calculating the RCOND
267  rnorm = 0
268  rinorm = 0
269  IF ( lsamen( 2, c2, 'PO' ) .OR. lsamen( 2, c2, 'SY' ) .OR.
270  $ lsamen( 2, c2, 'HE' ) ) THEN
271  DO i = 1, n
272  sumr = 0
273  sumri = 0
274  DO j = 1, n
275  sumr = sumr + s(i) * cabs1(a(i,j)) * s(j)
276  sumri = sumri + cabs1(invhilb(i, j)) / (s(j) * s(i))
277  END DO
278  rnorm = max(rnorm,sumr)
279  rinorm = max(rinorm,sumri)
280  END DO
281  ELSE IF ( lsamen( 2, c2, 'GE' ) .OR. lsamen( 2, c2, 'GB' ) )
282  $ THEN
283  DO i = 1, n
284  sumr = 0
285  sumri = 0
286  DO j = 1, n
287  sumr = sumr + r(i) * cabs1(a(i,j)) * c(j)
288  sumri = sumri + cabs1(invhilb(i, j)) / (r(j) * c(i))
289  END DO
290  rnorm = max(rnorm,sumr)
291  rinorm = max(rinorm,sumri)
292  END DO
293  END IF
294 
295  rnorm = rnorm / cabs1(a(1, 1))
296  rcond = 1.0d+0/(rnorm * rinorm)
297 
298 * Calculating the R for normwise rcond.
299  DO i = 1, n
300  rinv(i) = 0.0d+0
301  END DO
302  DO j = 1, n
303  DO i = 1, n
304  rinv(i) = rinv(i) + cabs1(a(i,j))
305  END DO
306  END DO
307 
308 * Calculating the Normwise rcond.
309  rinorm = 0.0d+0
310  DO i = 1, n
311  sumri = 0.0d+0
312  DO j = 1, n
313  sumri = sumri + cabs1(invhilb(i,j) * rinv(j))
314  END DO
315  rinorm = max(rinorm, sumri)
316  END DO
317 
318 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
319 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
320  ncond = cabs1(a(1,1)) / rinorm
321 
322  condthresh = m * eps
323  errthresh = m * eps
324 
325  DO k = 1, nrhs
326  normt = 0.0d+0
327  normdif = 0.0d+0
328  cwise_err = 0.0d+0
329  DO i = 1, n
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')
337  END IF
338  END DO
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')
343  ELSE
344  nwise_err = 0.0d+0
345  ENDIF
346 
347  DO i = 1, n
348  rinv(i) = 0.0d+0
349  END DO
350  DO j = 1, n
351  DO i = 1, n
352  rinv(i) = rinv(i) + cabs1(a(i, j) * invhilb(j, k))
353  END DO
354  END DO
355  rinorm = 0.0d+0
356  DO i = 1, n
357  sumri = 0.0d+0
358  DO j = 1, n
359  sumri = sumri
360  $ + cabs1(invhilb(i, j) * rinv(j) / invhilb(i, k))
361  END DO
362  rinorm = max(rinorm, sumri)
363  END DO
364 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
365 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
366  ccond = cabs1(a(1,1))/rinorm
367 
368 ! Forward error bound tests
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)
373 ! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
374 ! $ condthresh, ncond.ge.condthresh
375 ! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
376  IF (ncond .GE. condthresh) THEN
377  nguar = 'YES'
378  IF (nwise_bnd .GT. errthresh) THEN
379  tstrat(1) = 1/(2.0d+0*eps)
380  ELSE
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)
385  ELSE
386  tstrat(1) = 0.0d+0
387  END IF
388  IF (tstrat(1) .GT. 1.0d+0) THEN
389  tstrat(1) = 1/(4.0d+0*eps)
390  END IF
391  END IF
392  ELSE
393  nguar = 'NO'
394  IF (nwise_bnd .LT. 1.0d+0) THEN
395  tstrat(1) = 1/(8.0d+0*eps)
396  ELSE
397  tstrat(1) = 1.0d+0
398  END IF
399  END IF
400 ! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
401 ! $ condthresh, ccond.ge.condthresh
402 ! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
403  IF (ccond .GE. condthresh) THEN
404  cguar = 'YES'
405  IF (cwise_bnd .GT. errthresh) THEN
406  tstrat(2) = 1/(2.0d+0*eps)
407  ELSE
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)
412  ELSE
413  tstrat(2) = 0.0d+0
414  END IF
415  IF (tstrat(2) .GT. 1.0d+0) tstrat(2) = 1/(4.0d+0*eps)
416  END IF
417  ELSE
418  cguar = 'NO'
419  IF (cwise_bnd .LT. 1.0d+0) THEN
420  tstrat(2) = 1/(8.0d+0*eps)
421  ELSE
422  tstrat(2) = 1.0d+0
423  END IF
424  END IF
425 
426 ! Backwards error test
427  tstrat(3) = berr(k)/eps
428 
429 ! Condition number tests
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)
433 
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)
437 
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)
441 
442  DO i = 1, ntests
443  IF (tstrat(i) .GT. thresh) THEN
444  IF (.NOT.printed_guide) THEN
445  WRITE(*,*)
446  WRITE( *, 9996) 1
447  WRITE( *, 9995) 2
448  WRITE( *, 9994) 3
449  WRITE( *, 9993) 4
450  WRITE( *, 9992) 5
451  WRITE( *, 9991) 6
452  WRITE( *, 9990) 7
453  WRITE( *, 9989) 8
454  WRITE(*,*)
455  printed_guide = .true.
456  END IF
457  WRITE( *, 9999) c2, n, k, nguar, cguar, i, tstrat(i)
458  nfail = nfail + 1
459  END IF
460  END DO
461  END DO
462 
463 c$$$ WRITE(*,*)
464 c$$$ WRITE(*,*) 'Normwise Error Bounds'
465 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i)
466 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i)
467 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i)
468 c$$$ WRITE(*,*)
469 c$$$ WRITE(*,*) 'Componentwise Error Bounds'
470 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i)
471 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i)
472 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i)
473 c$$$ print *, 'Info: ', info
474 c$$$ WRITE(*,*)
475 * WRITE(*,*) 'TSTRAT: ',TSTRAT
476 
477  END DO
478 
479  WRITE(*,*)
480  IF( nfail .GT. 0 ) THEN
481  WRITE(*,9998) c2, nfail, ntests*n+n_aux_tests
482  ELSE
483  WRITE(*,9997) c2
484  END IF
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' )
491 * Test ratios.
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',
495  $ / 5x,
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' )
504 
505  8000 format( ' Z', a2, 'SVXX: N =', i2, ', INFO = ', i3,
506  $ ', ORCOND = ', g12.5, ', real RCOND = ', g12.5 )
507 
508  END