1 SUBROUTINE DCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
2 $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
3 $ XACT, WORK, RWORK, IWORK, NOUT )
4 *
5 * -- LAPACK test routine (version 3.3.0) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2010
8 *
9 * .. Scalar Arguments ..
10 LOGICAL TSTERR
11 INTEGER NMAX, NN, NNB, NNS, NOUT
12 DOUBLE PRECISION THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL DOTYPE( * )
16 INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
17 DOUBLE PRECISION A( * ), AFAC( * ), AINV( * ), B( * ),
18 $ RWORK( * ), WORK( * ), X( * ), XACT( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DCHKSY tests DSYTRF, -TRI2, -TRS, -TRS2, -RFS, and -CON.
25 *
26 * Arguments
27 * =========
28 *
29 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
30 * The matrix types to be used for testing. Matrices of type j
31 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
32 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
33 *
34 * NN (input) INTEGER
35 * The number of values of N contained in the vector NVAL.
36 *
37 * NVAL (input) INTEGER array, dimension (NN)
38 * The values of the matrix dimension N.
39 *
40 * NNB (input) INTEGER
41 * The number of values of NB contained in the vector NBVAL.
42 *
43 * NBVAL (input) INTEGER array, dimension (NBVAL)
44 * The values of the blocksize NB.
45 *
46 * NNS (input) INTEGER
47 * The number of values of NRHS contained in the vector NSVAL.
48 *
49 * NSVAL (input) INTEGER array, dimension (NNS)
50 * The values of the number of right hand sides NRHS.
51 *
52 * THRESH (input) DOUBLE PRECISION
53 * The threshold value for the test ratios. A result is
54 * included in the output file if RESULT >= THRESH. To have
55 * every test ratio printed, use THRESH = 0.
56 *
57 * TSTERR (input) LOGICAL
58 * Flag that indicates whether error exits are to be tested.
59 *
60 * NMAX (input) INTEGER
61 * The maximum value permitted for N, used in dimensioning the
62 * work arrays.
63 *
64 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
65 *
66 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
67 *
68 * AINV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
69 *
70 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
71 * where NSMAX is the largest entry in NSVAL.
72 *
73 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
74 *
75 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
76 *
77 * WORK (workspace) DOUBLE PRECISION array, dimension
78 * (NMAX*max(3,NSMAX))
79 *
80 * RWORK (workspace) DOUBLE PRECISION array, dimension
81 * (max(NMAX,2*NSMAX))
82 *
83 * IWORK (workspace) INTEGER array, dimension (2*NMAX)
84 *
85 * NOUT (input) INTEGER
86 * The unit number for output.
87 *
88 * =====================================================================
89 *
90 * .. Parameters ..
91 DOUBLE PRECISION ZERO
92 PARAMETER ( ZERO = 0.0D+0 )
93 INTEGER NTYPES
94 PARAMETER ( NTYPES = 10 )
95 INTEGER NTESTS
96 PARAMETER ( NTESTS = 9 )
97 * ..
98 * .. Local Scalars ..
99 LOGICAL TRFCON, ZEROT
100 CHARACTER DIST, TYPE, UPLO, XTYPE
101 CHARACTER*3 PATH
102 INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS,
103 $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE,
104 $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT
105 DOUBLE PRECISION ANORM, CNDNUM, RCOND, RCONDC
106 * ..
107 * .. Local Arrays ..
108 CHARACTER UPLOS( 2 )
109 INTEGER ISEED( 4 ), ISEEDY( 4 )
110 DOUBLE PRECISION RESULT( NTESTS )
111 DOUBLE PRECISION MYWORK( NTESTS )
112 * ..
113 * .. External Functions ..
114 DOUBLE PRECISION DGET06, DLANSY
115 EXTERNAL DGET06, DLANSY
116 * ..
117 * .. External Subroutines ..
118 EXTERNAL ALAERH, ALAHD, ALASUM, DERRSY, DGET04, DLACPY,
119 $ DLARHS, DLATB4, DLATMS, DPOT02, DPOT03, DPOT05,
120 $ DSYCON, DSYRFS, DSYT01, DSYTRF,
121 $ DSYTRI2, DSYTRS, DSYTRS2, XLAENV
122 * ..
123 * .. Intrinsic Functions ..
124 INTRINSIC MAX, MIN
125 * ..
126 * .. Scalars in Common ..
127 LOGICAL LERR, OK
128 CHARACTER*32 SRNAMT
129 INTEGER INFOT, NUNIT
130 * ..
131 * .. Common blocks ..
132 COMMON / INFOC / INFOT, NUNIT, OK, LERR
133 COMMON / SRNAMC / SRNAMT
134 * ..
135 * .. Data statements ..
136 DATA ISEEDY / 1988, 1989, 1990, 1991 /
137 DATA UPLOS / 'U', 'L' /
138 * ..
139 * .. Executable Statements ..
140 *
141 * Initialize constants and the random number seed.
142 *
143 PATH( 1: 1 ) = 'Double precision'
144 PATH( 2: 3 ) = 'SY'
145 NRUN = 0
146 NFAIL = 0
147 NERRS = 0
148 DO 10 I = 1, 4
149 ISEED( I ) = ISEEDY( I )
150 10 CONTINUE
151 *
152 * Test the error exits
153 *
154 IF( TSTERR )
155 $ CALL DERRSY( PATH, NOUT )
156 INFOT = 0
157 CALL XLAENV( 2, 2 )
158 *
159 * Do for each value of N in NVAL
160 *
161 DO 180 IN = 1, NN
162 N = NVAL( IN )
163 LDA = MAX( N, 1 )
164 XTYPE = 'N'
165 NIMAT = NTYPES
166 IF( N.LE.0 )
167 $ NIMAT = 1
168 *
169 IZERO = 0
170 DO 170 IMAT = 1, NIMAT
171 *
172 * Do the tests only if DOTYPE( IMAT ) is true.
173 *
174 IF( .NOT.DOTYPE( IMAT ) )
175 $ GO TO 170
176 *
177 * Skip types 3, 4, 5, or 6 if the matrix size is too small.
178 *
179 ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
180 IF( ZEROT .AND. N.LT.IMAT-2 )
181 $ GO TO 170
182 *
183 * Do first for UPLO = 'U', then for UPLO = 'L'
184 *
185 DO 160 IUPLO = 1, 2
186 UPLO = UPLOS( IUPLO )
187 *
188 * Set up parameters with DLATB4 and generate a test matrix
189 * with DLATMS.
190 *
191 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
192 $ CNDNUM, DIST )
193 *
194 SRNAMT = 'DLATMS'
195 CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
196 $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
197 $ INFO )
198 *
199 * Check error code from DLATMS.
200 *
201 IF( INFO.NE.0 ) THEN
202 CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1,
203 $ -1, -1, IMAT, NFAIL, NERRS, NOUT )
204 GO TO 160
205 END IF
206 *
207 * For types 3-6, zero one or more rows and columns of
208 * the matrix to test that INFO is returned correctly.
209 *
210 IF( ZEROT ) THEN
211 IF( IMAT.EQ.3 ) THEN
212 IZERO = 1
213 ELSE IF( IMAT.EQ.4 ) THEN
214 IZERO = N
215 ELSE
216 IZERO = N / 2 + 1
217 END IF
218 *
219 IF( IMAT.LT.6 ) THEN
220 *
221 * Set row and column IZERO to zero.
222 *
223 IF( IUPLO.EQ.1 ) THEN
224 IOFF = ( IZERO-1 )*LDA
225 DO 20 I = 1, IZERO - 1
226 A( IOFF+I ) = ZERO
227 20 CONTINUE
228 IOFF = IOFF + IZERO
229 DO 30 I = IZERO, N
230 A( IOFF ) = ZERO
231 IOFF = IOFF + LDA
232 30 CONTINUE
233 ELSE
234 IOFF = IZERO
235 DO 40 I = 1, IZERO - 1
236 A( IOFF ) = ZERO
237 IOFF = IOFF + LDA
238 40 CONTINUE
239 IOFF = IOFF - IZERO
240 DO 50 I = IZERO, N
241 A( IOFF+I ) = ZERO
242 50 CONTINUE
243 END IF
244 ELSE
245 IOFF = 0
246 IF( IUPLO.EQ.1 ) THEN
247 *
248 * Set the first IZERO rows and columns to zero.
249 *
250 DO 70 J = 1, N
251 I2 = MIN( J, IZERO )
252 DO 60 I = 1, I2
253 A( IOFF+I ) = ZERO
254 60 CONTINUE
255 IOFF = IOFF + LDA
256 70 CONTINUE
257 ELSE
258 *
259 * Set the last IZERO rows and columns to zero.
260 *
261 DO 90 J = 1, N
262 I1 = MAX( J, IZERO )
263 DO 80 I = I1, N
264 A( IOFF+I ) = ZERO
265 80 CONTINUE
266 IOFF = IOFF + LDA
267 90 CONTINUE
268 END IF
269 END IF
270 ELSE
271 IZERO = 0
272 END IF
273 *
274 * Do for each value of NB in NBVAL
275 *
276 DO 150 INB = 1, NNB
277 NB = NBVAL( INB )
278 CALL XLAENV( 1, NB )
279 *
280 * Compute the L*D*L' or U*D*U' factorization of the
281 * matrix.
282 *
283 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
284 LWORK = MAX( 2, NB )*LDA
285 SRNAMT = 'DSYTRF'
286 CALL DSYTRF( UPLO, N, AFAC, LDA, IWORK, AINV, LWORK,
287 $ INFO )
288 *
289 * Adjust the expected value of INFO to account for
290 * pivoting.
291 *
292 K = IZERO
293 IF( K.GT.0 ) THEN
294 100 CONTINUE
295 IF( IWORK( K ).LT.0 ) THEN
296 IF( IWORK( K ).NE.-K ) THEN
297 K = -IWORK( K )
298 GO TO 100
299 END IF
300 ELSE IF( IWORK( K ).NE.K ) THEN
301 K = IWORK( K )
302 GO TO 100
303 END IF
304 END IF
305 *
306 * Check error code from DSYTRF.
307 *
308 IF( INFO.NE.K )
309 $ CALL ALAERH( PATH, 'DSYTRF', INFO, K, UPLO, N, N,
310 $ -1, -1, NB, IMAT, NFAIL, NERRS, NOUT )
311 IF( INFO.NE.0 ) THEN
312 TRFCON = .TRUE.
313 ELSE
314 TRFCON = .FALSE.
315 END IF
316 *
317 *+ TEST 1
318 * Reconstruct matrix from factors and compute residual.
319 *
320 CALL DSYT01( UPLO, N, A, LDA, AFAC, LDA, IWORK, AINV,
321 $ LDA, RWORK, RESULT( 1 ) )
322 NT = 1
323 *
324 *+ TEST 2
325 * Form the inverse and compute the residual.
326 *
327 IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN
328 CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
329 SRNAMT = 'DSYTRI2'
330 LWORK = (N+NB+1)*(NB+3)
331 CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
332 $ LWORK, INFO )
333 *
334 * Check error code from DSYTRI2.
335 *
336 IF( INFO.NE.0 )
337 $ CALL ALAERH( PATH, 'DSYTRI2', INFO, -1, UPLO, N,
338 $ N, -1, -1, -1, IMAT, NFAIL, NERRS,
339 $ NOUT )
340 *
341 CALL DPOT03( UPLO, N, A, LDA, AINV, LDA, WORK, LDA,
342 $ RWORK, RCONDC, RESULT( 2 ) )
343 NT = 2
344 END IF
345 *
346 * Print information about the tests that did not pass
347 * the threshold.
348 *
349 DO 110 K = 1, NT
350 IF( RESULT( K ).GE.THRESH ) THEN
351 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
352 $ CALL ALAHD( NOUT, PATH )
353 WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K,
354 $ RESULT( K )
355 NFAIL = NFAIL + 1
356 END IF
357 110 CONTINUE
358 NRUN = NRUN + NT
359 *
360 * Skip the other tests if this is not the first block
361 * size.
362 *
363 IF( INB.GT.1 )
364 $ GO TO 150
365 *
366 * Do only the condition estimate if INFO is not 0.
367 *
368 IF( TRFCON ) THEN
369 RCONDC = ZERO
370 GO TO 140
371 END IF
372 *
373 DO 130 IRHS = 1, NNS
374 NRHS = NSVAL( IRHS )
375 *
376 *+ TEST 3 ( Using TRS)
377 * Solve and compute residual for A * X = B.
378 *
379 SRNAMT = 'DLARHS'
380 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
381 $ NRHS, A, LDA, XACT, LDA, B, LDA,
382 $ ISEED, INFO )
383 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
384 *
385 SRNAMT = 'DSYTRS'
386 CALL DSYTRS( UPLO, N, NRHS, AFAC, LDA, IWORK, X,
387 $ LDA, INFO )
388 *
389 * Check error code from DSYTRS.
390 *
391 IF( INFO.NE.0 )
392 $ CALL ALAERH( PATH, 'DSYTRS', INFO, 0, UPLO, N,
393 $ N, -1, -1, NRHS, IMAT, NFAIL,
394 $ NERRS, NOUT )
395 *
396 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
397 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
398 $ LDA, RWORK, RESULT( 3 ) )
399 *
400 *+ TEST 4 (Using TRS2)
401 *
402 * Solve and compute residual for A * X = B.
403 *
404 SRNAMT = 'DLARHS'
405 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
406 $ NRHS, A, LDA, XACT, LDA, B, LDA,
407 $ ISEED, INFO )
408 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
409 *
410 SRNAMT = 'DSYTRS2'
411 CALL DSYTRS2( UPLO, N, NRHS, AFAC, LDA, IWORK, X,
412 $ LDA, WORK, INFO )
413 *
414 * Check error code from DSYTRS2.
415 *
416 IF( INFO.NE.0 )
417 $ CALL ALAERH( PATH, 'DSYTRS2', INFO, 0, UPLO, N,
418 $ N, -1, -1, NRHS, IMAT, NFAIL,
419 $ NERRS, NOUT )
420 *
421 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
422 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
423 $ LDA, RWORK, RESULT( 4 ) )
424 *
425 *+ TEST 5
426 * Check solution from generated exact solution.
427 *
428 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
429 $ RESULT( 5 ) )
430 *
431 *+ TESTS 6, 7, and 8
432 * Use iterative refinement to improve the solution.
433 *
434 SRNAMT = 'DSYRFS'
435 CALL DSYRFS( UPLO, N, NRHS, A, LDA, AFAC, LDA,
436 $ IWORK, B, LDA, X, LDA, RWORK,
437 $ RWORK( NRHS+1 ), WORK, IWORK( N+1 ),
438 $ INFO )
439 *
440 * Check error code from DSYRFS.
441 *
442 IF( INFO.NE.0 )
443 $ CALL ALAERH( PATH, 'DSYRFS', INFO, 0, UPLO, N,
444 $ N, -1, -1, NRHS, IMAT, NFAIL,
445 $ NERRS, NOUT )
446 *
447 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
448 $ RESULT( 6 ) )
449 CALL DPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
450 $ XACT, LDA, RWORK, RWORK( NRHS+1 ),
451 $ RESULT( 7 ) )
452 *
453 * Print information about the tests that did not pass
454 * the threshold.
455 *
456 DO 120 K = 3, 8
457 IF( RESULT( K ).GE.THRESH ) THEN
458 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
459 $ CALL ALAHD( NOUT, PATH )
460 WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS,
461 $ IMAT, K, RESULT( K )
462 NFAIL = NFAIL + 1
463 END IF
464 120 CONTINUE
465 NRUN = NRUN + 5
466 130 CONTINUE
467 *
468 *+ TEST 9
469 * Get an estimate of RCOND = 1/CNDNUM.
470 *
471 140 CONTINUE
472 ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK )
473 SRNAMT = 'DSYCON'
474 CALL DSYCON( UPLO, N, AFAC, LDA, IWORK, ANORM, RCOND,
475 $ WORK, IWORK( N+1 ), INFO )
476 *
477 * Check error code from DSYCON.
478 *
479 IF( INFO.NE.0 )
480 $ CALL ALAERH( PATH, 'DSYCON', INFO, 0, UPLO, N, N,
481 $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
482 *
483 RESULT( 9 ) = DGET06( RCOND, RCONDC )
484 *
485 * Print information about the tests that did not pass
486 * the threshold.
487 *
488 IF( RESULT( 9 ).GE.THRESH ) THEN
489 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
490 $ CALL ALAHD( NOUT, PATH )
491 WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 9,
492 $ RESULT( 9 )
493 NFAIL = NFAIL + 1
494 END IF
495 NRUN = NRUN + 1
496 150 CONTINUE
497 *
498 160 CONTINUE
499 170 CONTINUE
500 180 CONTINUE
501 *
502 * Print a summary of the results.
503 *
504 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
505 *
506 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ',
507 $ I2, ', test ', I2, ', ratio =', G12.5 )
508 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
509 $ I2, ', test(', I2, ') =', G12.5 )
510 9997 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ',', 10X, ' type ', I2,
511 $ ', test(', I2, ') =', G12.5 )
512 RETURN
513 *
514 * End of DCHKSY
515 *
516 END
2 $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
3 $ XACT, WORK, RWORK, IWORK, NOUT )
4 *
5 * -- LAPACK test routine (version 3.3.0) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2010
8 *
9 * .. Scalar Arguments ..
10 LOGICAL TSTERR
11 INTEGER NMAX, NN, NNB, NNS, NOUT
12 DOUBLE PRECISION THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL DOTYPE( * )
16 INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
17 DOUBLE PRECISION A( * ), AFAC( * ), AINV( * ), B( * ),
18 $ RWORK( * ), WORK( * ), X( * ), XACT( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DCHKSY tests DSYTRF, -TRI2, -TRS, -TRS2, -RFS, and -CON.
25 *
26 * Arguments
27 * =========
28 *
29 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
30 * The matrix types to be used for testing. Matrices of type j
31 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
32 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
33 *
34 * NN (input) INTEGER
35 * The number of values of N contained in the vector NVAL.
36 *
37 * NVAL (input) INTEGER array, dimension (NN)
38 * The values of the matrix dimension N.
39 *
40 * NNB (input) INTEGER
41 * The number of values of NB contained in the vector NBVAL.
42 *
43 * NBVAL (input) INTEGER array, dimension (NBVAL)
44 * The values of the blocksize NB.
45 *
46 * NNS (input) INTEGER
47 * The number of values of NRHS contained in the vector NSVAL.
48 *
49 * NSVAL (input) INTEGER array, dimension (NNS)
50 * The values of the number of right hand sides NRHS.
51 *
52 * THRESH (input) DOUBLE PRECISION
53 * The threshold value for the test ratios. A result is
54 * included in the output file if RESULT >= THRESH. To have
55 * every test ratio printed, use THRESH = 0.
56 *
57 * TSTERR (input) LOGICAL
58 * Flag that indicates whether error exits are to be tested.
59 *
60 * NMAX (input) INTEGER
61 * The maximum value permitted for N, used in dimensioning the
62 * work arrays.
63 *
64 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
65 *
66 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
67 *
68 * AINV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
69 *
70 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
71 * where NSMAX is the largest entry in NSVAL.
72 *
73 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
74 *
75 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
76 *
77 * WORK (workspace) DOUBLE PRECISION array, dimension
78 * (NMAX*max(3,NSMAX))
79 *
80 * RWORK (workspace) DOUBLE PRECISION array, dimension
81 * (max(NMAX,2*NSMAX))
82 *
83 * IWORK (workspace) INTEGER array, dimension (2*NMAX)
84 *
85 * NOUT (input) INTEGER
86 * The unit number for output.
87 *
88 * =====================================================================
89 *
90 * .. Parameters ..
91 DOUBLE PRECISION ZERO
92 PARAMETER ( ZERO = 0.0D+0 )
93 INTEGER NTYPES
94 PARAMETER ( NTYPES = 10 )
95 INTEGER NTESTS
96 PARAMETER ( NTESTS = 9 )
97 * ..
98 * .. Local Scalars ..
99 LOGICAL TRFCON, ZEROT
100 CHARACTER DIST, TYPE, UPLO, XTYPE
101 CHARACTER*3 PATH
102 INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS,
103 $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE,
104 $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT
105 DOUBLE PRECISION ANORM, CNDNUM, RCOND, RCONDC
106 * ..
107 * .. Local Arrays ..
108 CHARACTER UPLOS( 2 )
109 INTEGER ISEED( 4 ), ISEEDY( 4 )
110 DOUBLE PRECISION RESULT( NTESTS )
111 DOUBLE PRECISION MYWORK( NTESTS )
112 * ..
113 * .. External Functions ..
114 DOUBLE PRECISION DGET06, DLANSY
115 EXTERNAL DGET06, DLANSY
116 * ..
117 * .. External Subroutines ..
118 EXTERNAL ALAERH, ALAHD, ALASUM, DERRSY, DGET04, DLACPY,
119 $ DLARHS, DLATB4, DLATMS, DPOT02, DPOT03, DPOT05,
120 $ DSYCON, DSYRFS, DSYT01, DSYTRF,
121 $ DSYTRI2, DSYTRS, DSYTRS2, XLAENV
122 * ..
123 * .. Intrinsic Functions ..
124 INTRINSIC MAX, MIN
125 * ..
126 * .. Scalars in Common ..
127 LOGICAL LERR, OK
128 CHARACTER*32 SRNAMT
129 INTEGER INFOT, NUNIT
130 * ..
131 * .. Common blocks ..
132 COMMON / INFOC / INFOT, NUNIT, OK, LERR
133 COMMON / SRNAMC / SRNAMT
134 * ..
135 * .. Data statements ..
136 DATA ISEEDY / 1988, 1989, 1990, 1991 /
137 DATA UPLOS / 'U', 'L' /
138 * ..
139 * .. Executable Statements ..
140 *
141 * Initialize constants and the random number seed.
142 *
143 PATH( 1: 1 ) = 'Double precision'
144 PATH( 2: 3 ) = 'SY'
145 NRUN = 0
146 NFAIL = 0
147 NERRS = 0
148 DO 10 I = 1, 4
149 ISEED( I ) = ISEEDY( I )
150 10 CONTINUE
151 *
152 * Test the error exits
153 *
154 IF( TSTERR )
155 $ CALL DERRSY( PATH, NOUT )
156 INFOT = 0
157 CALL XLAENV( 2, 2 )
158 *
159 * Do for each value of N in NVAL
160 *
161 DO 180 IN = 1, NN
162 N = NVAL( IN )
163 LDA = MAX( N, 1 )
164 XTYPE = 'N'
165 NIMAT = NTYPES
166 IF( N.LE.0 )
167 $ NIMAT = 1
168 *
169 IZERO = 0
170 DO 170 IMAT = 1, NIMAT
171 *
172 * Do the tests only if DOTYPE( IMAT ) is true.
173 *
174 IF( .NOT.DOTYPE( IMAT ) )
175 $ GO TO 170
176 *
177 * Skip types 3, 4, 5, or 6 if the matrix size is too small.
178 *
179 ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
180 IF( ZEROT .AND. N.LT.IMAT-2 )
181 $ GO TO 170
182 *
183 * Do first for UPLO = 'U', then for UPLO = 'L'
184 *
185 DO 160 IUPLO = 1, 2
186 UPLO = UPLOS( IUPLO )
187 *
188 * Set up parameters with DLATB4 and generate a test matrix
189 * with DLATMS.
190 *
191 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
192 $ CNDNUM, DIST )
193 *
194 SRNAMT = 'DLATMS'
195 CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
196 $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
197 $ INFO )
198 *
199 * Check error code from DLATMS.
200 *
201 IF( INFO.NE.0 ) THEN
202 CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1,
203 $ -1, -1, IMAT, NFAIL, NERRS, NOUT )
204 GO TO 160
205 END IF
206 *
207 * For types 3-6, zero one or more rows and columns of
208 * the matrix to test that INFO is returned correctly.
209 *
210 IF( ZEROT ) THEN
211 IF( IMAT.EQ.3 ) THEN
212 IZERO = 1
213 ELSE IF( IMAT.EQ.4 ) THEN
214 IZERO = N
215 ELSE
216 IZERO = N / 2 + 1
217 END IF
218 *
219 IF( IMAT.LT.6 ) THEN
220 *
221 * Set row and column IZERO to zero.
222 *
223 IF( IUPLO.EQ.1 ) THEN
224 IOFF = ( IZERO-1 )*LDA
225 DO 20 I = 1, IZERO - 1
226 A( IOFF+I ) = ZERO
227 20 CONTINUE
228 IOFF = IOFF + IZERO
229 DO 30 I = IZERO, N
230 A( IOFF ) = ZERO
231 IOFF = IOFF + LDA
232 30 CONTINUE
233 ELSE
234 IOFF = IZERO
235 DO 40 I = 1, IZERO - 1
236 A( IOFF ) = ZERO
237 IOFF = IOFF + LDA
238 40 CONTINUE
239 IOFF = IOFF - IZERO
240 DO 50 I = IZERO, N
241 A( IOFF+I ) = ZERO
242 50 CONTINUE
243 END IF
244 ELSE
245 IOFF = 0
246 IF( IUPLO.EQ.1 ) THEN
247 *
248 * Set the first IZERO rows and columns to zero.
249 *
250 DO 70 J = 1, N
251 I2 = MIN( J, IZERO )
252 DO 60 I = 1, I2
253 A( IOFF+I ) = ZERO
254 60 CONTINUE
255 IOFF = IOFF + LDA
256 70 CONTINUE
257 ELSE
258 *
259 * Set the last IZERO rows and columns to zero.
260 *
261 DO 90 J = 1, N
262 I1 = MAX( J, IZERO )
263 DO 80 I = I1, N
264 A( IOFF+I ) = ZERO
265 80 CONTINUE
266 IOFF = IOFF + LDA
267 90 CONTINUE
268 END IF
269 END IF
270 ELSE
271 IZERO = 0
272 END IF
273 *
274 * Do for each value of NB in NBVAL
275 *
276 DO 150 INB = 1, NNB
277 NB = NBVAL( INB )
278 CALL XLAENV( 1, NB )
279 *
280 * Compute the L*D*L' or U*D*U' factorization of the
281 * matrix.
282 *
283 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
284 LWORK = MAX( 2, NB )*LDA
285 SRNAMT = 'DSYTRF'
286 CALL DSYTRF( UPLO, N, AFAC, LDA, IWORK, AINV, LWORK,
287 $ INFO )
288 *
289 * Adjust the expected value of INFO to account for
290 * pivoting.
291 *
292 K = IZERO
293 IF( K.GT.0 ) THEN
294 100 CONTINUE
295 IF( IWORK( K ).LT.0 ) THEN
296 IF( IWORK( K ).NE.-K ) THEN
297 K = -IWORK( K )
298 GO TO 100
299 END IF
300 ELSE IF( IWORK( K ).NE.K ) THEN
301 K = IWORK( K )
302 GO TO 100
303 END IF
304 END IF
305 *
306 * Check error code from DSYTRF.
307 *
308 IF( INFO.NE.K )
309 $ CALL ALAERH( PATH, 'DSYTRF', INFO, K, UPLO, N, N,
310 $ -1, -1, NB, IMAT, NFAIL, NERRS, NOUT )
311 IF( INFO.NE.0 ) THEN
312 TRFCON = .TRUE.
313 ELSE
314 TRFCON = .FALSE.
315 END IF
316 *
317 *+ TEST 1
318 * Reconstruct matrix from factors and compute residual.
319 *
320 CALL DSYT01( UPLO, N, A, LDA, AFAC, LDA, IWORK, AINV,
321 $ LDA, RWORK, RESULT( 1 ) )
322 NT = 1
323 *
324 *+ TEST 2
325 * Form the inverse and compute the residual.
326 *
327 IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN
328 CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
329 SRNAMT = 'DSYTRI2'
330 LWORK = (N+NB+1)*(NB+3)
331 CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
332 $ LWORK, INFO )
333 *
334 * Check error code from DSYTRI2.
335 *
336 IF( INFO.NE.0 )
337 $ CALL ALAERH( PATH, 'DSYTRI2', INFO, -1, UPLO, N,
338 $ N, -1, -1, -1, IMAT, NFAIL, NERRS,
339 $ NOUT )
340 *
341 CALL DPOT03( UPLO, N, A, LDA, AINV, LDA, WORK, LDA,
342 $ RWORK, RCONDC, RESULT( 2 ) )
343 NT = 2
344 END IF
345 *
346 * Print information about the tests that did not pass
347 * the threshold.
348 *
349 DO 110 K = 1, NT
350 IF( RESULT( K ).GE.THRESH ) THEN
351 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
352 $ CALL ALAHD( NOUT, PATH )
353 WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K,
354 $ RESULT( K )
355 NFAIL = NFAIL + 1
356 END IF
357 110 CONTINUE
358 NRUN = NRUN + NT
359 *
360 * Skip the other tests if this is not the first block
361 * size.
362 *
363 IF( INB.GT.1 )
364 $ GO TO 150
365 *
366 * Do only the condition estimate if INFO is not 0.
367 *
368 IF( TRFCON ) THEN
369 RCONDC = ZERO
370 GO TO 140
371 END IF
372 *
373 DO 130 IRHS = 1, NNS
374 NRHS = NSVAL( IRHS )
375 *
376 *+ TEST 3 ( Using TRS)
377 * Solve and compute residual for A * X = B.
378 *
379 SRNAMT = 'DLARHS'
380 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
381 $ NRHS, A, LDA, XACT, LDA, B, LDA,
382 $ ISEED, INFO )
383 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
384 *
385 SRNAMT = 'DSYTRS'
386 CALL DSYTRS( UPLO, N, NRHS, AFAC, LDA, IWORK, X,
387 $ LDA, INFO )
388 *
389 * Check error code from DSYTRS.
390 *
391 IF( INFO.NE.0 )
392 $ CALL ALAERH( PATH, 'DSYTRS', INFO, 0, UPLO, N,
393 $ N, -1, -1, NRHS, IMAT, NFAIL,
394 $ NERRS, NOUT )
395 *
396 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
397 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
398 $ LDA, RWORK, RESULT( 3 ) )
399 *
400 *+ TEST 4 (Using TRS2)
401 *
402 * Solve and compute residual for A * X = B.
403 *
404 SRNAMT = 'DLARHS'
405 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
406 $ NRHS, A, LDA, XACT, LDA, B, LDA,
407 $ ISEED, INFO )
408 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
409 *
410 SRNAMT = 'DSYTRS2'
411 CALL DSYTRS2( UPLO, N, NRHS, AFAC, LDA, IWORK, X,
412 $ LDA, WORK, INFO )
413 *
414 * Check error code from DSYTRS2.
415 *
416 IF( INFO.NE.0 )
417 $ CALL ALAERH( PATH, 'DSYTRS2', INFO, 0, UPLO, N,
418 $ N, -1, -1, NRHS, IMAT, NFAIL,
419 $ NERRS, NOUT )
420 *
421 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
422 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
423 $ LDA, RWORK, RESULT( 4 ) )
424 *
425 *+ TEST 5
426 * Check solution from generated exact solution.
427 *
428 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
429 $ RESULT( 5 ) )
430 *
431 *+ TESTS 6, 7, and 8
432 * Use iterative refinement to improve the solution.
433 *
434 SRNAMT = 'DSYRFS'
435 CALL DSYRFS( UPLO, N, NRHS, A, LDA, AFAC, LDA,
436 $ IWORK, B, LDA, X, LDA, RWORK,
437 $ RWORK( NRHS+1 ), WORK, IWORK( N+1 ),
438 $ INFO )
439 *
440 * Check error code from DSYRFS.
441 *
442 IF( INFO.NE.0 )
443 $ CALL ALAERH( PATH, 'DSYRFS', INFO, 0, UPLO, N,
444 $ N, -1, -1, NRHS, IMAT, NFAIL,
445 $ NERRS, NOUT )
446 *
447 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
448 $ RESULT( 6 ) )
449 CALL DPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
450 $ XACT, LDA, RWORK, RWORK( NRHS+1 ),
451 $ RESULT( 7 ) )
452 *
453 * Print information about the tests that did not pass
454 * the threshold.
455 *
456 DO 120 K = 3, 8
457 IF( RESULT( K ).GE.THRESH ) THEN
458 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
459 $ CALL ALAHD( NOUT, PATH )
460 WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS,
461 $ IMAT, K, RESULT( K )
462 NFAIL = NFAIL + 1
463 END IF
464 120 CONTINUE
465 NRUN = NRUN + 5
466 130 CONTINUE
467 *
468 *+ TEST 9
469 * Get an estimate of RCOND = 1/CNDNUM.
470 *
471 140 CONTINUE
472 ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK )
473 SRNAMT = 'DSYCON'
474 CALL DSYCON( UPLO, N, AFAC, LDA, IWORK, ANORM, RCOND,
475 $ WORK, IWORK( N+1 ), INFO )
476 *
477 * Check error code from DSYCON.
478 *
479 IF( INFO.NE.0 )
480 $ CALL ALAERH( PATH, 'DSYCON', INFO, 0, UPLO, N, N,
481 $ -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
482 *
483 RESULT( 9 ) = DGET06( RCOND, RCONDC )
484 *
485 * Print information about the tests that did not pass
486 * the threshold.
487 *
488 IF( RESULT( 9 ).GE.THRESH ) THEN
489 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
490 $ CALL ALAHD( NOUT, PATH )
491 WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 9,
492 $ RESULT( 9 )
493 NFAIL = NFAIL + 1
494 END IF
495 NRUN = NRUN + 1
496 150 CONTINUE
497 *
498 160 CONTINUE
499 170 CONTINUE
500 180 CONTINUE
501 *
502 * Print a summary of the results.
503 *
504 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
505 *
506 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ',
507 $ I2, ', test ', I2, ', ratio =', G12.5 )
508 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
509 $ I2, ', test(', I2, ') =', G12.5 )
510 9997 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ',', 10X, ' type ', I2,
511 $ ', test(', I2, ') =', G12.5 )
512 RETURN
513 *
514 * End of DCHKSY
515 *
516 END