LAPACK  3.8.0
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.
Date
December 2016

Definition at line 98 of file debchvxx.f.

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 
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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:562
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:507
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:496
subroutine dlahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
DLAHILB
Definition: dlahilb.f:126
logical function lsamen(N, CA, CB)
LSAMEN
Definition: lsamen.f:76
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:542
Here is the call graph for this function: