1 SUBROUTINE DDRVSY( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
2 $ A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK,
3 $ NOUT )
4 *
5 * -- LAPACK test routine (version 3.2) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 LOGICAL TSTERR
11 INTEGER NMAX, NN, NOUT, NRHS
12 DOUBLE PRECISION THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL DOTYPE( * )
16 INTEGER IWORK( * ), NVAL( * )
17 DOUBLE PRECISION A( * ), AFAC( * ), AINV( * ), B( * ),
18 $ RWORK( * ), WORK( * ), X( * ), XACT( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DDRVSY tests the driver routines DSYSV, -SVX, and -SVXX.
25 *
26 * Note that this file is used only when the XBLAS are available,
27 * otherwise ddrvsy.f defines this subroutine.
28 *
29 * Arguments
30 * =========
31 *
32 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
33 * The matrix types to be used for testing. Matrices of type j
34 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
35 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
36 *
37 * NN (input) INTEGER
38 * The number of values of N contained in the vector NVAL.
39 *
40 * NVAL (input) INTEGER array, dimension (NN)
41 * The values of the matrix dimension N.
42 *
43 * NRHS (input) INTEGER
44 * The number of right hand side vectors to be generated for
45 * each linear system.
46 *
47 * THRESH (input) DOUBLE PRECISION
48 * The threshold value for the test ratios. A result is
49 * included in the output file if RESULT >= THRESH. To have
50 * every test ratio printed, use THRESH = 0.
51 *
52 * TSTERR (input) LOGICAL
53 * Flag that indicates whether error exits are to be tested.
54 *
55 * NMAX (input) INTEGER
56 * The maximum value permitted for N, used in dimensioning the
57 * work arrays.
58 *
59 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
60 *
61 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
62 *
63 * AINV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
64 *
65 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS)
66 *
67 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS)
68 *
69 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS)
70 *
71 * WORK (workspace) DOUBLE PRECISION array, dimension
72 * (NMAX*max(2,NRHS))
73 *
74 * RWORK (workspace) DOUBLE PRECISION array, dimension (NMAX+2*NRHS)
75 *
76 * IWORK (workspace) INTEGER array, dimension (2*NMAX)
77 *
78 * NOUT (input) INTEGER
79 * The unit number for output.
80 *
81 * =====================================================================
82 *
83 * .. Parameters ..
84 DOUBLE PRECISION ONE, ZERO
85 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
86 INTEGER NTYPES, NTESTS
87 PARAMETER ( NTYPES = 10, NTESTS = 6 )
88 INTEGER NFACT
89 PARAMETER ( NFACT = 2 )
90 * ..
91 * .. Local Scalars ..
92 LOGICAL ZEROT
93 CHARACTER DIST, EQUED, FACT, TYPE, UPLO, XTYPE
94 CHARACTER*3 PATH
95 INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
96 $ IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N,
97 $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT,
98 $ N_ERR_BNDS
99 DOUBLE PRECISION AINVNM, ANORM, CNDNUM, RCOND, RCONDC,
100 $ RPVGRW_SVXX
101 * ..
102 * .. Local Arrays ..
103 CHARACTER FACTS( NFACT ), UPLOS( 2 )
104 INTEGER ISEED( 4 ), ISEEDY( 4 )
105 DOUBLE PRECISION RESULT( NTESTS ), BERR( NRHS ),
106 $ ERRBNDS_N( NRHS, 3 ), ERRBNDS_C( NRHS, 3 )
107 * ..
108 * .. External Functions ..
109 DOUBLE PRECISION DGET06, DLANSY
110 EXTERNAL DGET06, DLANSY
111 * ..
112 * .. External Subroutines ..
113 EXTERNAL ALADHD, ALAERH, ALASVM, DERRVX, DGET04, DLACPY,
114 $ DLARHS, DLASET, DLATB4, DLATMS, DPOT02, DPOT05,
115 $ DSYSV, DSYSVX, DSYT01, DSYTRF, DSYTRI2, XLAENV,
116 $ DSYSVXX
117 * ..
118 * .. Scalars in Common ..
119 LOGICAL LERR, OK
120 CHARACTER*32 SRNAMT
121 INTEGER INFOT, NUNIT
122 * ..
123 * .. Common blocks ..
124 COMMON / INFOC / INFOT, NUNIT, OK, LERR
125 COMMON / SRNAMC / SRNAMT
126 * ..
127 * .. Intrinsic Functions ..
128 INTRINSIC MAX, MIN
129 * ..
130 * .. Data statements ..
131 DATA ISEEDY / 1988, 1989, 1990, 1991 /
132 DATA UPLOS / 'U', 'L' / , FACTS / 'F', 'N' /
133 * ..
134 * .. Executable Statements ..
135 *
136 * Initialize constants and the random number seed.
137 *
138 PATH( 1: 1 ) = 'Double precision'
139 PATH( 2: 3 ) = 'SY'
140 NRUN = 0
141 NFAIL = 0
142 NERRS = 0
143 DO 10 I = 1, 4
144 ISEED( I ) = ISEEDY( I )
145 10 CONTINUE
146 LWORK = MAX( 2*NMAX, NMAX*NRHS )
147 *
148 * Test the error exits
149 *
150 IF( TSTERR )
151 $ CALL DERRVX( PATH, NOUT )
152 INFOT = 0
153 *
154 * Set the block size and minimum block size for testing.
155 *
156 NB = 1
157 NBMIN = 2
158 CALL XLAENV( 1, NB )
159 CALL XLAENV( 2, NBMIN )
160 *
161 * Do for each value of N in NVAL
162 *
163 DO 180 IN = 1, NN
164 N = NVAL( IN )
165 LDA = MAX( N, 1 )
166 XTYPE = 'N'
167 NIMAT = NTYPES
168 IF( N.LE.0 )
169 $ NIMAT = 1
170 *
171 DO 170 IMAT = 1, NIMAT
172 *
173 * Do the tests only if DOTYPE( IMAT ) is true.
174 *
175 IF( .NOT.DOTYPE( IMAT ) )
176 $ GO TO 170
177 *
178 * Skip types 3, 4, 5, or 6 if the matrix size is too small.
179 *
180 ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
181 IF( ZEROT .AND. N.LT.IMAT-2 )
182 $ GO TO 170
183 *
184 * Do first for UPLO = 'U', then for UPLO = 'L'
185 *
186 DO 160 IUPLO = 1, 2
187 UPLO = UPLOS( IUPLO )
188 *
189 * Set up parameters with DLATB4 and generate a test matrix
190 * with DLATMS.
191 *
192 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
193 $ CNDNUM, DIST )
194 *
195 SRNAMT = 'DLATMS'
196 CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
197 $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
198 $ INFO )
199 *
200 * Check error code from DLATMS.
201 *
202 IF( INFO.NE.0 ) THEN
203 CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1,
204 $ -1, -1, IMAT, NFAIL, NERRS, NOUT )
205 GO TO 160
206 END IF
207 *
208 * For types 3-6, zero one or more rows and columns of the
209 * matrix to test that INFO is returned correctly.
210 *
211 IF( ZEROT ) THEN
212 IF( IMAT.EQ.3 ) THEN
213 IZERO = 1
214 ELSE IF( IMAT.EQ.4 ) THEN
215 IZERO = N
216 ELSE
217 IZERO = N / 2 + 1
218 END IF
219 *
220 IF( IMAT.LT.6 ) THEN
221 *
222 * Set row and column IZERO to zero.
223 *
224 IF( IUPLO.EQ.1 ) THEN
225 IOFF = ( IZERO-1 )*LDA
226 DO 20 I = 1, IZERO - 1
227 A( IOFF+I ) = ZERO
228 20 CONTINUE
229 IOFF = IOFF + IZERO
230 DO 30 I = IZERO, N
231 A( IOFF ) = ZERO
232 IOFF = IOFF + LDA
233 30 CONTINUE
234 ELSE
235 IOFF = IZERO
236 DO 40 I = 1, IZERO - 1
237 A( IOFF ) = ZERO
238 IOFF = IOFF + LDA
239 40 CONTINUE
240 IOFF = IOFF - IZERO
241 DO 50 I = IZERO, N
242 A( IOFF+I ) = ZERO
243 50 CONTINUE
244 END IF
245 ELSE
246 IOFF = 0
247 IF( IUPLO.EQ.1 ) THEN
248 *
249 * Set the first IZERO rows and columns to zero.
250 *
251 DO 70 J = 1, N
252 I2 = MIN( J, IZERO )
253 DO 60 I = 1, I2
254 A( IOFF+I ) = ZERO
255 60 CONTINUE
256 IOFF = IOFF + LDA
257 70 CONTINUE
258 ELSE
259 *
260 * Set the last IZERO rows and columns to zero.
261 *
262 DO 90 J = 1, N
263 I1 = MAX( J, IZERO )
264 DO 80 I = I1, N
265 A( IOFF+I ) = ZERO
266 80 CONTINUE
267 IOFF = IOFF + LDA
268 90 CONTINUE
269 END IF
270 END IF
271 ELSE
272 IZERO = 0
273 END IF
274 *
275 DO 150 IFACT = 1, NFACT
276 *
277 * Do first for FACT = 'F', then for other values.
278 *
279 FACT = FACTS( IFACT )
280 *
281 * Compute the condition number for comparison with
282 * the value returned by DSYSVX.
283 *
284 IF( ZEROT ) THEN
285 IF( IFACT.EQ.1 )
286 $ GO TO 150
287 RCONDC = ZERO
288 *
289 ELSE IF( IFACT.EQ.1 ) THEN
290 *
291 * Compute the 1-norm of A.
292 *
293 ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK )
294 *
295 * Factor the matrix A.
296 *
297 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
298 CALL DSYTRF( UPLO, N, AFAC, LDA, IWORK, WORK,
299 $ LWORK, INFO )
300 *
301 * Compute inv(A) and take its norm.
302 *
303 CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
304 LWORK = (N+NB+1)*(NB+3)
305 CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
306 $ LWORK, INFO )
307 AINVNM = DLANSY( '1', UPLO, N, AINV, LDA, RWORK )
308 *
309 * Compute the 1-norm condition number of A.
310 *
311 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
312 RCONDC = ONE
313 ELSE
314 RCONDC = ( ONE / ANORM ) / AINVNM
315 END IF
316 END IF
317 *
318 * Form an exact solution and set the right hand side.
319 *
320 SRNAMT = 'DLARHS'
321 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
322 $ NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
323 $ INFO )
324 XTYPE = 'C'
325 *
326 * --- Test DSYSV ---
327 *
328 IF( IFACT.EQ.2 ) THEN
329 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
330 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
331 *
332 * Factor the matrix and solve the system using DSYSV.
333 *
334 SRNAMT = 'DSYSV '
335 CALL DSYSV( UPLO, N, NRHS, AFAC, LDA, IWORK, X,
336 $ LDA, WORK, LWORK, INFO )
337 *
338 * Adjust the expected value of INFO to account for
339 * pivoting.
340 *
341 K = IZERO
342 IF( K.GT.0 ) THEN
343 100 CONTINUE
344 IF( IWORK( K ).LT.0 ) THEN
345 IF( IWORK( K ).NE.-K ) THEN
346 K = -IWORK( K )
347 GO TO 100
348 END IF
349 ELSE IF( IWORK( K ).NE.K ) THEN
350 K = IWORK( K )
351 GO TO 100
352 END IF
353 END IF
354 *
355 * Check error code from DSYSV .
356 *
357 IF( INFO.NE.K ) THEN
358 CALL ALAERH( PATH, 'DSYSV ', INFO, K, UPLO, N,
359 $ N, -1, -1, NRHS, IMAT, NFAIL,
360 $ NERRS, NOUT )
361 GO TO 120
362 ELSE IF( INFO.NE.0 ) THEN
363 GO TO 120
364 END IF
365 *
366 * Reconstruct matrix from factors and compute
367 * residual.
368 *
369 CALL DSYT01( UPLO, N, A, LDA, AFAC, LDA, IWORK,
370 $ AINV, LDA, RWORK, RESULT( 1 ) )
371 *
372 * Compute residual of the computed solution.
373 *
374 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
375 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
376 $ LDA, RWORK, RESULT( 2 ) )
377 *
378 * Check solution from generated exact solution.
379 *
380 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
381 $ RESULT( 3 ) )
382 NT = 3
383 *
384 * Print information about the tests that did not pass
385 * the threshold.
386 *
387 DO 110 K = 1, NT
388 IF( RESULT( K ).GE.THRESH ) THEN
389 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
390 $ CALL ALADHD( NOUT, PATH )
391 WRITE( NOUT, FMT = 9999 )'DSYSV ', UPLO, N,
392 $ IMAT, K, RESULT( K )
393 NFAIL = NFAIL + 1
394 END IF
395 110 CONTINUE
396 NRUN = NRUN + NT
397 120 CONTINUE
398 END IF
399 *
400 * --- Test DSYSVX ---
401 *
402 IF( IFACT.EQ.2 )
403 $ CALL DLASET( UPLO, N, N, ZERO, ZERO, AFAC, LDA )
404 CALL DLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDA )
405 *
406 * Solve the system and compute the condition number and
407 * error bounds using DSYSVX.
408 *
409 SRNAMT = 'DSYSVX'
410 CALL DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AFAC, LDA,
411 $ IWORK, B, LDA, X, LDA, RCOND, RWORK,
412 $ RWORK( NRHS+1 ), WORK, LWORK,
413 $ IWORK( N+1 ), INFO )
414 *
415 * Adjust the expected value of INFO to account for
416 * pivoting.
417 *
418 K = IZERO
419 IF( K.GT.0 ) THEN
420 130 CONTINUE
421 IF( IWORK( K ).LT.0 ) THEN
422 IF( IWORK( K ).NE.-K ) THEN
423 K = -IWORK( K )
424 GO TO 130
425 END IF
426 ELSE IF( IWORK( K ).NE.K ) THEN
427 K = IWORK( K )
428 GO TO 130
429 END IF
430 END IF
431 *
432 * Check the error code from DSYSVX.
433 *
434 IF( INFO.NE.K ) THEN
435 CALL ALAERH( PATH, 'DSYSVX', INFO, K, FACT // UPLO,
436 $ N, N, -1, -1, NRHS, IMAT, NFAIL,
437 $ NERRS, NOUT )
438 GO TO 150
439 END IF
440 *
441 IF( INFO.EQ.0 ) THEN
442 IF( IFACT.GE.2 ) THEN
443 *
444 * Reconstruct matrix from factors and compute
445 * residual.
446 *
447 CALL DSYT01( UPLO, N, A, LDA, AFAC, LDA, IWORK,
448 $ AINV, LDA, RWORK( 2*NRHS+1 ),
449 $ RESULT( 1 ) )
450 K1 = 1
451 ELSE
452 K1 = 2
453 END IF
454 *
455 * Compute residual of the computed solution.
456 *
457 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
458 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
459 $ LDA, RWORK( 2*NRHS+1 ), RESULT( 2 ) )
460 *
461 * Check solution from generated exact solution.
462 *
463 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
464 $ RESULT( 3 ) )
465 *
466 * Check the error bounds from iterative refinement.
467 *
468 CALL DPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
469 $ XACT, LDA, RWORK, RWORK( NRHS+1 ),
470 $ RESULT( 4 ) )
471 ELSE
472 K1 = 6
473 END IF
474 *
475 * Compare RCOND from DSYSVX with the computed value
476 * in RCONDC.
477 *
478 RESULT( 6 ) = DGET06( RCOND, RCONDC )
479 *
480 * Print information about the tests that did not pass
481 * the threshold.
482 *
483 DO 140 K = K1, 6
484 IF( RESULT( K ).GE.THRESH ) THEN
485 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
486 $ CALL ALADHD( NOUT, PATH )
487 WRITE( NOUT, FMT = 9998 )'DSYSVX', FACT, UPLO,
488 $ N, IMAT, K, RESULT( K )
489 NFAIL = NFAIL + 1
490 END IF
491 140 CONTINUE
492 NRUN = NRUN + 7 - K1
493 *
494 * --- Test DSYSVXX ---
495 *
496 * Restore the matrices A and B.
497 *
498 IF( IFACT.EQ.2 )
499 $ CALL DLASET( UPLO, N, N, ZERO, ZERO, AFAC, LDA )
500 CALL DLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDA )
501 *
502 * Solve the system and compute the condition number
503 * and error bounds using DSYSVXX.
504 *
505 SRNAMT = 'DSYSVXX'
506 N_ERR_BNDS = 3
507 EQUED = 'N'
508 CALL DSYSVXX( FACT, UPLO, N, NRHS, A, LDA, AFAC,
509 $ LDA, IWORK, EQUED, WORK( N+1 ), B, LDA, X,
510 $ LDA, RCOND, RPVGRW_SVXX, BERR, N_ERR_BNDS,
511 $ ERRBNDS_N, ERRBNDS_C, 0, ZERO, WORK,
512 $ IWORK( N+1 ), INFO )
513 *
514 * Adjust the expected value of INFO to account for
515 * pivoting.
516 *
517 K = IZERO
518 IF( K.GT.0 ) THEN
519 135 CONTINUE
520 IF( IWORK( K ).LT.0 ) THEN
521 IF( IWORK( K ).NE.-K ) THEN
522 K = -IWORK( K )
523 GO TO 135
524 END IF
525 ELSE IF( IWORK( K ).NE.K ) THEN
526 K = IWORK( K )
527 GO TO 135
528 END IF
529 END IF
530 *
531 * Check the error code from DSYSVXX.
532 *
533 IF( INFO.NE.K ) THEN
534 CALL ALAERH( PATH, 'DSYSVXX', INFO, K,
535 $ FACT // UPLO, N, N, -1, -1, NRHS, IMAT, NFAIL,
536 $ NERRS, NOUT )
537 GO TO 150
538 END IF
539 *
540 IF( INFO.EQ.0 ) THEN
541 IF( IFACT.GE.2 ) THEN
542 *
543 * Reconstruct matrix from factors and compute
544 * residual.
545 *
546 CALL DSYT01( UPLO, N, A, LDA, AFAC, LDA, IWORK,
547 $ AINV, LDA, RWORK(2*NRHS+1),
548 $ RESULT( 1 ) )
549 K1 = 1
550 ELSE
551 K1 = 2
552 END IF
553 *
554 * Compute residual of the computed solution.
555 *
556 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
557 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
558 $ LDA, RWORK( 2*NRHS+1 ), RESULT( 2 ) )
559 *
560 * Check solution from generated exact solution.
561 *
562 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
563 $ RESULT( 3 ) )
564 *
565 * Check the error bounds from iterative refinement.
566 *
567 CALL DPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
568 $ XACT, LDA, RWORK, RWORK( NRHS+1 ),
569 $ RESULT( 4 ) )
570 ELSE
571 K1 = 6
572 END IF
573 *
574 * Compare RCOND from DSYSVXX with the computed value
575 * in RCONDC.
576 *
577 RESULT( 6 ) = DGET06( RCOND, RCONDC )
578 *
579 * Print information about the tests that did not pass
580 * the threshold.
581 *
582 DO 85 K = K1, 6
583 IF( RESULT( K ).GE.THRESH ) THEN
584 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
585 $ CALL ALADHD( NOUT, PATH )
586 WRITE( NOUT, FMT = 9998 )'DSYSVXX',
587 $ FACT, UPLO, N, IMAT, K,
588 $ RESULT( K )
589 NFAIL = NFAIL + 1
590 END IF
591 85 CONTINUE
592 NRUN = NRUN + 7 - K1
593 *
594 150 CONTINUE
595 *
596 160 CONTINUE
597 170 CONTINUE
598 180 CONTINUE
599 *
600 * Print a summary of the results.
601 *
602 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
603 *
604
605 * Test Error Bounds from DSYSVXX
606
607 CALL DEBCHVXX(THRESH, PATH)
608
609 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2,
610 $ ', test ', I2, ', ratio =', G12.5 )
611 9998 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N =', I5,
612 $ ', type ', I2, ', test ', I2, ', ratio =', G12.5 )
613 RETURN
614 *
615 * End of DDRVSY
616 *
617 END
2 $ A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK,
3 $ NOUT )
4 *
5 * -- LAPACK test routine (version 3.2) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 LOGICAL TSTERR
11 INTEGER NMAX, NN, NOUT, NRHS
12 DOUBLE PRECISION THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL DOTYPE( * )
16 INTEGER IWORK( * ), NVAL( * )
17 DOUBLE PRECISION A( * ), AFAC( * ), AINV( * ), B( * ),
18 $ RWORK( * ), WORK( * ), X( * ), XACT( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DDRVSY tests the driver routines DSYSV, -SVX, and -SVXX.
25 *
26 * Note that this file is used only when the XBLAS are available,
27 * otherwise ddrvsy.f defines this subroutine.
28 *
29 * Arguments
30 * =========
31 *
32 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
33 * The matrix types to be used for testing. Matrices of type j
34 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
35 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
36 *
37 * NN (input) INTEGER
38 * The number of values of N contained in the vector NVAL.
39 *
40 * NVAL (input) INTEGER array, dimension (NN)
41 * The values of the matrix dimension N.
42 *
43 * NRHS (input) INTEGER
44 * The number of right hand side vectors to be generated for
45 * each linear system.
46 *
47 * THRESH (input) DOUBLE PRECISION
48 * The threshold value for the test ratios. A result is
49 * included in the output file if RESULT >= THRESH. To have
50 * every test ratio printed, use THRESH = 0.
51 *
52 * TSTERR (input) LOGICAL
53 * Flag that indicates whether error exits are to be tested.
54 *
55 * NMAX (input) INTEGER
56 * The maximum value permitted for N, used in dimensioning the
57 * work arrays.
58 *
59 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
60 *
61 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
62 *
63 * AINV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
64 *
65 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS)
66 *
67 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS)
68 *
69 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS)
70 *
71 * WORK (workspace) DOUBLE PRECISION array, dimension
72 * (NMAX*max(2,NRHS))
73 *
74 * RWORK (workspace) DOUBLE PRECISION array, dimension (NMAX+2*NRHS)
75 *
76 * IWORK (workspace) INTEGER array, dimension (2*NMAX)
77 *
78 * NOUT (input) INTEGER
79 * The unit number for output.
80 *
81 * =====================================================================
82 *
83 * .. Parameters ..
84 DOUBLE PRECISION ONE, ZERO
85 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
86 INTEGER NTYPES, NTESTS
87 PARAMETER ( NTYPES = 10, NTESTS = 6 )
88 INTEGER NFACT
89 PARAMETER ( NFACT = 2 )
90 * ..
91 * .. Local Scalars ..
92 LOGICAL ZEROT
93 CHARACTER DIST, EQUED, FACT, TYPE, UPLO, XTYPE
94 CHARACTER*3 PATH
95 INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
96 $ IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N,
97 $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT,
98 $ N_ERR_BNDS
99 DOUBLE PRECISION AINVNM, ANORM, CNDNUM, RCOND, RCONDC,
100 $ RPVGRW_SVXX
101 * ..
102 * .. Local Arrays ..
103 CHARACTER FACTS( NFACT ), UPLOS( 2 )
104 INTEGER ISEED( 4 ), ISEEDY( 4 )
105 DOUBLE PRECISION RESULT( NTESTS ), BERR( NRHS ),
106 $ ERRBNDS_N( NRHS, 3 ), ERRBNDS_C( NRHS, 3 )
107 * ..
108 * .. External Functions ..
109 DOUBLE PRECISION DGET06, DLANSY
110 EXTERNAL DGET06, DLANSY
111 * ..
112 * .. External Subroutines ..
113 EXTERNAL ALADHD, ALAERH, ALASVM, DERRVX, DGET04, DLACPY,
114 $ DLARHS, DLASET, DLATB4, DLATMS, DPOT02, DPOT05,
115 $ DSYSV, DSYSVX, DSYT01, DSYTRF, DSYTRI2, XLAENV,
116 $ DSYSVXX
117 * ..
118 * .. Scalars in Common ..
119 LOGICAL LERR, OK
120 CHARACTER*32 SRNAMT
121 INTEGER INFOT, NUNIT
122 * ..
123 * .. Common blocks ..
124 COMMON / INFOC / INFOT, NUNIT, OK, LERR
125 COMMON / SRNAMC / SRNAMT
126 * ..
127 * .. Intrinsic Functions ..
128 INTRINSIC MAX, MIN
129 * ..
130 * .. Data statements ..
131 DATA ISEEDY / 1988, 1989, 1990, 1991 /
132 DATA UPLOS / 'U', 'L' / , FACTS / 'F', 'N' /
133 * ..
134 * .. Executable Statements ..
135 *
136 * Initialize constants and the random number seed.
137 *
138 PATH( 1: 1 ) = 'Double precision'
139 PATH( 2: 3 ) = 'SY'
140 NRUN = 0
141 NFAIL = 0
142 NERRS = 0
143 DO 10 I = 1, 4
144 ISEED( I ) = ISEEDY( I )
145 10 CONTINUE
146 LWORK = MAX( 2*NMAX, NMAX*NRHS )
147 *
148 * Test the error exits
149 *
150 IF( TSTERR )
151 $ CALL DERRVX( PATH, NOUT )
152 INFOT = 0
153 *
154 * Set the block size and minimum block size for testing.
155 *
156 NB = 1
157 NBMIN = 2
158 CALL XLAENV( 1, NB )
159 CALL XLAENV( 2, NBMIN )
160 *
161 * Do for each value of N in NVAL
162 *
163 DO 180 IN = 1, NN
164 N = NVAL( IN )
165 LDA = MAX( N, 1 )
166 XTYPE = 'N'
167 NIMAT = NTYPES
168 IF( N.LE.0 )
169 $ NIMAT = 1
170 *
171 DO 170 IMAT = 1, NIMAT
172 *
173 * Do the tests only if DOTYPE( IMAT ) is true.
174 *
175 IF( .NOT.DOTYPE( IMAT ) )
176 $ GO TO 170
177 *
178 * Skip types 3, 4, 5, or 6 if the matrix size is too small.
179 *
180 ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
181 IF( ZEROT .AND. N.LT.IMAT-2 )
182 $ GO TO 170
183 *
184 * Do first for UPLO = 'U', then for UPLO = 'L'
185 *
186 DO 160 IUPLO = 1, 2
187 UPLO = UPLOS( IUPLO )
188 *
189 * Set up parameters with DLATB4 and generate a test matrix
190 * with DLATMS.
191 *
192 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
193 $ CNDNUM, DIST )
194 *
195 SRNAMT = 'DLATMS'
196 CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
197 $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
198 $ INFO )
199 *
200 * Check error code from DLATMS.
201 *
202 IF( INFO.NE.0 ) THEN
203 CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1,
204 $ -1, -1, IMAT, NFAIL, NERRS, NOUT )
205 GO TO 160
206 END IF
207 *
208 * For types 3-6, zero one or more rows and columns of the
209 * matrix to test that INFO is returned correctly.
210 *
211 IF( ZEROT ) THEN
212 IF( IMAT.EQ.3 ) THEN
213 IZERO = 1
214 ELSE IF( IMAT.EQ.4 ) THEN
215 IZERO = N
216 ELSE
217 IZERO = N / 2 + 1
218 END IF
219 *
220 IF( IMAT.LT.6 ) THEN
221 *
222 * Set row and column IZERO to zero.
223 *
224 IF( IUPLO.EQ.1 ) THEN
225 IOFF = ( IZERO-1 )*LDA
226 DO 20 I = 1, IZERO - 1
227 A( IOFF+I ) = ZERO
228 20 CONTINUE
229 IOFF = IOFF + IZERO
230 DO 30 I = IZERO, N
231 A( IOFF ) = ZERO
232 IOFF = IOFF + LDA
233 30 CONTINUE
234 ELSE
235 IOFF = IZERO
236 DO 40 I = 1, IZERO - 1
237 A( IOFF ) = ZERO
238 IOFF = IOFF + LDA
239 40 CONTINUE
240 IOFF = IOFF - IZERO
241 DO 50 I = IZERO, N
242 A( IOFF+I ) = ZERO
243 50 CONTINUE
244 END IF
245 ELSE
246 IOFF = 0
247 IF( IUPLO.EQ.1 ) THEN
248 *
249 * Set the first IZERO rows and columns to zero.
250 *
251 DO 70 J = 1, N
252 I2 = MIN( J, IZERO )
253 DO 60 I = 1, I2
254 A( IOFF+I ) = ZERO
255 60 CONTINUE
256 IOFF = IOFF + LDA
257 70 CONTINUE
258 ELSE
259 *
260 * Set the last IZERO rows and columns to zero.
261 *
262 DO 90 J = 1, N
263 I1 = MAX( J, IZERO )
264 DO 80 I = I1, N
265 A( IOFF+I ) = ZERO
266 80 CONTINUE
267 IOFF = IOFF + LDA
268 90 CONTINUE
269 END IF
270 END IF
271 ELSE
272 IZERO = 0
273 END IF
274 *
275 DO 150 IFACT = 1, NFACT
276 *
277 * Do first for FACT = 'F', then for other values.
278 *
279 FACT = FACTS( IFACT )
280 *
281 * Compute the condition number for comparison with
282 * the value returned by DSYSVX.
283 *
284 IF( ZEROT ) THEN
285 IF( IFACT.EQ.1 )
286 $ GO TO 150
287 RCONDC = ZERO
288 *
289 ELSE IF( IFACT.EQ.1 ) THEN
290 *
291 * Compute the 1-norm of A.
292 *
293 ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK )
294 *
295 * Factor the matrix A.
296 *
297 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
298 CALL DSYTRF( UPLO, N, AFAC, LDA, IWORK, WORK,
299 $ LWORK, INFO )
300 *
301 * Compute inv(A) and take its norm.
302 *
303 CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
304 LWORK = (N+NB+1)*(NB+3)
305 CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
306 $ LWORK, INFO )
307 AINVNM = DLANSY( '1', UPLO, N, AINV, LDA, RWORK )
308 *
309 * Compute the 1-norm condition number of A.
310 *
311 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
312 RCONDC = ONE
313 ELSE
314 RCONDC = ( ONE / ANORM ) / AINVNM
315 END IF
316 END IF
317 *
318 * Form an exact solution and set the right hand side.
319 *
320 SRNAMT = 'DLARHS'
321 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
322 $ NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
323 $ INFO )
324 XTYPE = 'C'
325 *
326 * --- Test DSYSV ---
327 *
328 IF( IFACT.EQ.2 ) THEN
329 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
330 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
331 *
332 * Factor the matrix and solve the system using DSYSV.
333 *
334 SRNAMT = 'DSYSV '
335 CALL DSYSV( UPLO, N, NRHS, AFAC, LDA, IWORK, X,
336 $ LDA, WORK, LWORK, INFO )
337 *
338 * Adjust the expected value of INFO to account for
339 * pivoting.
340 *
341 K = IZERO
342 IF( K.GT.0 ) THEN
343 100 CONTINUE
344 IF( IWORK( K ).LT.0 ) THEN
345 IF( IWORK( K ).NE.-K ) THEN
346 K = -IWORK( K )
347 GO TO 100
348 END IF
349 ELSE IF( IWORK( K ).NE.K ) THEN
350 K = IWORK( K )
351 GO TO 100
352 END IF
353 END IF
354 *
355 * Check error code from DSYSV .
356 *
357 IF( INFO.NE.K ) THEN
358 CALL ALAERH( PATH, 'DSYSV ', INFO, K, UPLO, N,
359 $ N, -1, -1, NRHS, IMAT, NFAIL,
360 $ NERRS, NOUT )
361 GO TO 120
362 ELSE IF( INFO.NE.0 ) THEN
363 GO TO 120
364 END IF
365 *
366 * Reconstruct matrix from factors and compute
367 * residual.
368 *
369 CALL DSYT01( UPLO, N, A, LDA, AFAC, LDA, IWORK,
370 $ AINV, LDA, RWORK, RESULT( 1 ) )
371 *
372 * Compute residual of the computed solution.
373 *
374 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
375 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
376 $ LDA, RWORK, RESULT( 2 ) )
377 *
378 * Check solution from generated exact solution.
379 *
380 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
381 $ RESULT( 3 ) )
382 NT = 3
383 *
384 * Print information about the tests that did not pass
385 * the threshold.
386 *
387 DO 110 K = 1, NT
388 IF( RESULT( K ).GE.THRESH ) THEN
389 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
390 $ CALL ALADHD( NOUT, PATH )
391 WRITE( NOUT, FMT = 9999 )'DSYSV ', UPLO, N,
392 $ IMAT, K, RESULT( K )
393 NFAIL = NFAIL + 1
394 END IF
395 110 CONTINUE
396 NRUN = NRUN + NT
397 120 CONTINUE
398 END IF
399 *
400 * --- Test DSYSVX ---
401 *
402 IF( IFACT.EQ.2 )
403 $ CALL DLASET( UPLO, N, N, ZERO, ZERO, AFAC, LDA )
404 CALL DLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDA )
405 *
406 * Solve the system and compute the condition number and
407 * error bounds using DSYSVX.
408 *
409 SRNAMT = 'DSYSVX'
410 CALL DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AFAC, LDA,
411 $ IWORK, B, LDA, X, LDA, RCOND, RWORK,
412 $ RWORK( NRHS+1 ), WORK, LWORK,
413 $ IWORK( N+1 ), INFO )
414 *
415 * Adjust the expected value of INFO to account for
416 * pivoting.
417 *
418 K = IZERO
419 IF( K.GT.0 ) THEN
420 130 CONTINUE
421 IF( IWORK( K ).LT.0 ) THEN
422 IF( IWORK( K ).NE.-K ) THEN
423 K = -IWORK( K )
424 GO TO 130
425 END IF
426 ELSE IF( IWORK( K ).NE.K ) THEN
427 K = IWORK( K )
428 GO TO 130
429 END IF
430 END IF
431 *
432 * Check the error code from DSYSVX.
433 *
434 IF( INFO.NE.K ) THEN
435 CALL ALAERH( PATH, 'DSYSVX', INFO, K, FACT // UPLO,
436 $ N, N, -1, -1, NRHS, IMAT, NFAIL,
437 $ NERRS, NOUT )
438 GO TO 150
439 END IF
440 *
441 IF( INFO.EQ.0 ) THEN
442 IF( IFACT.GE.2 ) THEN
443 *
444 * Reconstruct matrix from factors and compute
445 * residual.
446 *
447 CALL DSYT01( UPLO, N, A, LDA, AFAC, LDA, IWORK,
448 $ AINV, LDA, RWORK( 2*NRHS+1 ),
449 $ RESULT( 1 ) )
450 K1 = 1
451 ELSE
452 K1 = 2
453 END IF
454 *
455 * Compute residual of the computed solution.
456 *
457 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
458 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
459 $ LDA, RWORK( 2*NRHS+1 ), RESULT( 2 ) )
460 *
461 * Check solution from generated exact solution.
462 *
463 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
464 $ RESULT( 3 ) )
465 *
466 * Check the error bounds from iterative refinement.
467 *
468 CALL DPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
469 $ XACT, LDA, RWORK, RWORK( NRHS+1 ),
470 $ RESULT( 4 ) )
471 ELSE
472 K1 = 6
473 END IF
474 *
475 * Compare RCOND from DSYSVX with the computed value
476 * in RCONDC.
477 *
478 RESULT( 6 ) = DGET06( RCOND, RCONDC )
479 *
480 * Print information about the tests that did not pass
481 * the threshold.
482 *
483 DO 140 K = K1, 6
484 IF( RESULT( K ).GE.THRESH ) THEN
485 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
486 $ CALL ALADHD( NOUT, PATH )
487 WRITE( NOUT, FMT = 9998 )'DSYSVX', FACT, UPLO,
488 $ N, IMAT, K, RESULT( K )
489 NFAIL = NFAIL + 1
490 END IF
491 140 CONTINUE
492 NRUN = NRUN + 7 - K1
493 *
494 * --- Test DSYSVXX ---
495 *
496 * Restore the matrices A and B.
497 *
498 IF( IFACT.EQ.2 )
499 $ CALL DLASET( UPLO, N, N, ZERO, ZERO, AFAC, LDA )
500 CALL DLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDA )
501 *
502 * Solve the system and compute the condition number
503 * and error bounds using DSYSVXX.
504 *
505 SRNAMT = 'DSYSVXX'
506 N_ERR_BNDS = 3
507 EQUED = 'N'
508 CALL DSYSVXX( FACT, UPLO, N, NRHS, A, LDA, AFAC,
509 $ LDA, IWORK, EQUED, WORK( N+1 ), B, LDA, X,
510 $ LDA, RCOND, RPVGRW_SVXX, BERR, N_ERR_BNDS,
511 $ ERRBNDS_N, ERRBNDS_C, 0, ZERO, WORK,
512 $ IWORK( N+1 ), INFO )
513 *
514 * Adjust the expected value of INFO to account for
515 * pivoting.
516 *
517 K = IZERO
518 IF( K.GT.0 ) THEN
519 135 CONTINUE
520 IF( IWORK( K ).LT.0 ) THEN
521 IF( IWORK( K ).NE.-K ) THEN
522 K = -IWORK( K )
523 GO TO 135
524 END IF
525 ELSE IF( IWORK( K ).NE.K ) THEN
526 K = IWORK( K )
527 GO TO 135
528 END IF
529 END IF
530 *
531 * Check the error code from DSYSVXX.
532 *
533 IF( INFO.NE.K ) THEN
534 CALL ALAERH( PATH, 'DSYSVXX', INFO, K,
535 $ FACT // UPLO, N, N, -1, -1, NRHS, IMAT, NFAIL,
536 $ NERRS, NOUT )
537 GO TO 150
538 END IF
539 *
540 IF( INFO.EQ.0 ) THEN
541 IF( IFACT.GE.2 ) THEN
542 *
543 * Reconstruct matrix from factors and compute
544 * residual.
545 *
546 CALL DSYT01( UPLO, N, A, LDA, AFAC, LDA, IWORK,
547 $ AINV, LDA, RWORK(2*NRHS+1),
548 $ RESULT( 1 ) )
549 K1 = 1
550 ELSE
551 K1 = 2
552 END IF
553 *
554 * Compute residual of the computed solution.
555 *
556 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
557 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
558 $ LDA, RWORK( 2*NRHS+1 ), RESULT( 2 ) )
559 *
560 * Check solution from generated exact solution.
561 *
562 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
563 $ RESULT( 3 ) )
564 *
565 * Check the error bounds from iterative refinement.
566 *
567 CALL DPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
568 $ XACT, LDA, RWORK, RWORK( NRHS+1 ),
569 $ RESULT( 4 ) )
570 ELSE
571 K1 = 6
572 END IF
573 *
574 * Compare RCOND from DSYSVXX with the computed value
575 * in RCONDC.
576 *
577 RESULT( 6 ) = DGET06( RCOND, RCONDC )
578 *
579 * Print information about the tests that did not pass
580 * the threshold.
581 *
582 DO 85 K = K1, 6
583 IF( RESULT( K ).GE.THRESH ) THEN
584 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
585 $ CALL ALADHD( NOUT, PATH )
586 WRITE( NOUT, FMT = 9998 )'DSYSVXX',
587 $ FACT, UPLO, N, IMAT, K,
588 $ RESULT( K )
589 NFAIL = NFAIL + 1
590 END IF
591 85 CONTINUE
592 NRUN = NRUN + 7 - K1
593 *
594 150 CONTINUE
595 *
596 160 CONTINUE
597 170 CONTINUE
598 180 CONTINUE
599 *
600 * Print a summary of the results.
601 *
602 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
603 *
604
605 * Test Error Bounds from DSYSVXX
606
607 CALL DEBCHVXX(THRESH, PATH)
608
609 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2,
610 $ ', test ', I2, ', ratio =', G12.5 )
611 9998 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N =', I5,
612 $ ', type ', I2, ', test ', I2, ', ratio =', G12.5 )
613 RETURN
614 *
615 * End of DDRVSY
616 *
617 END