LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ debchvxx()

subroutine debchvxx ( double precision  THRESH,
character*3  PATH 
)

DEBCHVXX

Purpose:
  DEBCHVXX will run D**SVXX on a series of Hilbert matrices and then
  compare the error bounds returned by D**SVXX to see if the returned
  answer indeed falls within those bounds.

  Eight test ratios will be computed.  The tests will pass if they are .LT.
  THRESH.  There are two cases that are determined by 1 / (SQRT( N ) * EPS).
  If that value is .LE. to the component wise reciprocal condition number,
  it uses the guaranteed case, other wise it uses the unguaranteed case.

  Test ratios:
     Let Xc be X_computed and Xt be X_truth.
     The norm used is the infinity norm.

     Let A be the guaranteed case and B be the unguaranteed case.

       1. Normwise guaranteed forward error bound.
       A: norm ( abs( Xc - Xt ) / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ) and
          ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N),10) * EPS.
          If these conditions are met, the test ratio is set to be
          ERRBND( *, nwise_i, bnd_i ) / MAX(SQRT(N), 10).  Otherwise it is 1/EPS.
       B: For this case, CGESVXX should just return 1.  If it is less than
          one, treat it the same as in 1A.  Otherwise it fails. (Set test
          ratio to ERRBND( *, nwise_i, bnd_i ) * THRESH?)

       2. Componentwise guaranteed forward error bound.
       A: norm ( abs( Xc(j) - Xt(j) ) ) / norm (Xt(j)) .LE. ERRBND( *, cwise_i, bnd_i )
          for all j .AND. ERRBND( *, cwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS.
          If these conditions are met, the test ratio is set to be
          ERRBND( *, cwise_i, bnd_i ) / MAX(SQRT(N), 10).  Otherwise it is 1/EPS.
       B: Same as normwise test ratio.

       3. Backwards error.
       A: The test ratio is set to BERR/EPS.
       B: Same test ratio.

       4. Reciprocal condition number.
       A: A condition number is computed with Xt and compared with the one
          returned from CGESVXX.  Let RCONDc be the RCOND returned by D**SVXX
          and RCONDt be the RCOND from the truth value.  Test ratio is set to
          MAX(RCONDc/RCONDt, RCONDt/RCONDc).
       B: Test ratio is set to 1 / (EPS * RCONDc).

       5. Reciprocal normwise condition number.
       A: The test ratio is set to
          MAX(ERRBND( *, nwise_i, cond_i ) / NCOND, NCOND / ERRBND( *, nwise_i, cond_i )).
       B: Test ratio is set to 1 / (EPS * ERRBND( *, nwise_i, cond_i )).

       6. Reciprocal componentwise condition number.
       A: Test ratio is set to
          MAX(ERRBND( *, cwise_i, cond_i ) / CCOND, CCOND / ERRBND( *, cwise_i, cond_i )).
       B: Test ratio is set to 1 / (EPS * ERRBND( *, cwise_i, cond_i )).

     .. Parameters ..
     NMAX is determined by the largest number in the inverse of the hilbert
     matrix.  Precision is exhausted when the largest entry in it is greater
     than 2 to the power of the number of bits in the fraction of the data
     type used plus one, which is 24 for single precision.
     NMAX should be 6 for single and 11 for double.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 95 of file debchvxx.f.

96  IMPLICIT NONE
97 * .. Scalar Arguments ..
98  DOUBLE PRECISION THRESH
99  CHARACTER*3 PATH
100 
101  INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU
102  parameter(nmax = 10, nparams = 2, nerrbnd = 3,
103  $ ntests = 6)
104 
105 * .. Local Scalars ..
106  INTEGER N, NRHS, INFO, I ,J, k, NFAIL, LDA,
107  $ N_AUX_TESTS, LDAB, LDAFB
108  CHARACTER FACT, TRANS, UPLO, EQUED
109  CHARACTER*2 C2
110  CHARACTER(3) NGUAR, CGUAR
111  LOGICAL printed_guide
112  DOUBLE PRECISION NCOND, CCOND, M, NORMDIF, NORMT, RCOND,
113  $ RNORM, RINORM, SUMR, SUMRI, EPS,
114  $ BERR(NMAX), RPVGRW, ORCOND,
115  $ CWISE_ERR, NWISE_ERR, CWISE_BND, NWISE_BND,
116  $ CWISE_RCOND, NWISE_RCOND,
117  $ CONDTHRESH, ERRTHRESH
118 
119 * .. Local Arrays ..
120  DOUBLE PRECISION TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
121  $ S(NMAX),R(NMAX),C(NMAX), DIFF(NMAX, NMAX),
122  $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3),
123  $ A(NMAX,NMAX),INVHILB(NMAX,NMAX),X(NMAX,NMAX),
124  $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
125  $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
126  $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX ),
127  $ WORK(NMAX*3*5), AF(NMAX, NMAX),B(NMAX, NMAX),
128  $ ACOPY(NMAX, NMAX)
129  INTEGER IPIV(NMAX), IWORK(3*NMAX)
130 
131 * .. External Functions ..
132  DOUBLE PRECISION DLAMCH
133 
134 * .. External Subroutines ..
135  EXTERNAL dlahilb, dgesvxx, dposvxx, dsysvxx,
136  $ dgbsvxx, dlacpy, lsamen
137  LOGICAL LSAMEN
138 
139 * .. Intrinsic Functions ..
140  INTRINSIC sqrt, max, abs, dble
141 
142 * .. Parameters ..
143  INTEGER NWISE_I, CWISE_I
144  parameter(nwise_i = 1, cwise_i = 1)
145  INTEGER BND_I, COND_I
146  parameter(bnd_i = 2, cond_i = 3)
147 
148 * Create the loop to test out the Hilbert matrices
149 
150  fact = 'E'
151  uplo = 'U'
152  trans = 'N'
153  equed = 'N'
154  eps = dlamch('Epsilon')
155  nfail = 0
156  n_aux_tests = 0
157  lda = nmax
158  ldab = (nmax-1)+(nmax-1)+1
159  ldafb = 2*(nmax-1)+(nmax-1)+1
160  c2 = path( 2: 3 )
161 
162 * Main loop to test the different Hilbert Matrices.
163 
164  printed_guide = .false.
165 
166  DO n = 1 , nmax
167  params(1) = -1
168  params(2) = -1
169 
170  kl = n-1
171  ku = n-1
172  nrhs = n
173  m = max(sqrt(dble(n)), 10.0d+0)
174 
175 * Generate the Hilbert matrix, its inverse, and the
176 * right hand side, all scaled by the LCM(1,..,2N-1).
177  CALL dlahilb(n, n, a, lda, invhilb, lda, b, lda, work, info)
178 
179 * Copy A into ACOPY.
180  CALL dlacpy('ALL', n, n, a, nmax, acopy, nmax)
181 
182 * Store A in band format for GB tests
183  DO j = 1, n
184  DO i = 1, kl+ku+1
185  ab( i, j ) = 0.0d+0
186  END DO
187  END DO
188  DO j = 1, n
189  DO i = max( 1, j-ku ), min( n, j+kl )
190  ab( ku+1+i-j, j ) = a( i, j )
191  END DO
192  END DO
193 
194 * Copy AB into ABCOPY.
195  DO j = 1, n
196  DO i = 1, kl+ku+1
197  abcopy( i, j ) = 0.0d+0
198  END DO
199  END DO
200  CALL dlacpy('ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)
201 
202 * Call D**SVXX with default PARAMS and N_ERR_BND = 3.
203  IF ( lsamen( 2, c2, 'SY' ) ) THEN
204  CALL dsysvxx(fact, uplo, n, nrhs, acopy, lda, af, lda,
205  $ ipiv, equed, s, b, lda, x, lda, orcond,
206  $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
207  $ params, work, iwork, info)
208  ELSE IF ( lsamen( 2, c2, 'PO' ) ) THEN
209  CALL dposvxx(fact, uplo, n, nrhs, acopy, lda, af, lda,
210  $ equed, s, b, lda, x, lda, orcond,
211  $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
212  $ params, work, iwork, info)
213  ELSE IF ( lsamen( 2, c2, 'GB' ) ) THEN
214  CALL dgbsvxx(fact, trans, n, kl, ku, nrhs, abcopy,
215  $ ldab, afb, ldafb, ipiv, equed, r, c, b,
216  $ lda, x, lda, orcond, rpvgrw, berr, nerrbnd,
217  $ errbnd_n, errbnd_c, nparams, params, work, iwork,
218  $ info)
219  ELSE
220  CALL dgesvxx(fact, trans, n, nrhs, acopy, lda, af, lda,
221  $ ipiv, equed, r, c, b, lda, x, lda, orcond,
222  $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
223  $ params, work, iwork, info)
224  END IF
225 
226  n_aux_tests = n_aux_tests + 1
227  IF (orcond .LT. eps) THEN
228 ! Either factorization failed or the matrix is flagged, and 1 <=
229 ! INFO <= N+1. We don't decide based on rcond anymore.
230 ! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN
231 ! NFAIL = NFAIL + 1
232 ! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND
233 ! END IF
234  ELSE
235 ! Either everything succeeded (INFO == 0) or some solution failed
236 ! to converge (INFO > N+1).
237  IF (info .GT. 0 .AND. info .LE. n+1) THEN
238  nfail = nfail + 1
239  WRITE (*, fmt=8000) c2, n, info, orcond, rcond
240  END IF
241  END IF
242 
243 * Calculating the difference between D**SVXX's X and the true X.
244  DO i = 1,n
245  DO j =1,nrhs
246  diff(i,j) = x(i,j) - invhilb(i,j)
247  END DO
248  END DO
249 
250 * Calculating the RCOND
251  rnorm = 0.0d+0
252  rinorm = 0.0d+0
253  IF ( lsamen( 2, c2, 'PO' ) .OR. lsamen( 2, c2, 'SY' ) ) THEN
254  DO i = 1, n
255  sumr = 0.0d+0
256  sumri = 0.0d+0
257  DO j = 1, n
258  sumr = sumr + s(i) * abs(a(i,j)) * s(j)
259  sumri = sumri + abs(invhilb(i, j)) / (s(j) * s(i))
260 
261  END DO
262  rnorm = max(rnorm,sumr)
263  rinorm = max(rinorm,sumri)
264  END DO
265  ELSE IF ( lsamen( 2, c2, 'GE' ) .OR. lsamen( 2, c2, 'GB' ) )
266  $ THEN
267  DO i = 1, n
268  sumr = 0.0d+0
269  sumri = 0.0d+0
270  DO j = 1, n
271  sumr = sumr + r(i) * abs(a(i,j)) * c(j)
272  sumri = sumri + abs(invhilb(i, j)) / (r(j) * c(i))
273  END DO
274  rnorm = max(rnorm,sumr)
275  rinorm = max(rinorm,sumri)
276  END DO
277  END IF
278 
279  rnorm = rnorm / abs(a(1, 1))
280  rcond = 1.0d+0/(rnorm * rinorm)
281 
282 * Calculating the R for normwise rcond.
283  DO i = 1, n
284  rinv(i) = 0.0d+0
285  END DO
286  DO j = 1, n
287  DO i = 1, n
288  rinv(i) = rinv(i) + abs(a(i,j))
289  END DO
290  END DO
291 
292 * Calculating the Normwise rcond.
293  rinorm = 0.0d+0
294  DO i = 1, n
295  sumri = 0.0d+0
296  DO j = 1, n
297  sumri = sumri + abs(invhilb(i,j) * rinv(j))
298  END DO
299  rinorm = max(rinorm, sumri)
300  END DO
301 
302 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
303 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
304  ncond = abs(a(1,1)) / rinorm
305 
306  condthresh = m * eps
307  errthresh = m * eps
308 
309  DO k = 1, nrhs
310  normt = 0.0d+0
311  normdif = 0.0d+0
312  cwise_err = 0.0d+0
313  DO i = 1, n
314  normt = max(abs(invhilb(i, k)), normt)
315  normdif = max(abs(x(i,k) - invhilb(i,k)), normdif)
316  IF (invhilb(i,k) .NE. 0.0d+0) THEN
317  cwise_err = max(abs(x(i,k) - invhilb(i,k))
318  $ /abs(invhilb(i,k)), cwise_err)
319  ELSE IF (x(i, k) .NE. 0.0d+0) THEN
320  cwise_err = dlamch('OVERFLOW')
321  END IF
322  END DO
323  IF (normt .NE. 0.0d+0) THEN
324  nwise_err = normdif / normt
325  ELSE IF (normdif .NE. 0.0d+0) THEN
326  nwise_err = dlamch('OVERFLOW')
327  ELSE
328  nwise_err = 0.0d+0
329  ENDIF
330 
331  DO i = 1, n
332  rinv(i) = 0.0d+0
333  END DO
334  DO j = 1, n
335  DO i = 1, n
336  rinv(i) = rinv(i) + abs(a(i, j) * invhilb(j, k))
337  END DO
338  END DO
339  rinorm = 0.0d+0
340  DO i = 1, n
341  sumri = 0.0d+0
342  DO j = 1, n
343  sumri = sumri
344  $ + abs(invhilb(i, j) * rinv(j) / invhilb(i, k))
345  END DO
346  rinorm = max(rinorm, sumri)
347  END DO
348 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
349 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
350  ccond = abs(a(1,1))/rinorm
351 
352 ! Forward error bound tests
353  nwise_bnd = errbnd_n(k + (bnd_i-1)*nrhs)
354  cwise_bnd = errbnd_c(k + (bnd_i-1)*nrhs)
355  nwise_rcond = errbnd_n(k + (cond_i-1)*nrhs)
356  cwise_rcond = errbnd_c(k + (cond_i-1)*nrhs)
357 ! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
358 ! $ condthresh, ncond.ge.condthresh
359 ! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
360  IF (ncond .GE. condthresh) THEN
361  nguar = 'YES'
362  IF (nwise_bnd .GT. errthresh) THEN
363  tstrat(1) = 1/(2.0d+0*eps)
364  ELSE
365  IF (nwise_bnd .NE. 0.0d+0) THEN
366  tstrat(1) = nwise_err / nwise_bnd
367  ELSE IF (nwise_err .NE. 0.0d+0) THEN
368  tstrat(1) = 1/(16.0*eps)
369  ELSE
370  tstrat(1) = 0.0d+0
371  END IF
372  IF (tstrat(1) .GT. 1.0d+0) THEN
373  tstrat(1) = 1/(4.0d+0*eps)
374  END IF
375  END IF
376  ELSE
377  nguar = 'NO'
378  IF (nwise_bnd .LT. 1.0d+0) THEN
379  tstrat(1) = 1/(8.0d+0*eps)
380  ELSE
381  tstrat(1) = 1.0d+0
382  END IF
383  END IF
384 ! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
385 ! $ condthresh, ccond.ge.condthresh
386 ! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
387  IF (ccond .GE. condthresh) THEN
388  cguar = 'YES'
389  IF (cwise_bnd .GT. errthresh) THEN
390  tstrat(2) = 1/(2.0d+0*eps)
391  ELSE
392  IF (cwise_bnd .NE. 0.0d+0) THEN
393  tstrat(2) = cwise_err / cwise_bnd
394  ELSE IF (cwise_err .NE. 0.0d+0) THEN
395  tstrat(2) = 1/(16.0d+0*eps)
396  ELSE
397  tstrat(2) = 0.0d+0
398  END IF
399  IF (tstrat(2) .GT. 1.0d+0) tstrat(2) = 1/(4.0d+0*eps)
400  END IF
401  ELSE
402  cguar = 'NO'
403  IF (cwise_bnd .LT. 1.0d+0) THEN
404  tstrat(2) = 1/(8.0d+0*eps)
405  ELSE
406  tstrat(2) = 1.0d+0
407  END IF
408  END IF
409 
410 ! Backwards error test
411  tstrat(3) = berr(k)/eps
412 
413 ! Condition number tests
414  tstrat(4) = rcond / orcond
415  IF (rcond .GE. condthresh .AND. tstrat(4) .LT. 1.0d+0)
416  $ tstrat(4) = 1.0d+0 / tstrat(4)
417 
418  tstrat(5) = ncond / nwise_rcond
419  IF (ncond .GE. condthresh .AND. tstrat(5) .LT. 1.0d+0)
420  $ tstrat(5) = 1.0d+0 / tstrat(5)
421 
422  tstrat(6) = ccond / nwise_rcond
423  IF (ccond .GE. condthresh .AND. tstrat(6) .LT. 1.0d+0)
424  $ tstrat(6) = 1.0d+0 / tstrat(6)
425 
426  DO i = 1, ntests
427  IF (tstrat(i) .GT. thresh) THEN
428  IF (.NOT.printed_guide) THEN
429  WRITE(*,*)
430  WRITE( *, 9996) 1
431  WRITE( *, 9995) 2
432  WRITE( *, 9994) 3
433  WRITE( *, 9993) 4
434  WRITE( *, 9992) 5
435  WRITE( *, 9991) 6
436  WRITE( *, 9990) 7
437  WRITE( *, 9989) 8
438  WRITE(*,*)
439  printed_guide = .true.
440  END IF
441  WRITE( *, 9999) c2, n, k, nguar, cguar, i, tstrat(i)
442  nfail = nfail + 1
443  END IF
444  END DO
445  END DO
446 
447 c$$$ WRITE(*,*)
448 c$$$ WRITE(*,*) 'Normwise Error Bounds'
449 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i)
450 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i)
451 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i)
452 c$$$ WRITE(*,*)
453 c$$$ WRITE(*,*) 'Componentwise Error Bounds'
454 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i)
455 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i)
456 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i)
457 c$$$ print *, 'Info: ', info
458 c$$$ WRITE(*,*)
459 * WRITE(*,*) 'TSTRAT: ',TSTRAT
460 
461  END DO
462 
463  WRITE(*,*)
464  IF( nfail .GT. 0 ) THEN
465  WRITE(*,9998) c2, nfail, ntests*n+n_aux_tests
466  ELSE
467  WRITE(*,9997) c2
468  END IF
469  9999 FORMAT( ' D', a2, 'SVXX: N =', i2, ', RHS = ', i2,
470  $ ', NWISE GUAR. = ', a, ', CWISE GUAR. = ', a,
471  $ ' test(',i1,') =', g12.5 )
472  9998 FORMAT( ' D', a2, 'SVXX: ', i6, ' out of ', i6,
473  $ ' tests failed to pass the threshold' )
474  9997 FORMAT( ' D', a2, 'SVXX passed the tests of error bounds' )
475 * Test ratios.
476  9996 FORMAT( 3x, i2, ': Normwise guaranteed forward error', / 5x,
477  $ 'Guaranteed case: if norm ( abs( Xc - Xt )',
478  $ .LE.' / norm ( Xt ) ERRBND( *, nwise_i, bnd_i ), then',
479  $ / 5x,
480  $ .LE.'ERRBND( *, nwise_i, bnd_i ) MAX(SQRT(N), 10) * EPS')
481  9995 FORMAT( 3x, i2, ': Componentwise guaranteed forward error' )
482  9994 FORMAT( 3x, i2, ': Backwards error' )
483  9993 FORMAT( 3x, i2, ': Reciprocal condition number' )
484  9992 FORMAT( 3x, i2, ': Reciprocal normwise condition number' )
485  9991 FORMAT( 3x, i2, ': Raw normwise error estimate' )
486  9990 FORMAT( 3x, i2, ': Reciprocal componentwise condition number' )
487  9989 FORMAT( 3x, i2, ': Raw componentwise error estimate' )
488 
489  8000 FORMAT( ' D', a2, 'SVXX: N =', i2, ', INFO = ', i3,
490  $ ', ORCOND = ', g12.5, ', real RCOND = ', g12.5 )
491 *
492 * End of DEBCHVXX
493 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
logical function lsamen(N, CA, CB)
LSAMEN
Definition: lsamen.f:74
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
DLAHILB
Definition: dlahilb.f:124
subroutine dgbsvxx(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, IWORK, INFO)
DGBSVXX computes the solution to system of linear equations A * X = B for GB matrices
Definition: dgbsvxx.f:560
subroutine dgesvxx(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, IWORK, INFO)
DGESVXX computes the solution to system of linear equations A * X = B for GE matrices
Definition: dgesvxx.f:540
subroutine dposvxx(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, IWORK, INFO)
DPOSVXX computes the solution to system of linear equations A * X = B for PO matrices
Definition: dposvxx.f:494
subroutine dsysvxx(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, IWORK, INFO)
DSYSVXX
Definition: dsysvxx.f:505
Here is the call graph for this function: