LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ sebchvxx()

subroutine sebchvxx ( real  THRESH,
character*3  PATH 
)

SEBCHVXX

Purpose:
  SEBCHVXX will run S**SVXX on a series of Hilbert matrices and then
  compare the error bounds returned by SGESVXX 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, SGESVXX 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 SGESVXX.  Let RCONDc be the RCOND returned by SGESVXX
          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 )).

       7. 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 sebchvxx.f.

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 
subroutine sposvxx(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)
SPOSVXX computes the solution to system of linear equations A * X = B for PO matrices ...
Definition: sposvxx.f:499
subroutine slahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
SLAHILB
Definition: slahilb.f:126
logical function lsamen(N, CA, CB)
LSAMEN
Definition: lsamen.f:76
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine sgesvxx(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)
SGESVXX computes the solution to system of linear equations A * X = B for GE matrices ...
Definition: sgesvxx.f:545
subroutine ssysvxx(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)
SSYSVXX
Definition: ssysvxx.f:510
subroutine sgbsvxx(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)
SGBSVXX computes the solution to system of linear equations A * X = B for GB matrices ...
Definition: sgbsvxx.f:565
Here is the call graph for this function: