1 SUBROUTINE SEBCHVXX( THRESH, PATH )
2 IMPLICIT NONE
3 * .. Scalar Arguments ..
4 REAL THRESH
5 CHARACTER*3 PATH
6 *
7 * Purpose
8 * ======
9 *
10 * SEBCHVXX will run S**SVXX on a series of Hilbert matrices and then
11 * compare the error bounds returned by SGESVXX to see if the returned
12 * answer indeed falls within those bounds.
13 *
14 * Eight test ratios will be computed. The tests will pass if they are .LT.
15 * THRESH. There are two cases that are determined by 1 / (SQRT( N ) * EPS).
16 * If that value is .LE. to the component wise reciprocal condition number,
17 * it uses the guaranteed case, other wise it uses the unguaranteed case.
18 *
19 * Test ratios:
20 * Let Xc be X_computed and Xt be X_truth.
21 * The norm used is the infinity norm.
22
23 * Let A be the guaranteed case and B be the unguaranteed case.
24 *
25 * 1. Normwise guaranteed forward error bound.
26 * A: norm ( abs( Xc - Xt ) / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ) and
27 * ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N),10) * EPS.
28 * If these conditions are met, the test ratio is set to be
29 * ERRBND( *, nwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS.
30 * B: For this case, SGESVXX should just return 1. If it is less than
31 * one, treat it the same as in 1A. Otherwise it fails. (Set test
32 * ratio to ERRBND( *, nwise_i, bnd_i ) * THRESH?)
33 *
34 * 2. Componentwise guaranteed forward error bound.
35 * A: norm ( abs( Xc(j) - Xt(j) ) ) / norm (Xt(j)) .LE. ERRBND( *, cwise_i, bnd_i )
36 * for all j .AND. ERRBND( *, cwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS.
37 * If these conditions are met, the test ratio is set to be
38 * ERRBND( *, cwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS.
39 * B: Same as normwise test ratio.
40 *
41 * 3. Backwards error.
42 * A: The test ratio is set to BERR/EPS.
43 * B: Same test ratio.
44 *
45 * 4. Reciprocal condition number.
46 * A: A condition number is computed with Xt and compared with the one
47 * returned from SGESVXX. Let RCONDc be the RCOND returned by SGESVXX
48 * and RCONDt be the RCOND from the truth value. Test ratio is set to
49 * MAX(RCONDc/RCONDt, RCONDt/RCONDc).
50 * B: Test ratio is set to 1 / (EPS * RCONDc).
51 *
52 * 5. Reciprocal normwise condition number.
53 * A: The test ratio is set to
54 * MAX(ERRBND( *, nwise_i, cond_i ) / NCOND, NCOND / ERRBND( *, nwise_i, cond_i )).
55 * B: Test ratio is set to 1 / (EPS * ERRBND( *, nwise_i, cond_i )).
56 *
57 * 7. Reciprocal componentwise condition number.
58 * A: Test ratio is set to
59 * MAX(ERRBND( *, cwise_i, cond_i ) / CCOND, CCOND / ERRBND( *, cwise_i, cond_i )).
60 * B: Test ratio is set to 1 / (EPS * ERRBND( *, cwise_i, cond_i )).
61 *
62 * .. Parameters ..
63 * NMAX is determined by the largest number in the inverse of the Hilbert
64 * matrix. Precision is exhausted when the largest entry in it is greater
65 * than 2 to the power of the number of bits in the fraction of the data
66 * type used plus one, which is 24 for single precision.
67 * NMAX should be 6 for single and 11 for double.
68
69 INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU
70 PARAMETER (NMAX = 6, NPARAMS = 2, NERRBND = 3,
71 $ NTESTS = 6)
72
73 * .. Local Scalars ..
74 INTEGER N, NRHS, INFO, I ,J, k, NFAIL, LDA, LDAB,
75 $ LDAFB, N_AUX_TESTS
76 CHARACTER FACT, TRANS, UPLO, EQUED
77 CHARACTER*2 C2
78 CHARACTER(3) NGUAR, CGUAR
79 LOGICAL printed_guide
80 REAL NCOND, CCOND, M, NORMDIF, NORMT, RCOND,
81 $ RNORM, RINORM, SUMR, SUMRI, EPS,
82 $ BERR(NMAX), RPVGRW, ORCOND,
83 $ CWISE_ERR, NWISE_ERR, CWISE_BND, NWISE_BND,
84 $ CWISE_RCOND, NWISE_RCOND,
85 $ CONDTHRESH, ERRTHRESH
86
87 * .. Local Arrays ..
88 REAL TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
89 $ A(NMAX, NMAX), ACOPY(NMAX, NMAX),
90 $ INVHILB(NMAX, NMAX), R(NMAX), C(NMAX), S(NMAX),
91 $ WORK(NMAX * 5), B(NMAX, NMAX), X(NMAX, NMAX),
92 $ DIFF(NMAX, NMAX), AF(NMAX, NMAX),
93 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
94 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
95 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX ),
96 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3)
97 INTEGER IWORK(NMAX), IPIV(NMAX)
98
99 * .. External Functions ..
100 REAL SLAMCH
101
102 * .. External Subroutines ..
103 EXTERNAL SLAHILB, SGESVXX, SSYSVXX, SPOSVXX, SGBSVXX,
104 $ SLACPY, LSAMEN
105 LOGICAL LSAMEN
106
107 * .. Intrinsic Functions ..
108 INTRINSIC SQRT, MAX, ABS
109
110 * .. Parameters ..
111 INTEGER NWISE_I, CWISE_I
112 PARAMETER (NWISE_I = 1, CWISE_I = 1)
113 INTEGER BND_I, COND_I
114 PARAMETER (BND_I = 2, COND_I = 3)
115
116 * Create the loop to test out the Hilbert matrices
117
118 FACT = 'E'
119 UPLO = 'U'
120 TRANS = 'N'
121 EQUED = 'N'
122 EPS = SLAMCH('Epsilon')
123 NFAIL = 0
124 N_AUX_TESTS = 0
125 LDA = NMAX
126 LDAB = (NMAX-1)+(NMAX-1)+1
127 LDAFB = 2*(NMAX-1)+(NMAX-1)+1
128 C2 = PATH( 2: 3 )
129
130 * Main loop to test the different Hilbert Matrices.
131
132 printed_guide = .false.
133
134 DO N = 1 , NMAX
135 PARAMS(1) = -1
136 PARAMS(2) = -1
137
138 KL = N-1
139 KU = N-1
140 NRHS = n
141 M = MAX(SQRT(REAL(N)), 10.0)
142
143 * Generate the Hilbert matrix, its inverse, and the
144 * right hand side, all scaled by the LCM(1,..,2N-1).
145 CALL SLAHILB(N, N, A, LDA, INVHILB, LDA, B, LDA, WORK, INFO)
146
147 * Copy A into ACOPY.
148 CALL SLACPY('ALL', N, N, A, NMAX, ACOPY, NMAX)
149
150 * Store A in band format for GB tests
151 DO J = 1, N
152 DO I = 1, KL+KU+1
153 AB( I, J ) = 0.0E+0
154 END DO
155 END DO
156 DO J = 1, N
157 DO I = MAX( 1, J-KU ), MIN( N, J+KL )
158 AB( KU+1+I-J, J ) = A( I, J )
159 END DO
160 END DO
161
162 * Copy AB into ABCOPY.
163 DO J = 1, N
164 DO I = 1, KL+KU+1
165 ABCOPY( I, J ) = 0.0E+0
166 END DO
167 END DO
168 CALL SLACPY('ALL', KL+KU+1, N, AB, LDAB, ABCOPY, LDAB)
169
170 * Call S**SVXX with default PARAMS and N_ERR_BND = 3.
171 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
172 CALL SSYSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
173 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
174 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
175 $ PARAMS, WORK, IWORK, INFO)
176 ELSE IF ( LSAMEN( 2, C2, 'PO' ) ) THEN
177 CALL SPOSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
178 $ EQUED, S, B, LDA, X, LDA, ORCOND,
179 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
180 $ PARAMS, WORK, IWORK, INFO)
181 ELSE IF ( LSAMEN( 2, C2, 'GB' ) ) THEN
182 CALL SGBSVXX(FACT, TRANS, N, KL, KU, NRHS, ABCOPY,
183 $ LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B,
184 $ LDA, X, LDA, ORCOND, RPVGRW, BERR, NERRBND,
185 $ ERRBND_N, ERRBND_C, NPARAMS, PARAMS, WORK, IWORK,
186 $ INFO)
187 ELSE
188 CALL SGESVXX(FACT, TRANS, N, NRHS, ACOPY, LDA, AF, LDA,
189 $ IPIV, EQUED, R, C, B, LDA, X, LDA, ORCOND,
190 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
191 $ PARAMS, WORK, IWORK, INFO)
192 END IF
193
194 N_AUX_TESTS = N_AUX_TESTS + 1
195 IF (ORCOND .LT. EPS) THEN
196 ! Either factorization failed or the matrix is flagged, and 1 <=
197 ! INFO <= N+1. We don't decide based on rcond anymore.
198 ! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN
199 ! NFAIL = NFAIL + 1
200 ! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND
201 ! END IF
202 ELSE
203 ! Either everything succeeded (INFO == 0) or some solution failed
204 ! to converge (INFO > N+1).
205 IF (INFO .GT. 0 .AND. INFO .LE. N+1) THEN
206 NFAIL = NFAIL + 1
207 WRITE (*, FMT=8000) C2, N, INFO, ORCOND, RCOND
208 END IF
209 END IF
210
211 * Calculating the difference between S**SVXX's X and the true X.
212 DO I = 1, N
213 DO J = 1, NRHS
214 DIFF( I, J ) = X( I, J ) - INVHILB( I, J )
215 END DO
216 END DO
217
218 * Calculating the RCOND
219 RNORM = 0
220 RINORM = 0
221 IF ( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'SY' ) ) THEN
222 DO I = 1, N
223 SUMR = 0
224 SUMRI = 0
225 DO J = 1, N
226 SUMR = SUMR + ABS(S(I) * A(I,J) * S(J))
227 SUMRI = SUMRI + ABS(INVHILB(I, J) / S(J) / S(I))
228 END DO
229 RNORM = MAX(RNORM,SUMR)
230 RINORM = MAX(RINORM,SUMRI)
231 END DO
232 ELSE IF ( LSAMEN( 2, C2, 'GE' ) .OR. LSAMEN( 2, C2, 'GB' ) )
233 $ THEN
234 DO I = 1, N
235 SUMR = 0
236 SUMRI = 0
237 DO J = 1, N
238 SUMR = SUMR + ABS(R(I) * A(I,J) * C(J))
239 SUMRI = SUMRI + ABS(INVHILB(I, J) / R(J) / C(I))
240 END DO
241 RNORM = MAX(RNORM,SUMR)
242 RINORM = MAX(RINORM,SUMRI)
243 END DO
244 END IF
245
246 RNORM = RNORM / A(1, 1)
247 RCOND = 1.0/(RNORM * RINORM)
248
249 * Calculating the R for normwise rcond.
250 DO I = 1, N
251 RINV(I) = 0.0
252 END DO
253 DO J = 1, N
254 DO I = 1, N
255 RINV(I) = RINV(I) + ABS(A(I,J))
256 END DO
257 END DO
258
259 * Calculating the Normwise rcond.
260 RINORM = 0.0
261 DO I = 1, N
262 SUMRI = 0.0
263 DO J = 1, N
264 SUMRI = SUMRI + ABS(INVHILB(I,J) * RINV(J))
265 END DO
266 RINORM = MAX(RINORM, SUMRI)
267 END DO
268
269 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
270 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
271 NCOND = A(1,1) / RINORM
272
273 CONDTHRESH = M * EPS
274 ERRTHRESH = M * EPS
275
276 DO K = 1, NRHS
277 NORMT = 0.0
278 NORMDIF = 0.0
279 CWISE_ERR = 0.0
280 DO I = 1, N
281 NORMT = MAX(ABS(INVHILB(I, K)), NORMT)
282 NORMDIF = MAX(ABS(X(I,K) - INVHILB(I,K)), NORMDIF)
283 IF (INVHILB(I,K) .NE. 0.0) THEN
284 CWISE_ERR = MAX(ABS(X(I,K) - INVHILB(I,K))
285 $ /ABS(INVHILB(I,K)), CWISE_ERR)
286 ELSE IF (X(I, K) .NE. 0.0) THEN
287 CWISE_ERR = SLAMCH('OVERFLOW')
288 END IF
289 END DO
290 IF (NORMT .NE. 0.0) THEN
291 NWISE_ERR = NORMDIF / NORMT
292 ELSE IF (NORMDIF .NE. 0.0) THEN
293 NWISE_ERR = SLAMCH('OVERFLOW')
294 ELSE
295 NWISE_ERR = 0.0
296 ENDIF
297
298 DO I = 1, N
299 RINV(I) = 0.0
300 END DO
301 DO J = 1, N
302 DO I = 1, N
303 RINV(I) = RINV(I) + ABS(A(I, J) * INVHILB(J, K))
304 END DO
305 END DO
306 RINORM = 0.0
307 DO I = 1, N
308 SUMRI = 0.0
309 DO J = 1, N
310 SUMRI = SUMRI
311 $ + ABS(INVHILB(I, J) * RINV(J) / INVHILB(I, K))
312 END DO
313 RINORM = MAX(RINORM, SUMRI)
314 END DO
315 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
316 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
317 CCOND = A(1,1)/RINORM
318
319 ! Forward error bound tests
320 NWISE_BND = ERRBND_N(K + (BND_I-1)*NRHS)
321 CWISE_BND = ERRBND_C(K + (BND_I-1)*NRHS)
322 NWISE_RCOND = ERRBND_N(K + (COND_I-1)*NRHS)
323 CWISE_RCOND = ERRBND_C(K + (COND_I-1)*NRHS)
324 ! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
325 ! $ condthresh, ncond.ge.condthresh
326 ! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
327
328 IF (NCOND .GE. CONDTHRESH) THEN
329 NGUAR = 'YES'
330 IF (NWISE_BND .GT. ERRTHRESH) THEN
331 TSTRAT(1) = 1/(2.0*EPS)
332 ELSE
333
334 IF (NWISE_BND .NE. 0.0) THEN
335 TSTRAT(1) = NWISE_ERR / NWISE_BND
336 ELSE IF (NWISE_ERR .NE. 0.0) THEN
337 TSTRAT(1) = 1/(16.0*EPS)
338 ELSE
339 TSTRAT(1) = 0.0
340 END IF
341 IF (TSTRAT(1) .GT. 1.0) THEN
342 TSTRAT(1) = 1/(4.0*EPS)
343 END IF
344 END IF
345 ELSE
346 NGUAR = 'NO'
347 IF (NWISE_BND .LT. 1.0) THEN
348 TSTRAT(1) = 1/(8.0*EPS)
349 ELSE
350 TSTRAT(1) = 1.0
351 END IF
352 END IF
353 ! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
354 ! $ condthresh, ccond.ge.condthresh
355 ! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
356 IF (CCOND .GE. CONDTHRESH) THEN
357 CGUAR = 'YES'
358
359 IF (CWISE_BND .GT. ERRTHRESH) THEN
360 TSTRAT(2) = 1/(2.0*EPS)
361 ELSE
362 IF (CWISE_BND .NE. 0.0) THEN
363 TSTRAT(2) = CWISE_ERR / CWISE_BND
364 ELSE IF (CWISE_ERR .NE. 0.0) THEN
365 TSTRAT(2) = 1/(16.0*EPS)
366 ELSE
367 TSTRAT(2) = 0.0
368 END IF
369 IF (TSTRAT(2) .GT. 1.0) TSTRAT(2) = 1/(4.0*EPS)
370 END IF
371 ELSE
372 CGUAR = 'NO'
373 IF (CWISE_BND .LT. 1.0) THEN
374 TSTRAT(2) = 1/(8.0*EPS)
375 ELSE
376 TSTRAT(2) = 1.0
377 END IF
378 END IF
379
380 ! Backwards error test
381 TSTRAT(3) = BERR(K)/EPS
382
383 ! Condition number tests
384 TSTRAT(4) = RCOND / ORCOND
385 IF (RCOND .GE. CONDTHRESH .AND. TSTRAT(4) .LT. 1.0)
386 $ TSTRAT(4) = 1.0 / TSTRAT(4)
387
388 TSTRAT(5) = NCOND / NWISE_RCOND
389 IF (NCOND .GE. CONDTHRESH .AND. TSTRAT(5) .LT. 1.0)
390 $ TSTRAT(5) = 1.0 / TSTRAT(5)
391
392 TSTRAT(6) = CCOND / NWISE_RCOND
393 IF (CCOND .GE. CONDTHRESH .AND. TSTRAT(6) .LT. 1.0)
394 $ TSTRAT(6) = 1.0 / TSTRAT(6)
395
396 DO I = 1, NTESTS
397 IF (TSTRAT(I) .GT. THRESH) THEN
398 IF (.NOT.PRINTED_GUIDE) THEN
399 WRITE(*,*)
400 WRITE( *, 9996) 1
401 WRITE( *, 9995) 2
402 WRITE( *, 9994) 3
403 WRITE( *, 9993) 4
404 WRITE( *, 9992) 5
405 WRITE( *, 9991) 6
406 WRITE( *, 9990) 7
407 WRITE( *, 9989) 8
408 WRITE(*,*)
409 PRINTED_GUIDE = .TRUE.
410 END IF
411 WRITE( *, 9999) C2, N, K, NGUAR, CGUAR, I, TSTRAT(I)
412 NFAIL = NFAIL + 1
413 END IF
414 END DO
415 END DO
416
417 c$$$ WRITE(*,*)
418 c$$$ WRITE(*,*) 'Normwise Error Bounds'
419 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i)
420 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i)
421 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i)
422 c$$$ WRITE(*,*)
423 c$$$ WRITE(*,*) 'Componentwise Error Bounds'
424 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i)
425 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i)
426 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i)
427 c$$$ print *, 'Info: ', info
428 c$$$ WRITE(*,*)
429 * WRITE(*,*) 'TSTRAT: ',TSTRAT
430
431 END DO
432
433 WRITE(*,*)
434 IF( NFAIL .GT. 0 ) THEN
435 WRITE(*,9998) C2, NFAIL, NTESTS*N+N_AUX_TESTS
436 ELSE
437 WRITE(*,9997) C2
438 END IF
439 9999 FORMAT( ' S', A2, 'SVXX: N =', I2, ', RHS = ', I2,
440 $ ', NWISE GUAR. = ', A, ', CWISE GUAR. = ', A,
441 $ ' test(',I1,') =', G12.5 )
442 9998 FORMAT( ' S', A2, 'SVXX: ', I6, ' out of ', I6,
443 $ ' tests failed to pass the threshold' )
444 9997 FORMAT( ' S', A2, 'SVXX passed the tests of error bounds' )
445 * Test ratios.
446 9996 FORMAT( 3X, I2, ': Normwise guaranteed forward error', / 5X,
447 $ 'Guaranteed case: if norm ( abs( Xc - Xt )',
448 $ ' / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then',
449 $ / 5X,
450 $ 'ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS')
451 9995 FORMAT( 3X, I2, ': Componentwise guaranteed forward error' )
452 9994 FORMAT( 3X, I2, ': Backwards error' )
453 9993 FORMAT( 3X, I2, ': Reciprocal condition number' )
454 9992 FORMAT( 3X, I2, ': Reciprocal normwise condition number' )
455 9991 FORMAT( 3X, I2, ': Raw normwise error estimate' )
456 9990 FORMAT( 3X, I2, ': Reciprocal componentwise condition number' )
457 9989 FORMAT( 3X, I2, ': Raw componentwise error estimate' )
458
459 8000 FORMAT( ' S', A2, 'SVXX: N =', I2, ', INFO = ', I3,
460 $ ', ORCOND = ', G12.5, ', real RCOND = ', G12.5 )
461
462 END
2 IMPLICIT NONE
3 * .. Scalar Arguments ..
4 REAL THRESH
5 CHARACTER*3 PATH
6 *
7 * Purpose
8 * ======
9 *
10 * SEBCHVXX will run S**SVXX on a series of Hilbert matrices and then
11 * compare the error bounds returned by SGESVXX to see if the returned
12 * answer indeed falls within those bounds.
13 *
14 * Eight test ratios will be computed. The tests will pass if they are .LT.
15 * THRESH. There are two cases that are determined by 1 / (SQRT( N ) * EPS).
16 * If that value is .LE. to the component wise reciprocal condition number,
17 * it uses the guaranteed case, other wise it uses the unguaranteed case.
18 *
19 * Test ratios:
20 * Let Xc be X_computed and Xt be X_truth.
21 * The norm used is the infinity norm.
22
23 * Let A be the guaranteed case and B be the unguaranteed case.
24 *
25 * 1. Normwise guaranteed forward error bound.
26 * A: norm ( abs( Xc - Xt ) / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ) and
27 * ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N),10) * EPS.
28 * If these conditions are met, the test ratio is set to be
29 * ERRBND( *, nwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS.
30 * B: For this case, SGESVXX should just return 1. If it is less than
31 * one, treat it the same as in 1A. Otherwise it fails. (Set test
32 * ratio to ERRBND( *, nwise_i, bnd_i ) * THRESH?)
33 *
34 * 2. Componentwise guaranteed forward error bound.
35 * A: norm ( abs( Xc(j) - Xt(j) ) ) / norm (Xt(j)) .LE. ERRBND( *, cwise_i, bnd_i )
36 * for all j .AND. ERRBND( *, cwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS.
37 * If these conditions are met, the test ratio is set to be
38 * ERRBND( *, cwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS.
39 * B: Same as normwise test ratio.
40 *
41 * 3. Backwards error.
42 * A: The test ratio is set to BERR/EPS.
43 * B: Same test ratio.
44 *
45 * 4. Reciprocal condition number.
46 * A: A condition number is computed with Xt and compared with the one
47 * returned from SGESVXX. Let RCONDc be the RCOND returned by SGESVXX
48 * and RCONDt be the RCOND from the truth value. Test ratio is set to
49 * MAX(RCONDc/RCONDt, RCONDt/RCONDc).
50 * B: Test ratio is set to 1 / (EPS * RCONDc).
51 *
52 * 5. Reciprocal normwise condition number.
53 * A: The test ratio is set to
54 * MAX(ERRBND( *, nwise_i, cond_i ) / NCOND, NCOND / ERRBND( *, nwise_i, cond_i )).
55 * B: Test ratio is set to 1 / (EPS * ERRBND( *, nwise_i, cond_i )).
56 *
57 * 7. Reciprocal componentwise condition number.
58 * A: Test ratio is set to
59 * MAX(ERRBND( *, cwise_i, cond_i ) / CCOND, CCOND / ERRBND( *, cwise_i, cond_i )).
60 * B: Test ratio is set to 1 / (EPS * ERRBND( *, cwise_i, cond_i )).
61 *
62 * .. Parameters ..
63 * NMAX is determined by the largest number in the inverse of the Hilbert
64 * matrix. Precision is exhausted when the largest entry in it is greater
65 * than 2 to the power of the number of bits in the fraction of the data
66 * type used plus one, which is 24 for single precision.
67 * NMAX should be 6 for single and 11 for double.
68
69 INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU
70 PARAMETER (NMAX = 6, NPARAMS = 2, NERRBND = 3,
71 $ NTESTS = 6)
72
73 * .. Local Scalars ..
74 INTEGER N, NRHS, INFO, I ,J, k, NFAIL, LDA, LDAB,
75 $ LDAFB, N_AUX_TESTS
76 CHARACTER FACT, TRANS, UPLO, EQUED
77 CHARACTER*2 C2
78 CHARACTER(3) NGUAR, CGUAR
79 LOGICAL printed_guide
80 REAL NCOND, CCOND, M, NORMDIF, NORMT, RCOND,
81 $ RNORM, RINORM, SUMR, SUMRI, EPS,
82 $ BERR(NMAX), RPVGRW, ORCOND,
83 $ CWISE_ERR, NWISE_ERR, CWISE_BND, NWISE_BND,
84 $ CWISE_RCOND, NWISE_RCOND,
85 $ CONDTHRESH, ERRTHRESH
86
87 * .. Local Arrays ..
88 REAL TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
89 $ A(NMAX, NMAX), ACOPY(NMAX, NMAX),
90 $ INVHILB(NMAX, NMAX), R(NMAX), C(NMAX), S(NMAX),
91 $ WORK(NMAX * 5), B(NMAX, NMAX), X(NMAX, NMAX),
92 $ DIFF(NMAX, NMAX), AF(NMAX, NMAX),
93 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
94 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
95 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX ),
96 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3)
97 INTEGER IWORK(NMAX), IPIV(NMAX)
98
99 * .. External Functions ..
100 REAL SLAMCH
101
102 * .. External Subroutines ..
103 EXTERNAL SLAHILB, SGESVXX, SSYSVXX, SPOSVXX, SGBSVXX,
104 $ SLACPY, LSAMEN
105 LOGICAL LSAMEN
106
107 * .. Intrinsic Functions ..
108 INTRINSIC SQRT, MAX, ABS
109
110 * .. Parameters ..
111 INTEGER NWISE_I, CWISE_I
112 PARAMETER (NWISE_I = 1, CWISE_I = 1)
113 INTEGER BND_I, COND_I
114 PARAMETER (BND_I = 2, COND_I = 3)
115
116 * Create the loop to test out the Hilbert matrices
117
118 FACT = 'E'
119 UPLO = 'U'
120 TRANS = 'N'
121 EQUED = 'N'
122 EPS = SLAMCH('Epsilon')
123 NFAIL = 0
124 N_AUX_TESTS = 0
125 LDA = NMAX
126 LDAB = (NMAX-1)+(NMAX-1)+1
127 LDAFB = 2*(NMAX-1)+(NMAX-1)+1
128 C2 = PATH( 2: 3 )
129
130 * Main loop to test the different Hilbert Matrices.
131
132 printed_guide = .false.
133
134 DO N = 1 , NMAX
135 PARAMS(1) = -1
136 PARAMS(2) = -1
137
138 KL = N-1
139 KU = N-1
140 NRHS = n
141 M = MAX(SQRT(REAL(N)), 10.0)
142
143 * Generate the Hilbert matrix, its inverse, and the
144 * right hand side, all scaled by the LCM(1,..,2N-1).
145 CALL SLAHILB(N, N, A, LDA, INVHILB, LDA, B, LDA, WORK, INFO)
146
147 * Copy A into ACOPY.
148 CALL SLACPY('ALL', N, N, A, NMAX, ACOPY, NMAX)
149
150 * Store A in band format for GB tests
151 DO J = 1, N
152 DO I = 1, KL+KU+1
153 AB( I, J ) = 0.0E+0
154 END DO
155 END DO
156 DO J = 1, N
157 DO I = MAX( 1, J-KU ), MIN( N, J+KL )
158 AB( KU+1+I-J, J ) = A( I, J )
159 END DO
160 END DO
161
162 * Copy AB into ABCOPY.
163 DO J = 1, N
164 DO I = 1, KL+KU+1
165 ABCOPY( I, J ) = 0.0E+0
166 END DO
167 END DO
168 CALL SLACPY('ALL', KL+KU+1, N, AB, LDAB, ABCOPY, LDAB)
169
170 * Call S**SVXX with default PARAMS and N_ERR_BND = 3.
171 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
172 CALL SSYSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
173 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
174 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
175 $ PARAMS, WORK, IWORK, INFO)
176 ELSE IF ( LSAMEN( 2, C2, 'PO' ) ) THEN
177 CALL SPOSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
178 $ EQUED, S, B, LDA, X, LDA, ORCOND,
179 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
180 $ PARAMS, WORK, IWORK, INFO)
181 ELSE IF ( LSAMEN( 2, C2, 'GB' ) ) THEN
182 CALL SGBSVXX(FACT, TRANS, N, KL, KU, NRHS, ABCOPY,
183 $ LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B,
184 $ LDA, X, LDA, ORCOND, RPVGRW, BERR, NERRBND,
185 $ ERRBND_N, ERRBND_C, NPARAMS, PARAMS, WORK, IWORK,
186 $ INFO)
187 ELSE
188 CALL SGESVXX(FACT, TRANS, N, NRHS, ACOPY, LDA, AF, LDA,
189 $ IPIV, EQUED, R, C, B, LDA, X, LDA, ORCOND,
190 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
191 $ PARAMS, WORK, IWORK, INFO)
192 END IF
193
194 N_AUX_TESTS = N_AUX_TESTS + 1
195 IF (ORCOND .LT. EPS) THEN
196 ! Either factorization failed or the matrix is flagged, and 1 <=
197 ! INFO <= N+1. We don't decide based on rcond anymore.
198 ! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN
199 ! NFAIL = NFAIL + 1
200 ! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND
201 ! END IF
202 ELSE
203 ! Either everything succeeded (INFO == 0) or some solution failed
204 ! to converge (INFO > N+1).
205 IF (INFO .GT. 0 .AND. INFO .LE. N+1) THEN
206 NFAIL = NFAIL + 1
207 WRITE (*, FMT=8000) C2, N, INFO, ORCOND, RCOND
208 END IF
209 END IF
210
211 * Calculating the difference between S**SVXX's X and the true X.
212 DO I = 1, N
213 DO J = 1, NRHS
214 DIFF( I, J ) = X( I, J ) - INVHILB( I, J )
215 END DO
216 END DO
217
218 * Calculating the RCOND
219 RNORM = 0
220 RINORM = 0
221 IF ( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'SY' ) ) THEN
222 DO I = 1, N
223 SUMR = 0
224 SUMRI = 0
225 DO J = 1, N
226 SUMR = SUMR + ABS(S(I) * A(I,J) * S(J))
227 SUMRI = SUMRI + ABS(INVHILB(I, J) / S(J) / S(I))
228 END DO
229 RNORM = MAX(RNORM,SUMR)
230 RINORM = MAX(RINORM,SUMRI)
231 END DO
232 ELSE IF ( LSAMEN( 2, C2, 'GE' ) .OR. LSAMEN( 2, C2, 'GB' ) )
233 $ THEN
234 DO I = 1, N
235 SUMR = 0
236 SUMRI = 0
237 DO J = 1, N
238 SUMR = SUMR + ABS(R(I) * A(I,J) * C(J))
239 SUMRI = SUMRI + ABS(INVHILB(I, J) / R(J) / C(I))
240 END DO
241 RNORM = MAX(RNORM,SUMR)
242 RINORM = MAX(RINORM,SUMRI)
243 END DO
244 END IF
245
246 RNORM = RNORM / A(1, 1)
247 RCOND = 1.0/(RNORM * RINORM)
248
249 * Calculating the R for normwise rcond.
250 DO I = 1, N
251 RINV(I) = 0.0
252 END DO
253 DO J = 1, N
254 DO I = 1, N
255 RINV(I) = RINV(I) + ABS(A(I,J))
256 END DO
257 END DO
258
259 * Calculating the Normwise rcond.
260 RINORM = 0.0
261 DO I = 1, N
262 SUMRI = 0.0
263 DO J = 1, N
264 SUMRI = SUMRI + ABS(INVHILB(I,J) * RINV(J))
265 END DO
266 RINORM = MAX(RINORM, SUMRI)
267 END DO
268
269 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
270 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
271 NCOND = A(1,1) / RINORM
272
273 CONDTHRESH = M * EPS
274 ERRTHRESH = M * EPS
275
276 DO K = 1, NRHS
277 NORMT = 0.0
278 NORMDIF = 0.0
279 CWISE_ERR = 0.0
280 DO I = 1, N
281 NORMT = MAX(ABS(INVHILB(I, K)), NORMT)
282 NORMDIF = MAX(ABS(X(I,K) - INVHILB(I,K)), NORMDIF)
283 IF (INVHILB(I,K) .NE. 0.0) THEN
284 CWISE_ERR = MAX(ABS(X(I,K) - INVHILB(I,K))
285 $ /ABS(INVHILB(I,K)), CWISE_ERR)
286 ELSE IF (X(I, K) .NE. 0.0) THEN
287 CWISE_ERR = SLAMCH('OVERFLOW')
288 END IF
289 END DO
290 IF (NORMT .NE. 0.0) THEN
291 NWISE_ERR = NORMDIF / NORMT
292 ELSE IF (NORMDIF .NE. 0.0) THEN
293 NWISE_ERR = SLAMCH('OVERFLOW')
294 ELSE
295 NWISE_ERR = 0.0
296 ENDIF
297
298 DO I = 1, N
299 RINV(I) = 0.0
300 END DO
301 DO J = 1, N
302 DO I = 1, N
303 RINV(I) = RINV(I) + ABS(A(I, J) * INVHILB(J, K))
304 END DO
305 END DO
306 RINORM = 0.0
307 DO I = 1, N
308 SUMRI = 0.0
309 DO J = 1, N
310 SUMRI = SUMRI
311 $ + ABS(INVHILB(I, J) * RINV(J) / INVHILB(I, K))
312 END DO
313 RINORM = MAX(RINORM, SUMRI)
314 END DO
315 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
316 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
317 CCOND = A(1,1)/RINORM
318
319 ! Forward error bound tests
320 NWISE_BND = ERRBND_N(K + (BND_I-1)*NRHS)
321 CWISE_BND = ERRBND_C(K + (BND_I-1)*NRHS)
322 NWISE_RCOND = ERRBND_N(K + (COND_I-1)*NRHS)
323 CWISE_RCOND = ERRBND_C(K + (COND_I-1)*NRHS)
324 ! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
325 ! $ condthresh, ncond.ge.condthresh
326 ! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
327
328 IF (NCOND .GE. CONDTHRESH) THEN
329 NGUAR = 'YES'
330 IF (NWISE_BND .GT. ERRTHRESH) THEN
331 TSTRAT(1) = 1/(2.0*EPS)
332 ELSE
333
334 IF (NWISE_BND .NE. 0.0) THEN
335 TSTRAT(1) = NWISE_ERR / NWISE_BND
336 ELSE IF (NWISE_ERR .NE. 0.0) THEN
337 TSTRAT(1) = 1/(16.0*EPS)
338 ELSE
339 TSTRAT(1) = 0.0
340 END IF
341 IF (TSTRAT(1) .GT. 1.0) THEN
342 TSTRAT(1) = 1/(4.0*EPS)
343 END IF
344 END IF
345 ELSE
346 NGUAR = 'NO'
347 IF (NWISE_BND .LT. 1.0) THEN
348 TSTRAT(1) = 1/(8.0*EPS)
349 ELSE
350 TSTRAT(1) = 1.0
351 END IF
352 END IF
353 ! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
354 ! $ condthresh, ccond.ge.condthresh
355 ! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
356 IF (CCOND .GE. CONDTHRESH) THEN
357 CGUAR = 'YES'
358
359 IF (CWISE_BND .GT. ERRTHRESH) THEN
360 TSTRAT(2) = 1/(2.0*EPS)
361 ELSE
362 IF (CWISE_BND .NE. 0.0) THEN
363 TSTRAT(2) = CWISE_ERR / CWISE_BND
364 ELSE IF (CWISE_ERR .NE. 0.0) THEN
365 TSTRAT(2) = 1/(16.0*EPS)
366 ELSE
367 TSTRAT(2) = 0.0
368 END IF
369 IF (TSTRAT(2) .GT. 1.0) TSTRAT(2) = 1/(4.0*EPS)
370 END IF
371 ELSE
372 CGUAR = 'NO'
373 IF (CWISE_BND .LT. 1.0) THEN
374 TSTRAT(2) = 1/(8.0*EPS)
375 ELSE
376 TSTRAT(2) = 1.0
377 END IF
378 END IF
379
380 ! Backwards error test
381 TSTRAT(3) = BERR(K)/EPS
382
383 ! Condition number tests
384 TSTRAT(4) = RCOND / ORCOND
385 IF (RCOND .GE. CONDTHRESH .AND. TSTRAT(4) .LT. 1.0)
386 $ TSTRAT(4) = 1.0 / TSTRAT(4)
387
388 TSTRAT(5) = NCOND / NWISE_RCOND
389 IF (NCOND .GE. CONDTHRESH .AND. TSTRAT(5) .LT. 1.0)
390 $ TSTRAT(5) = 1.0 / TSTRAT(5)
391
392 TSTRAT(6) = CCOND / NWISE_RCOND
393 IF (CCOND .GE. CONDTHRESH .AND. TSTRAT(6) .LT. 1.0)
394 $ TSTRAT(6) = 1.0 / TSTRAT(6)
395
396 DO I = 1, NTESTS
397 IF (TSTRAT(I) .GT. THRESH) THEN
398 IF (.NOT.PRINTED_GUIDE) THEN
399 WRITE(*,*)
400 WRITE( *, 9996) 1
401 WRITE( *, 9995) 2
402 WRITE( *, 9994) 3
403 WRITE( *, 9993) 4
404 WRITE( *, 9992) 5
405 WRITE( *, 9991) 6
406 WRITE( *, 9990) 7
407 WRITE( *, 9989) 8
408 WRITE(*,*)
409 PRINTED_GUIDE = .TRUE.
410 END IF
411 WRITE( *, 9999) C2, N, K, NGUAR, CGUAR, I, TSTRAT(I)
412 NFAIL = NFAIL + 1
413 END IF
414 END DO
415 END DO
416
417 c$$$ WRITE(*,*)
418 c$$$ WRITE(*,*) 'Normwise Error Bounds'
419 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i)
420 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i)
421 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i)
422 c$$$ WRITE(*,*)
423 c$$$ WRITE(*,*) 'Componentwise Error Bounds'
424 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i)
425 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i)
426 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i)
427 c$$$ print *, 'Info: ', info
428 c$$$ WRITE(*,*)
429 * WRITE(*,*) 'TSTRAT: ',TSTRAT
430
431 END DO
432
433 WRITE(*,*)
434 IF( NFAIL .GT. 0 ) THEN
435 WRITE(*,9998) C2, NFAIL, NTESTS*N+N_AUX_TESTS
436 ELSE
437 WRITE(*,9997) C2
438 END IF
439 9999 FORMAT( ' S', A2, 'SVXX: N =', I2, ', RHS = ', I2,
440 $ ', NWISE GUAR. = ', A, ', CWISE GUAR. = ', A,
441 $ ' test(',I1,') =', G12.5 )
442 9998 FORMAT( ' S', A2, 'SVXX: ', I6, ' out of ', I6,
443 $ ' tests failed to pass the threshold' )
444 9997 FORMAT( ' S', A2, 'SVXX passed the tests of error bounds' )
445 * Test ratios.
446 9996 FORMAT( 3X, I2, ': Normwise guaranteed forward error', / 5X,
447 $ 'Guaranteed case: if norm ( abs( Xc - Xt )',
448 $ ' / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then',
449 $ / 5X,
450 $ 'ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS')
451 9995 FORMAT( 3X, I2, ': Componentwise guaranteed forward error' )
452 9994 FORMAT( 3X, I2, ': Backwards error' )
453 9993 FORMAT( 3X, I2, ': Reciprocal condition number' )
454 9992 FORMAT( 3X, I2, ': Reciprocal normwise condition number' )
455 9991 FORMAT( 3X, I2, ': Raw normwise error estimate' )
456 9990 FORMAT( 3X, I2, ': Reciprocal componentwise condition number' )
457 9989 FORMAT( 3X, I2, ': Raw componentwise error estimate' )
458
459 8000 FORMAT( ' S', A2, 'SVXX: N =', I2, ', INFO = ', I3,
460 $ ', ORCOND = ', G12.5, ', real RCOND = ', G12.5 )
461
462 END