1 SUBROUTINE CEBCHVXX( THRESH, PATH )
2 IMPLICIT NONE
3 * .. Scalar Arguments ..
4 REAL THRESH
5 CHARACTER*3 PATH
6 *
7 * Purpose
8 * ======
9 *
10 * CEBCHVXX will run CGESVXX on a series of Hilbert matrices and then
11 * compare the error bounds returned by CGESVXX 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, CGESVXX 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 CGESVXX. Let RCONDc be the RCOND returned by CGESVXX
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 * 6. 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,
75 $ N_AUX_TESTS, LDAB, LDAFB
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 COMPLEX ZDUM
87
88 * .. Local Arrays ..
89 REAL TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
90 $ S(NMAX), R(NMAX),C(NMAX),RWORK(3*NMAX),
91 $ DIFF(NMAX, NMAX),
92 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3)
93 INTEGER IPIV(NMAX)
94 COMPLEX A(NMAX,NMAX),INVHILB(NMAX,NMAX),X(NMAX,NMAX),
95 $ WORK(NMAX*3*5), AF(NMAX, NMAX),B(NMAX, NMAX),
96 $ ACOPY(NMAX, NMAX),
97 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
98 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
99 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX )
100
101 * .. External Functions ..
102 REAL SLAMCH
103
104 * .. External Subroutines ..
105 EXTERNAL CLAHILB, CGESVXX, CSYSVXX, CPOSVXX,
106 $ CGBSVXX, CLACPY, LSAMEN
107 LOGICAL LSAMEN
108
109 * .. Intrinsic Functions ..
110 INTRINSIC SQRT, MAX, ABS, REAL, AIMAG
111
112 * .. Statement Functions ..
113 REAL CABS1
114 * ..
115 * .. Statement Function Definitions ..
116 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
117
118 * .. Parameters ..
119 INTEGER NWISE_I, CWISE_I
120 PARAMETER (NWISE_I = 1, CWISE_I = 1)
121 INTEGER BND_I, COND_I
122 PARAMETER (BND_I = 2, COND_I = 3)
123
124 * Create the loop to test out the Hilbert matrices
125
126 FACT = 'E'
127 UPLO = 'U'
128 TRANS = 'N'
129 EQUED = 'N'
130 EPS = SLAMCH('Epsilon')
131 NFAIL = 0
132 N_AUX_TESTS = 0
133 LDA = NMAX
134 LDAB = (NMAX-1)+(NMAX-1)+1
135 LDAFB = 2*(NMAX-1)+(NMAX-1)+1
136 C2 = PATH( 2: 3 )
137
138 * Main loop to test the different Hilbert Matrices.
139
140 printed_guide = .false.
141
142 DO N = 1 , NMAX
143 PARAMS(1) = -1
144 PARAMS(2) = -1
145
146 KL = N-1
147 KU = N-1
148 NRHS = n
149 M = MAX(SQRT(REAL(N)), 10.0)
150
151 * Generate the Hilbert matrix, its inverse, and the
152 * right hand side, all scaled by the LCM(1,..,2N-1).
153 CALL CLAHILB(N, N, A, LDA, INVHILB, LDA, B,
154 $ LDA, WORK, INFO, PATH)
155
156 * Copy A into ACOPY.
157 CALL CLACPY('ALL', N, N, A, NMAX, ACOPY, NMAX)
158
159 * Store A in band format for GB tests
160 DO J = 1, N
161 DO I = 1, KL+KU+1
162 AB( I, J ) = (0.0E+0,0.0E+0)
163 END DO
164 END DO
165 DO J = 1, N
166 DO I = MAX( 1, J-KU ), MIN( N, J+KL )
167 AB( KU+1+I-J, J ) = A( I, J )
168 END DO
169 END DO
170
171 * Copy AB into ABCOPY.
172 DO J = 1, N
173 DO I = 1, KL+KU+1
174 ABCOPY( I, J ) = (0.0E+0,0.0E+0)
175 END DO
176 END DO
177 CALL CLACPY('ALL', KL+KU+1, N, AB, LDAB, ABCOPY, LDAB)
178
179 * Call C**SVXX with default PARAMS and N_ERR_BND = 3.
180 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
181 CALL CSYSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
182 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
183 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
184 $ PARAMS, WORK, RWORK, INFO)
185 ELSE IF ( LSAMEN( 2, C2, 'PO' ) ) THEN
186 CALL CPOSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
187 $ EQUED, S, B, LDA, X, LDA, ORCOND,
188 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
189 $ PARAMS, WORK, RWORK, INFO)
190 ELSE IF ( LSAMEN( 2, C2, 'HE' ) ) THEN
191 CALL CHESVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
192 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
193 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
194 $ PARAMS, WORK, RWORK, INFO)
195 ELSE IF ( LSAMEN( 2, C2, 'GB' ) ) THEN
196 CALL CGBSVXX(FACT, TRANS, N, KL, KU, NRHS, ABCOPY,
197 $ LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B,
198 $ LDA, X, LDA, ORCOND, RPVGRW, BERR, NERRBND,
199 $ ERRBND_N, ERRBND_C, NPARAMS, PARAMS, WORK, RWORK,
200 $ INFO)
201 ELSE
202 CALL CGESVXX(FACT, TRANS, N, NRHS, ACOPY, LDA, AF, LDA,
203 $ IPIV, EQUED, R, C, B, LDA, X, LDA, ORCOND,
204 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
205 $ PARAMS, WORK, RWORK, INFO)
206 END IF
207
208 N_AUX_TESTS = N_AUX_TESTS + 1
209 IF (ORCOND .LT. EPS) THEN
210 ! Either factorization failed or the matrix is flagged, and 1 <=
211 ! INFO <= N+1. We don't decide based on rcond anymore.
212 ! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN
213 ! NFAIL = NFAIL + 1
214 ! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND
215 ! END IF
216 ELSE
217 ! Either everything succeeded (INFO == 0) or some solution failed
218 ! to converge (INFO > N+1).
219 IF (INFO .GT. 0 .AND. INFO .LE. N+1) THEN
220 NFAIL = NFAIL + 1
221 WRITE (*, FMT=8000) C2, N, INFO, ORCOND, RCOND
222 END IF
223 END IF
224
225 * Calculating the difference between C**SVXX's X and the true X.
226 DO I = 1,N
227 DO J =1,NRHS
228 DIFF(I,J) = X(I,J) - INVHILB(I,J)
229 END DO
230 END DO
231
232 * Calculating the RCOND
233 RNORM = 0
234 RINORM = 0
235 IF ( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'SY' ) .OR.
236 $ LSAMEN( 2, C2, 'HE' ) ) THEN
237 DO I = 1, N
238 SUMR = 0
239 SUMRI = 0
240 DO J = 1, N
241 SUMR = SUMR + S(I) * CABS1(A(I,J)) * S(J)
242 SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (S(J) * S(I))
243 END DO
244 RNORM = MAX(RNORM,SUMR)
245 RINORM = MAX(RINORM,SUMRI)
246 END DO
247 ELSE IF ( LSAMEN( 2, C2, 'GE' ) .OR. LSAMEN( 2, C2, 'GB' ) )
248 $ THEN
249 DO I = 1, N
250 SUMR = 0
251 SUMRI = 0
252 DO J = 1, N
253 SUMR = SUMR + R(I) * CABS1(A(I,J)) * C(J)
254 SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (R(J) * C(I))
255 END DO
256 RNORM = MAX(RNORM,SUMR)
257 RINORM = MAX(RINORM,SUMRI)
258 END DO
259 END IF
260
261 RNORM = RNORM / CABS1(A(1, 1))
262 RCOND = 1.0/(RNORM * RINORM)
263
264 * Calculating the R for normwise rcond.
265 DO I = 1, N
266 RINV(I) = 0.0
267 END DO
268 DO J = 1, N
269 DO I = 1, N
270 RINV(I) = RINV(I) + CABS1(A(I,J))
271 END DO
272 END DO
273
274 * Calculating the Normwise rcond.
275 RINORM = 0.0
276 DO I = 1, N
277 SUMRI = 0.0
278 DO J = 1, N
279 SUMRI = SUMRI + CABS1(INVHILB(I,J) * RINV(J))
280 END DO
281 RINORM = MAX(RINORM, SUMRI)
282 END DO
283
284 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
285 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
286 NCOND = CABS1(A(1,1)) / RINORM
287
288 CONDTHRESH = M * EPS
289 ERRTHRESH = M * EPS
290
291 DO K = 1, NRHS
292 NORMT = 0.0
293 NORMDIF = 0.0
294 CWISE_ERR = 0.0
295 DO I = 1, N
296 NORMT = MAX(CABS1(INVHILB(I, K)), NORMT)
297 NORMDIF = MAX(CABS1(X(I,K) - INVHILB(I,K)), NORMDIF)
298 IF (INVHILB(I,K) .NE. 0.0) THEN
299 CWISE_ERR = MAX(CABS1(X(I,K) - INVHILB(I,K))
300 $ /CABS1(INVHILB(I,K)), CWISE_ERR)
301 ELSE IF (X(I, K) .NE. 0.0) THEN
302 CWISE_ERR = SLAMCH('OVERFLOW')
303 END IF
304 END DO
305 IF (NORMT .NE. 0.0) THEN
306 NWISE_ERR = NORMDIF / NORMT
307 ELSE IF (NORMDIF .NE. 0.0) THEN
308 NWISE_ERR = SLAMCH('OVERFLOW')
309 ELSE
310 NWISE_ERR = 0.0
311 ENDIF
312
313 DO I = 1, N
314 RINV(I) = 0.0
315 END DO
316 DO J = 1, N
317 DO I = 1, N
318 RINV(I) = RINV(I) + CABS1(A(I, J) * INVHILB(J, K))
319 END DO
320 END DO
321 RINORM = 0.0
322 DO I = 1, N
323 SUMRI = 0.0
324 DO J = 1, N
325 SUMRI = SUMRI
326 $ + CABS1(INVHILB(I, J) * RINV(J) / INVHILB(I, K))
327 END DO
328 RINORM = MAX(RINORM, SUMRI)
329 END DO
330 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
331 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
332 CCOND = CABS1(A(1,1))/RINORM
333
334 ! Forward error bound tests
335 NWISE_BND = ERRBND_N(K + (BND_I-1)*NRHS)
336 CWISE_BND = ERRBND_C(K + (BND_I-1)*NRHS)
337 NWISE_RCOND = ERRBND_N(K + (COND_I-1)*NRHS)
338 CWISE_RCOND = ERRBND_C(K + (COND_I-1)*NRHS)
339 ! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
340 ! $ condthresh, ncond.ge.condthresh
341 ! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
342 IF (NCOND .GE. CONDTHRESH) THEN
343 NGUAR = 'YES'
344 IF (NWISE_BND .GT. ERRTHRESH) THEN
345 TSTRAT(1) = 1/(2.0*EPS)
346 ELSE
347 IF (NWISE_BND .NE. 0.0) THEN
348 TSTRAT(1) = NWISE_ERR / NWISE_BND
349 ELSE IF (NWISE_ERR .NE. 0.0) THEN
350 TSTRAT(1) = 1/(16.0*EPS)
351 ELSE
352 TSTRAT(1) = 0.0
353 END IF
354 IF (TSTRAT(1) .GT. 1.0) THEN
355 TSTRAT(1) = 1/(4.0*EPS)
356 END IF
357 END IF
358 ELSE
359 NGUAR = 'NO'
360 IF (NWISE_BND .LT. 1.0) THEN
361 TSTRAT(1) = 1/(8.0*EPS)
362 ELSE
363 TSTRAT(1) = 1.0
364 END IF
365 END IF
366 ! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
367 ! $ condthresh, ccond.ge.condthresh
368 ! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
369 IF (CCOND .GE. CONDTHRESH) THEN
370 CGUAR = 'YES'
371 IF (CWISE_BND .GT. ERRTHRESH) THEN
372 TSTRAT(2) = 1/(2.0*EPS)
373 ELSE
374 IF (CWISE_BND .NE. 0.0) THEN
375 TSTRAT(2) = CWISE_ERR / CWISE_BND
376 ELSE IF (CWISE_ERR .NE. 0.0) THEN
377 TSTRAT(2) = 1/(16.0*EPS)
378 ELSE
379 TSTRAT(2) = 0.0
380 END IF
381 IF (TSTRAT(2) .GT. 1.0) TSTRAT(2) = 1/(4.0*EPS)
382 END IF
383 ELSE
384 CGUAR = 'NO'
385 IF (CWISE_BND .LT. 1.0) THEN
386 TSTRAT(2) = 1/(8.0*EPS)
387 ELSE
388 TSTRAT(2) = 1.0
389 END IF
390 END IF
391
392 ! Backwards error test
393 TSTRAT(3) = BERR(K)/EPS
394
395 ! Condition number tests
396 TSTRAT(4) = RCOND / ORCOND
397 IF (RCOND .GE. CONDTHRESH .AND. TSTRAT(4) .LT. 1.0)
398 $ TSTRAT(4) = 1.0 / TSTRAT(4)
399
400 TSTRAT(5) = NCOND / NWISE_RCOND
401 IF (NCOND .GE. CONDTHRESH .AND. TSTRAT(5) .LT. 1.0)
402 $ TSTRAT(5) = 1.0 / TSTRAT(5)
403
404 TSTRAT(6) = CCOND / NWISE_RCOND
405 IF (CCOND .GE. CONDTHRESH .AND. TSTRAT(6) .LT. 1.0)
406 $ TSTRAT(6) = 1.0 / TSTRAT(6)
407
408 DO I = 1, NTESTS
409 IF (TSTRAT(I) .GT. THRESH) THEN
410 IF (.NOT.PRINTED_GUIDE) THEN
411 WRITE(*,*)
412 WRITE( *, 9996) 1
413 WRITE( *, 9995) 2
414 WRITE( *, 9994) 3
415 WRITE( *, 9993) 4
416 WRITE( *, 9992) 5
417 WRITE( *, 9991) 6
418 WRITE( *, 9990) 7
419 WRITE( *, 9989) 8
420 WRITE(*,*)
421 PRINTED_GUIDE = .TRUE.
422 END IF
423 WRITE( *, 9999) C2, N, K, NGUAR, CGUAR, I, TSTRAT(I)
424 NFAIL = NFAIL + 1
425 END IF
426 END DO
427 END DO
428
429 c$$$ WRITE(*,*)
430 c$$$ WRITE(*,*) 'Normwise Error Bounds'
431 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i)
432 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i)
433 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i)
434 c$$$ WRITE(*,*)
435 c$$$ WRITE(*,*) 'Componentwise Error Bounds'
436 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i)
437 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i)
438 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i)
439 c$$$ print *, 'Info: ', info
440 c$$$ WRITE(*,*)
441 * WRITE(*,*) 'TSTRAT: ',TSTRAT
442
443 END DO
444
445 WRITE(*,*)
446 IF( NFAIL .GT. 0 ) THEN
447 WRITE(*,9998) C2, NFAIL, NTESTS*N+N_AUX_TESTS
448 ELSE
449 WRITE(*,9997) C2
450 END IF
451 9999 FORMAT( ' C', A2, 'SVXX: N =', I2, ', RHS = ', I2,
452 $ ', NWISE GUAR. = ', A, ', CWISE GUAR. = ', A,
453 $ ' test(',I1,') =', G12.5 )
454 9998 FORMAT( ' C', A2, 'SVXX: ', I6, ' out of ', I6,
455 $ ' tests failed to pass the threshold' )
456 9997 FORMAT( ' C', A2, 'SVXX passed the tests of error bounds' )
457 * Test ratios.
458 9996 FORMAT( 3X, I2, ': Normwise guaranteed forward error', / 5X,
459 $ 'Guaranteed case: if norm ( abs( Xc - Xt )',
460 $ ' / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then',
461 $ / 5X,
462 $ 'ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS')
463 9995 FORMAT( 3X, I2, ': Componentwise guaranteed forward error' )
464 9994 FORMAT( 3X, I2, ': Backwards error' )
465 9993 FORMAT( 3X, I2, ': Reciprocal condition number' )
466 9992 FORMAT( 3X, I2, ': Reciprocal normwise condition number' )
467 9991 FORMAT( 3X, I2, ': Raw normwise error estimate' )
468 9990 FORMAT( 3X, I2, ': Reciprocal componentwise condition number' )
469 9989 FORMAT( 3X, I2, ': Raw componentwise error estimate' )
470
471 8000 FORMAT( ' C', A2, 'SVXX: N =', I2, ', INFO = ', I3,
472 $ ', ORCOND = ', G12.5, ', real RCOND = ', G12.5 )
473
474 END
2 IMPLICIT NONE
3 * .. Scalar Arguments ..
4 REAL THRESH
5 CHARACTER*3 PATH
6 *
7 * Purpose
8 * ======
9 *
10 * CEBCHVXX will run CGESVXX on a series of Hilbert matrices and then
11 * compare the error bounds returned by CGESVXX 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, CGESVXX 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 CGESVXX. Let RCONDc be the RCOND returned by CGESVXX
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 * 6. 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,
75 $ N_AUX_TESTS, LDAB, LDAFB
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 COMPLEX ZDUM
87
88 * .. Local Arrays ..
89 REAL TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
90 $ S(NMAX), R(NMAX),C(NMAX),RWORK(3*NMAX),
91 $ DIFF(NMAX, NMAX),
92 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3)
93 INTEGER IPIV(NMAX)
94 COMPLEX A(NMAX,NMAX),INVHILB(NMAX,NMAX),X(NMAX,NMAX),
95 $ WORK(NMAX*3*5), AF(NMAX, NMAX),B(NMAX, NMAX),
96 $ ACOPY(NMAX, NMAX),
97 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
98 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
99 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX )
100
101 * .. External Functions ..
102 REAL SLAMCH
103
104 * .. External Subroutines ..
105 EXTERNAL CLAHILB, CGESVXX, CSYSVXX, CPOSVXX,
106 $ CGBSVXX, CLACPY, LSAMEN
107 LOGICAL LSAMEN
108
109 * .. Intrinsic Functions ..
110 INTRINSIC SQRT, MAX, ABS, REAL, AIMAG
111
112 * .. Statement Functions ..
113 REAL CABS1
114 * ..
115 * .. Statement Function Definitions ..
116 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
117
118 * .. Parameters ..
119 INTEGER NWISE_I, CWISE_I
120 PARAMETER (NWISE_I = 1, CWISE_I = 1)
121 INTEGER BND_I, COND_I
122 PARAMETER (BND_I = 2, COND_I = 3)
123
124 * Create the loop to test out the Hilbert matrices
125
126 FACT = 'E'
127 UPLO = 'U'
128 TRANS = 'N'
129 EQUED = 'N'
130 EPS = SLAMCH('Epsilon')
131 NFAIL = 0
132 N_AUX_TESTS = 0
133 LDA = NMAX
134 LDAB = (NMAX-1)+(NMAX-1)+1
135 LDAFB = 2*(NMAX-1)+(NMAX-1)+1
136 C2 = PATH( 2: 3 )
137
138 * Main loop to test the different Hilbert Matrices.
139
140 printed_guide = .false.
141
142 DO N = 1 , NMAX
143 PARAMS(1) = -1
144 PARAMS(2) = -1
145
146 KL = N-1
147 KU = N-1
148 NRHS = n
149 M = MAX(SQRT(REAL(N)), 10.0)
150
151 * Generate the Hilbert matrix, its inverse, and the
152 * right hand side, all scaled by the LCM(1,..,2N-1).
153 CALL CLAHILB(N, N, A, LDA, INVHILB, LDA, B,
154 $ LDA, WORK, INFO, PATH)
155
156 * Copy A into ACOPY.
157 CALL CLACPY('ALL', N, N, A, NMAX, ACOPY, NMAX)
158
159 * Store A in band format for GB tests
160 DO J = 1, N
161 DO I = 1, KL+KU+1
162 AB( I, J ) = (0.0E+0,0.0E+0)
163 END DO
164 END DO
165 DO J = 1, N
166 DO I = MAX( 1, J-KU ), MIN( N, J+KL )
167 AB( KU+1+I-J, J ) = A( I, J )
168 END DO
169 END DO
170
171 * Copy AB into ABCOPY.
172 DO J = 1, N
173 DO I = 1, KL+KU+1
174 ABCOPY( I, J ) = (0.0E+0,0.0E+0)
175 END DO
176 END DO
177 CALL CLACPY('ALL', KL+KU+1, N, AB, LDAB, ABCOPY, LDAB)
178
179 * Call C**SVXX with default PARAMS and N_ERR_BND = 3.
180 IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
181 CALL CSYSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
182 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
183 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
184 $ PARAMS, WORK, RWORK, INFO)
185 ELSE IF ( LSAMEN( 2, C2, 'PO' ) ) THEN
186 CALL CPOSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
187 $ EQUED, S, B, LDA, X, LDA, ORCOND,
188 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
189 $ PARAMS, WORK, RWORK, INFO)
190 ELSE IF ( LSAMEN( 2, C2, 'HE' ) ) THEN
191 CALL CHESVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
192 $ IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
193 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
194 $ PARAMS, WORK, RWORK, INFO)
195 ELSE IF ( LSAMEN( 2, C2, 'GB' ) ) THEN
196 CALL CGBSVXX(FACT, TRANS, N, KL, KU, NRHS, ABCOPY,
197 $ LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B,
198 $ LDA, X, LDA, ORCOND, RPVGRW, BERR, NERRBND,
199 $ ERRBND_N, ERRBND_C, NPARAMS, PARAMS, WORK, RWORK,
200 $ INFO)
201 ELSE
202 CALL CGESVXX(FACT, TRANS, N, NRHS, ACOPY, LDA, AF, LDA,
203 $ IPIV, EQUED, R, C, B, LDA, X, LDA, ORCOND,
204 $ RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
205 $ PARAMS, WORK, RWORK, INFO)
206 END IF
207
208 N_AUX_TESTS = N_AUX_TESTS + 1
209 IF (ORCOND .LT. EPS) THEN
210 ! Either factorization failed or the matrix is flagged, and 1 <=
211 ! INFO <= N+1. We don't decide based on rcond anymore.
212 ! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN
213 ! NFAIL = NFAIL + 1
214 ! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND
215 ! END IF
216 ELSE
217 ! Either everything succeeded (INFO == 0) or some solution failed
218 ! to converge (INFO > N+1).
219 IF (INFO .GT. 0 .AND. INFO .LE. N+1) THEN
220 NFAIL = NFAIL + 1
221 WRITE (*, FMT=8000) C2, N, INFO, ORCOND, RCOND
222 END IF
223 END IF
224
225 * Calculating the difference between C**SVXX's X and the true X.
226 DO I = 1,N
227 DO J =1,NRHS
228 DIFF(I,J) = X(I,J) - INVHILB(I,J)
229 END DO
230 END DO
231
232 * Calculating the RCOND
233 RNORM = 0
234 RINORM = 0
235 IF ( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'SY' ) .OR.
236 $ LSAMEN( 2, C2, 'HE' ) ) THEN
237 DO I = 1, N
238 SUMR = 0
239 SUMRI = 0
240 DO J = 1, N
241 SUMR = SUMR + S(I) * CABS1(A(I,J)) * S(J)
242 SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (S(J) * S(I))
243 END DO
244 RNORM = MAX(RNORM,SUMR)
245 RINORM = MAX(RINORM,SUMRI)
246 END DO
247 ELSE IF ( LSAMEN( 2, C2, 'GE' ) .OR. LSAMEN( 2, C2, 'GB' ) )
248 $ THEN
249 DO I = 1, N
250 SUMR = 0
251 SUMRI = 0
252 DO J = 1, N
253 SUMR = SUMR + R(I) * CABS1(A(I,J)) * C(J)
254 SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (R(J) * C(I))
255 END DO
256 RNORM = MAX(RNORM,SUMR)
257 RINORM = MAX(RINORM,SUMRI)
258 END DO
259 END IF
260
261 RNORM = RNORM / CABS1(A(1, 1))
262 RCOND = 1.0/(RNORM * RINORM)
263
264 * Calculating the R for normwise rcond.
265 DO I = 1, N
266 RINV(I) = 0.0
267 END DO
268 DO J = 1, N
269 DO I = 1, N
270 RINV(I) = RINV(I) + CABS1(A(I,J))
271 END DO
272 END DO
273
274 * Calculating the Normwise rcond.
275 RINORM = 0.0
276 DO I = 1, N
277 SUMRI = 0.0
278 DO J = 1, N
279 SUMRI = SUMRI + CABS1(INVHILB(I,J) * RINV(J))
280 END DO
281 RINORM = MAX(RINORM, SUMRI)
282 END DO
283
284 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
285 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
286 NCOND = CABS1(A(1,1)) / RINORM
287
288 CONDTHRESH = M * EPS
289 ERRTHRESH = M * EPS
290
291 DO K = 1, NRHS
292 NORMT = 0.0
293 NORMDIF = 0.0
294 CWISE_ERR = 0.0
295 DO I = 1, N
296 NORMT = MAX(CABS1(INVHILB(I, K)), NORMT)
297 NORMDIF = MAX(CABS1(X(I,K) - INVHILB(I,K)), NORMDIF)
298 IF (INVHILB(I,K) .NE. 0.0) THEN
299 CWISE_ERR = MAX(CABS1(X(I,K) - INVHILB(I,K))
300 $ /CABS1(INVHILB(I,K)), CWISE_ERR)
301 ELSE IF (X(I, K) .NE. 0.0) THEN
302 CWISE_ERR = SLAMCH('OVERFLOW')
303 END IF
304 END DO
305 IF (NORMT .NE. 0.0) THEN
306 NWISE_ERR = NORMDIF / NORMT
307 ELSE IF (NORMDIF .NE. 0.0) THEN
308 NWISE_ERR = SLAMCH('OVERFLOW')
309 ELSE
310 NWISE_ERR = 0.0
311 ENDIF
312
313 DO I = 1, N
314 RINV(I) = 0.0
315 END DO
316 DO J = 1, N
317 DO I = 1, N
318 RINV(I) = RINV(I) + CABS1(A(I, J) * INVHILB(J, K))
319 END DO
320 END DO
321 RINORM = 0.0
322 DO I = 1, N
323 SUMRI = 0.0
324 DO J = 1, N
325 SUMRI = SUMRI
326 $ + CABS1(INVHILB(I, J) * RINV(J) / INVHILB(I, K))
327 END DO
328 RINORM = MAX(RINORM, SUMRI)
329 END DO
330 ! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
331 ! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
332 CCOND = CABS1(A(1,1))/RINORM
333
334 ! Forward error bound tests
335 NWISE_BND = ERRBND_N(K + (BND_I-1)*NRHS)
336 CWISE_BND = ERRBND_C(K + (BND_I-1)*NRHS)
337 NWISE_RCOND = ERRBND_N(K + (COND_I-1)*NRHS)
338 CWISE_RCOND = ERRBND_C(K + (COND_I-1)*NRHS)
339 ! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
340 ! $ condthresh, ncond.ge.condthresh
341 ! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
342 IF (NCOND .GE. CONDTHRESH) THEN
343 NGUAR = 'YES'
344 IF (NWISE_BND .GT. ERRTHRESH) THEN
345 TSTRAT(1) = 1/(2.0*EPS)
346 ELSE
347 IF (NWISE_BND .NE. 0.0) THEN
348 TSTRAT(1) = NWISE_ERR / NWISE_BND
349 ELSE IF (NWISE_ERR .NE. 0.0) THEN
350 TSTRAT(1) = 1/(16.0*EPS)
351 ELSE
352 TSTRAT(1) = 0.0
353 END IF
354 IF (TSTRAT(1) .GT. 1.0) THEN
355 TSTRAT(1) = 1/(4.0*EPS)
356 END IF
357 END IF
358 ELSE
359 NGUAR = 'NO'
360 IF (NWISE_BND .LT. 1.0) THEN
361 TSTRAT(1) = 1/(8.0*EPS)
362 ELSE
363 TSTRAT(1) = 1.0
364 END IF
365 END IF
366 ! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
367 ! $ condthresh, ccond.ge.condthresh
368 ! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
369 IF (CCOND .GE. CONDTHRESH) THEN
370 CGUAR = 'YES'
371 IF (CWISE_BND .GT. ERRTHRESH) THEN
372 TSTRAT(2) = 1/(2.0*EPS)
373 ELSE
374 IF (CWISE_BND .NE. 0.0) THEN
375 TSTRAT(2) = CWISE_ERR / CWISE_BND
376 ELSE IF (CWISE_ERR .NE. 0.0) THEN
377 TSTRAT(2) = 1/(16.0*EPS)
378 ELSE
379 TSTRAT(2) = 0.0
380 END IF
381 IF (TSTRAT(2) .GT. 1.0) TSTRAT(2) = 1/(4.0*EPS)
382 END IF
383 ELSE
384 CGUAR = 'NO'
385 IF (CWISE_BND .LT. 1.0) THEN
386 TSTRAT(2) = 1/(8.0*EPS)
387 ELSE
388 TSTRAT(2) = 1.0
389 END IF
390 END IF
391
392 ! Backwards error test
393 TSTRAT(3) = BERR(K)/EPS
394
395 ! Condition number tests
396 TSTRAT(4) = RCOND / ORCOND
397 IF (RCOND .GE. CONDTHRESH .AND. TSTRAT(4) .LT. 1.0)
398 $ TSTRAT(4) = 1.0 / TSTRAT(4)
399
400 TSTRAT(5) = NCOND / NWISE_RCOND
401 IF (NCOND .GE. CONDTHRESH .AND. TSTRAT(5) .LT. 1.0)
402 $ TSTRAT(5) = 1.0 / TSTRAT(5)
403
404 TSTRAT(6) = CCOND / NWISE_RCOND
405 IF (CCOND .GE. CONDTHRESH .AND. TSTRAT(6) .LT. 1.0)
406 $ TSTRAT(6) = 1.0 / TSTRAT(6)
407
408 DO I = 1, NTESTS
409 IF (TSTRAT(I) .GT. THRESH) THEN
410 IF (.NOT.PRINTED_GUIDE) THEN
411 WRITE(*,*)
412 WRITE( *, 9996) 1
413 WRITE( *, 9995) 2
414 WRITE( *, 9994) 3
415 WRITE( *, 9993) 4
416 WRITE( *, 9992) 5
417 WRITE( *, 9991) 6
418 WRITE( *, 9990) 7
419 WRITE( *, 9989) 8
420 WRITE(*,*)
421 PRINTED_GUIDE = .TRUE.
422 END IF
423 WRITE( *, 9999) C2, N, K, NGUAR, CGUAR, I, TSTRAT(I)
424 NFAIL = NFAIL + 1
425 END IF
426 END DO
427 END DO
428
429 c$$$ WRITE(*,*)
430 c$$$ WRITE(*,*) 'Normwise Error Bounds'
431 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i)
432 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i)
433 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i)
434 c$$$ WRITE(*,*)
435 c$$$ WRITE(*,*) 'Componentwise Error Bounds'
436 c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i)
437 c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i)
438 c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i)
439 c$$$ print *, 'Info: ', info
440 c$$$ WRITE(*,*)
441 * WRITE(*,*) 'TSTRAT: ',TSTRAT
442
443 END DO
444
445 WRITE(*,*)
446 IF( NFAIL .GT. 0 ) THEN
447 WRITE(*,9998) C2, NFAIL, NTESTS*N+N_AUX_TESTS
448 ELSE
449 WRITE(*,9997) C2
450 END IF
451 9999 FORMAT( ' C', A2, 'SVXX: N =', I2, ', RHS = ', I2,
452 $ ', NWISE GUAR. = ', A, ', CWISE GUAR. = ', A,
453 $ ' test(',I1,') =', G12.5 )
454 9998 FORMAT( ' C', A2, 'SVXX: ', I6, ' out of ', I6,
455 $ ' tests failed to pass the threshold' )
456 9997 FORMAT( ' C', A2, 'SVXX passed the tests of error bounds' )
457 * Test ratios.
458 9996 FORMAT( 3X, I2, ': Normwise guaranteed forward error', / 5X,
459 $ 'Guaranteed case: if norm ( abs( Xc - Xt )',
460 $ ' / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then',
461 $ / 5X,
462 $ 'ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS')
463 9995 FORMAT( 3X, I2, ': Componentwise guaranteed forward error' )
464 9994 FORMAT( 3X, I2, ': Backwards error' )
465 9993 FORMAT( 3X, I2, ': Reciprocal condition number' )
466 9992 FORMAT( 3X, I2, ': Reciprocal normwise condition number' )
467 9991 FORMAT( 3X, I2, ': Raw normwise error estimate' )
468 9990 FORMAT( 3X, I2, ': Reciprocal componentwise condition number' )
469 9989 FORMAT( 3X, I2, ': Raw componentwise error estimate' )
470
471 8000 FORMAT( ' C', A2, 'SVXX: N =', I2, ', INFO = ', I3,
472 $ ', ORCOND = ', G12.5, ', real RCOND = ', G12.5 )
473
474 END