LAPACK  3.10.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.

Definition at line 95 of file sebchvxx.f.

96  IMPLICIT NONE
97 * .. Scalar Arguments ..
98  REAL THRESH
99  CHARACTER*3 PATH
100 
101  INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU
102  parameter(nmax = 6, nparams = 2, nerrbnd = 3,
103  $ ntests = 6)
104 
105 * .. Local Scalars ..
106  INTEGER N, NRHS, INFO, I ,J, k, NFAIL, LDA, LDAB,
107  $ LDAFB, N_AUX_TESTS
108  CHARACTER FACT, TRANS, UPLO, EQUED
109  CHARACTER*2 C2
110  CHARACTER(3) NGUAR, CGUAR
111  LOGICAL printed_guide
112  REAL 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  REAL TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
121  $ A(NMAX, NMAX), ACOPY(NMAX, NMAX),
122  $ INVHILB(NMAX, NMAX), R(NMAX), C(NMAX), S(NMAX),
123  $ WORK(NMAX * 5), B(NMAX, NMAX), X(NMAX, NMAX),
124  $ DIFF(NMAX, NMAX), AF(NMAX, NMAX),
125  $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
126  $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
127  $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX ),
128  $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3)
129  INTEGER IWORK(NMAX), IPIV(NMAX)
130 
131 * .. External Functions ..
132  REAL SLAMCH
133 
134 * .. External Subroutines ..
135  EXTERNAL slahilb, sgesvxx, ssysvxx, sposvxx, sgbsvxx,
136  $ slacpy, lsamen
137  LOGICAL LSAMEN
138 
139 * .. Intrinsic Functions ..
140  INTRINSIC sqrt, max, abs
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 = slamch('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(real(n)), 10.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 slahilb(n, n, a, lda, invhilb, lda, b, lda, work, info)
178 
179 * Copy A into ACOPY.
180  CALL slacpy('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.0e+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.0e+0
198  END DO
199  END DO
200  CALL slacpy('ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)
201 
202 * Call S**SVXX with default PARAMS and N_ERR_BND = 3.
203  IF ( lsamen( 2, c2, 'SY' ) ) THEN
204  CALL ssysvxx(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 sposvxx(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 sgbsvxx(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 sgesvxx(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 S**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
252  rinorm = 0
253  IF ( lsamen( 2, c2, 'PO' ) .OR. lsamen( 2, c2, 'SY' ) ) THEN
254  DO i = 1, n
255  sumr = 0
256  sumri = 0
257  DO j = 1, n
258  sumr = sumr + abs(s(i) * a(i,j) * s(j))
259  sumri = sumri + abs(invhilb(i, j) / s(j) / s(i))
260  END DO
261  rnorm = max(rnorm,sumr)
262  rinorm = max(rinorm,sumri)
263  END DO
264  ELSE IF ( lsamen( 2, c2, 'GE' ) .OR. lsamen( 2, c2, 'GB' ) )
265  $ THEN
266  DO i = 1, n
267  sumr = 0
268  sumri = 0
269  DO j = 1, n
270  sumr = sumr + abs(r(i) * a(i,j) * c(j))
271  sumri = sumri + abs(invhilb(i, j) / r(j) / c(i))
272  END DO
273  rnorm = max(rnorm,sumr)
274  rinorm = max(rinorm,sumri)
275  END DO
276  END IF
277 
278  rnorm = rnorm / a(1, 1)
279  rcond = 1.0/(rnorm * rinorm)
280 
281 * Calculating the R for normwise rcond.
282  DO i = 1, n
283  rinv(i) = 0.0
284  END DO
285  DO j = 1, n
286  DO i = 1, n
287  rinv(i) = rinv(i) + abs(a(i,j))
288  END DO
289  END DO
290 
291 * Calculating the Normwise rcond.
292  rinorm = 0.0
293  DO i = 1, n
294  sumri = 0.0
295  DO j = 1, n
296  sumri = sumri + abs(invhilb(i,j) * rinv(j))
297  END DO
298  rinorm = max(rinorm, sumri)
299  END DO
300 
301 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
302 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
303  ncond = a(1,1) / rinorm
304 
305  condthresh = m * eps
306  errthresh = m * eps
307 
308  DO k = 1, nrhs
309  normt = 0.0
310  normdif = 0.0
311  cwise_err = 0.0
312  DO i = 1, n
313  normt = max(abs(invhilb(i, k)), normt)
314  normdif = max(abs(x(i,k) - invhilb(i,k)), normdif)
315  IF (invhilb(i,k) .NE. 0.0) THEN
316  cwise_err = max(abs(x(i,k) - invhilb(i,k))
317  $ /abs(invhilb(i,k)), cwise_err)
318  ELSE IF (x(i, k) .NE. 0.0) THEN
319  cwise_err = slamch('OVERFLOW')
320  END IF
321  END DO
322  IF (normt .NE. 0.0) THEN
323  nwise_err = normdif / normt
324  ELSE IF (normdif .NE. 0.0) THEN
325  nwise_err = slamch('OVERFLOW')
326  ELSE
327  nwise_err = 0.0
328  ENDIF
329 
330  DO i = 1, n
331  rinv(i) = 0.0
332  END DO
333  DO j = 1, n
334  DO i = 1, n
335  rinv(i) = rinv(i) + abs(a(i, j) * invhilb(j, k))
336  END DO
337  END DO
338  rinorm = 0.0
339  DO i = 1, n
340  sumri = 0.0
341  DO j = 1, n
342  sumri = sumri
343  $ + abs(invhilb(i, j) * rinv(j) / invhilb(i, k))
344  END DO
345  rinorm = max(rinorm, sumri)
346  END DO
347 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
348 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
349  ccond = a(1,1)/rinorm
350 
351 ! Forward error bound tests
352  nwise_bnd = errbnd_n(k + (bnd_i-1)*nrhs)
353  cwise_bnd = errbnd_c(k + (bnd_i-1)*nrhs)
354  nwise_rcond = errbnd_n(k + (cond_i-1)*nrhs)
355  cwise_rcond = errbnd_c(k + (cond_i-1)*nrhs)
356 ! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
357 ! $ condthresh, ncond.ge.condthresh
358 ! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
359 
360  IF (ncond .GE. condthresh) THEN
361  nguar = 'YES'
362  IF (nwise_bnd .GT. errthresh) THEN
363  tstrat(1) = 1/(2.0*eps)
364  ELSE
365 
366  IF (nwise_bnd .NE. 0.0) THEN
367  tstrat(1) = nwise_err / nwise_bnd
368  ELSE IF (nwise_err .NE. 0.0) THEN
369  tstrat(1) = 1/(16.0*eps)
370  ELSE
371  tstrat(1) = 0.0
372  END IF
373  IF (tstrat(1) .GT. 1.0) THEN
374  tstrat(1) = 1/(4.0*eps)
375  END IF
376  END IF
377  ELSE
378  nguar = 'NO'
379  IF (nwise_bnd .LT. 1.0) THEN
380  tstrat(1) = 1/(8.0*eps)
381  ELSE
382  tstrat(1) = 1.0
383  END IF
384  END IF
385 ! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
386 ! $ condthresh, ccond.ge.condthresh
387 ! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
388  IF (ccond .GE. condthresh) THEN
389  cguar = 'YES'
390 
391  IF (cwise_bnd .GT. errthresh) THEN
392  tstrat(2) = 1/(2.0*eps)
393  ELSE
394  IF (cwise_bnd .NE. 0.0) THEN
395  tstrat(2) = cwise_err / cwise_bnd
396  ELSE IF (cwise_err .NE. 0.0) THEN
397  tstrat(2) = 1/(16.0*eps)
398  ELSE
399  tstrat(2) = 0.0
400  END IF
401  IF (tstrat(2) .GT. 1.0) tstrat(2) = 1/(4.0*eps)
402  END IF
403  ELSE
404  cguar = 'NO'
405  IF (cwise_bnd .LT. 1.0) THEN
406  tstrat(2) = 1/(8.0*eps)
407  ELSE
408  tstrat(2) = 1.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.0)
418  $ tstrat(4) = 1.0 / tstrat(4)
419 
420  tstrat(5) = ncond / nwise_rcond
421  IF (ncond .GE. condthresh .AND. tstrat(5) .LT. 1.0)
422  $ tstrat(5) = 1.0 / tstrat(5)
423 
424  tstrat(6) = ccond / nwise_rcond
425  IF (ccond .GE. condthresh .AND. tstrat(6) .LT. 1.0)
426  $ tstrat(6) = 1.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( ' S', a2, 'SVXX: N =', i2, ', RHS = ', i2,
472  $ ', NWISE GUAR. = ', a, ', CWISE GUAR. = ', a,
473  $ ' test(',i1,') =', g12.5 )
474  9998 FORMAT( ' S', a2, 'SVXX: ', i6, ' out of ', i6,
475  $ ' tests failed to pass the threshold' )
476  9997 FORMAT( ' S', 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( ' S', a2, 'SVXX: N =', i2, ', INFO = ', i3,
492  $ ', ORCOND = ', g12.5, ', real RCOND = ', g12.5 )
493 *
494 * End of SEBCHVXX
495 *
logical function lsamen(N, CA, CB)
LSAMEN
Definition: lsamen.f:74
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
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:563
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:543
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:497
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:508
subroutine slahilb(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
SLAHILB
Definition: slahilb.f:124
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function: