1 SUBROUTINE DDRVPO( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
2 $ A, AFAC, ASAV, B, BSAV, X, XACT, S, WORK,
3 $ RWORK, IWORK, NOUT )
4 *
5 * -- LAPACK test routine (version 3.1) --
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( * ), ASAV( * ), B( * ),
18 $ BSAV( * ), RWORK( * ), S( * ), WORK( * ),
19 $ X( * ), XACT( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DDRVPO tests the driver routines DPOSV and -SVX.
26 *
27 * Arguments
28 * =========
29 *
30 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
31 * The matrix types to be used for testing. Matrices of type j
32 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
33 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
34 *
35 * NN (input) INTEGER
36 * The number of values of N contained in the vector NVAL.
37 *
38 * NVAL (input) INTEGER array, dimension (NN)
39 * The values of the matrix dimension N.
40 *
41 * NRHS (input) INTEGER
42 * The number of right hand side vectors to be generated for
43 * each linear system.
44 *
45 * THRESH (input) DOUBLE PRECISION
46 * The threshold value for the test ratios. A result is
47 * included in the output file if RESULT >= THRESH. To have
48 * every test ratio printed, use THRESH = 0.
49 *
50 * TSTERR (input) LOGICAL
51 * Flag that indicates whether error exits are to be tested.
52 *
53 * NMAX (input) INTEGER
54 * The maximum value permitted for N, used in dimensioning the
55 * work arrays.
56 *
57 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
58 *
59 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
60 *
61 * ASAV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
62 *
63 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS)
64 *
65 * BSAV (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 * S (workspace) DOUBLE PRECISION array, dimension (NMAX)
72 *
73 * WORK (workspace) DOUBLE PRECISION array, dimension
74 * (NMAX*max(3,NRHS))
75 *
76 * RWORK (workspace) DOUBLE PRECISION array, dimension (NMAX+2*NRHS)
77 *
78 * IWORK (workspace) INTEGER array, dimension (NMAX)
79 *
80 * NOUT (input) INTEGER
81 * The unit number for output.
82 *
83 * =====================================================================
84 *
85 * .. Parameters ..
86 DOUBLE PRECISION ONE, ZERO
87 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
88 INTEGER NTYPES
89 PARAMETER ( NTYPES = 9 )
90 INTEGER NTESTS
91 PARAMETER ( NTESTS = 6 )
92 * ..
93 * .. Local Scalars ..
94 LOGICAL EQUIL, NOFACT, PREFAC, ZEROT
95 CHARACTER DIST, EQUED, FACT, TYPE, UPLO, XTYPE
96 CHARACTER*3 PATH
97 INTEGER I, IEQUED, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
98 $ IZERO, K, K1, KL, KU, LDA, MODE, N, NB, NBMIN,
99 $ NERRS, NFACT, NFAIL, NIMAT, NRUN, NT
100 DOUBLE PRECISION AINVNM, AMAX, ANORM, CNDNUM, RCOND, RCONDC,
101 $ ROLDC, SCOND
102 * ..
103 * .. Local Arrays ..
104 CHARACTER EQUEDS( 2 ), FACTS( 3 ), UPLOS( 2 )
105 INTEGER ISEED( 4 ), ISEEDY( 4 )
106 DOUBLE PRECISION RESULT( NTESTS )
107 * ..
108 * .. External Functions ..
109 LOGICAL LSAME
110 DOUBLE PRECISION DGET06, DLANSY
111 EXTERNAL LSAME, DGET06, DLANSY
112 * ..
113 * .. External Subroutines ..
114 EXTERNAL ALADHD, ALAERH, ALASVM, DERRVX, DGET04, DLACPY,
115 $ DLAQSY, DLARHS, DLASET, DLATB4, DLATMS, DPOEQU,
116 $ DPOSV, DPOSVX, DPOT01, DPOT02, DPOT05, DPOTRF,
117 $ DPOTRI, XLAENV
118 * ..
119 * .. Intrinsic Functions ..
120 INTRINSIC MAX
121 * ..
122 * .. Scalars in Common ..
123 LOGICAL LERR, OK
124 CHARACTER*32 SRNAMT
125 INTEGER INFOT, NUNIT
126 * ..
127 * .. Common blocks ..
128 COMMON / INFOC / INFOT, NUNIT, OK, LERR
129 COMMON / SRNAMC / SRNAMT
130 * ..
131 * .. Data statements ..
132 DATA ISEEDY / 1988, 1989, 1990, 1991 /
133 DATA UPLOS / 'U', 'L' /
134 DATA FACTS / 'F', 'N', 'E' /
135 DATA EQUEDS / 'N', 'Y' /
136 * ..
137 * .. Executable Statements ..
138 *
139 * Initialize constants and the random number seed.
140 *
141 PATH( 1: 1 ) = 'Double precision'
142 PATH( 2: 3 ) = 'PO'
143 NRUN = 0
144 NFAIL = 0
145 NERRS = 0
146 DO 10 I = 1, 4
147 ISEED( I ) = ISEEDY( I )
148 10 CONTINUE
149 *
150 * Test the error exits
151 *
152 IF( TSTERR )
153 $ CALL DERRVX( PATH, NOUT )
154 INFOT = 0
155 *
156 * Set the block size and minimum block size for testing.
157 *
158 NB = 1
159 NBMIN = 2
160 CALL XLAENV( 1, NB )
161 CALL XLAENV( 2, NBMIN )
162 *
163 * Do for each value of N in NVAL
164 *
165 DO 130 IN = 1, NN
166 N = NVAL( IN )
167 LDA = MAX( N, 1 )
168 XTYPE = 'N'
169 NIMAT = NTYPES
170 IF( N.LE.0 )
171 $ NIMAT = 1
172 *
173 DO 120 IMAT = 1, NIMAT
174 *
175 * Do the tests only if DOTYPE( IMAT ) is true.
176 *
177 IF( .NOT.DOTYPE( IMAT ) )
178 $ GO TO 120
179 *
180 * Skip types 3, 4, or 5 if the matrix size is too small.
181 *
182 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5
183 IF( ZEROT .AND. N.LT.IMAT-2 )
184 $ GO TO 120
185 *
186 * Do first for UPLO = 'U', then for UPLO = 'L'
187 *
188 DO 110 IUPLO = 1, 2
189 UPLO = UPLOS( IUPLO )
190 *
191 * Set up parameters with DLATB4 and generate a test matrix
192 * with DLATMS.
193 *
194 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
195 $ CNDNUM, DIST )
196 *
197 SRNAMT = 'DLATMS'
198 CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
199 $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
200 $ INFO )
201 *
202 * Check error code from DLATMS.
203 *
204 IF( INFO.NE.0 ) THEN
205 CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1,
206 $ -1, -1, IMAT, NFAIL, NERRS, NOUT )
207 GO TO 110
208 END IF
209 *
210 * For types 3-5, zero one row and column of the matrix to
211 * test that INFO is returned correctly.
212 *
213 IF( ZEROT ) THEN
214 IF( IMAT.EQ.3 ) THEN
215 IZERO = 1
216 ELSE IF( IMAT.EQ.4 ) THEN
217 IZERO = N
218 ELSE
219 IZERO = N / 2 + 1
220 END IF
221 IOFF = ( IZERO-1 )*LDA
222 *
223 * Set row and column IZERO of A to 0.
224 *
225 IF( IUPLO.EQ.1 ) THEN
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 IZERO = 0
247 END IF
248 *
249 * Save a copy of the matrix A in ASAV.
250 *
251 CALL DLACPY( UPLO, N, N, A, LDA, ASAV, LDA )
252 *
253 DO 100 IEQUED = 1, 2
254 EQUED = EQUEDS( IEQUED )
255 IF( IEQUED.EQ.1 ) THEN
256 NFACT = 3
257 ELSE
258 NFACT = 1
259 END IF
260 *
261 DO 90 IFACT = 1, NFACT
262 FACT = FACTS( IFACT )
263 PREFAC = LSAME( FACT, 'F' )
264 NOFACT = LSAME( FACT, 'N' )
265 EQUIL = LSAME( FACT, 'E' )
266 *
267 IF( ZEROT ) THEN
268 IF( PREFAC )
269 $ GO TO 90
270 RCONDC = ZERO
271 *
272 ELSE IF( .NOT.LSAME( FACT, 'N' ) ) THEN
273 *
274 * Compute the condition number for comparison with
275 * the value returned by DPOSVX (FACT = 'N' reuses
276 * the condition number from the previous iteration
277 * with FACT = 'F').
278 *
279 CALL DLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA )
280 IF( EQUIL .OR. IEQUED.GT.1 ) THEN
281 *
282 * Compute row and column scale factors to
283 * equilibrate the matrix A.
284 *
285 CALL DPOEQU( N, AFAC, LDA, S, SCOND, AMAX,
286 $ INFO )
287 IF( INFO.EQ.0 .AND. N.GT.0 ) THEN
288 IF( IEQUED.GT.1 )
289 $ SCOND = ZERO
290 *
291 * Equilibrate the matrix.
292 *
293 CALL DLAQSY( UPLO, N, AFAC, LDA, S, SCOND,
294 $ AMAX, EQUED )
295 END IF
296 END IF
297 *
298 * Save the condition number of the
299 * non-equilibrated system for use in DGET04.
300 *
301 IF( EQUIL )
302 $ ROLDC = RCONDC
303 *
304 * Compute the 1-norm of A.
305 *
306 ANORM = DLANSY( '1', UPLO, N, AFAC, LDA, RWORK )
307 *
308 * Factor the matrix A.
309 *
310 CALL DPOTRF( UPLO, N, AFAC, LDA, INFO )
311 *
312 * Form the inverse of A.
313 *
314 CALL DLACPY( UPLO, N, N, AFAC, LDA, A, LDA )
315 CALL DPOTRI( UPLO, N, A, LDA, INFO )
316 *
317 * Compute the 1-norm condition number of A.
318 *
319 AINVNM = DLANSY( '1', UPLO, N, A, LDA, RWORK )
320 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
321 RCONDC = ONE
322 ELSE
323 RCONDC = ( ONE / ANORM ) / AINVNM
324 END IF
325 END IF
326 *
327 * Restore the matrix A.
328 *
329 CALL DLACPY( UPLO, N, N, ASAV, LDA, A, LDA )
330 *
331 * Form an exact solution and set the right hand side.
332 *
333 SRNAMT = 'DLARHS'
334 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
335 $ NRHS, A, LDA, XACT, LDA, B, LDA,
336 $ ISEED, INFO )
337 XTYPE = 'C'
338 CALL DLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA )
339 *
340 IF( NOFACT ) THEN
341 *
342 * --- Test DPOSV ---
343 *
344 * Compute the L*L' or U'*U factorization of the
345 * matrix and solve the system.
346 *
347 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
348 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
349 *
350 SRNAMT = 'DPOSV '
351 CALL DPOSV( UPLO, N, NRHS, AFAC, LDA, X, LDA,
352 $ INFO )
353 *
354 * Check error code from DPOSV .
355 *
356 IF( INFO.NE.IZERO ) THEN
357 CALL ALAERH( PATH, 'DPOSV ', INFO, IZERO,
358 $ UPLO, N, N, -1, -1, NRHS, IMAT,
359 $ NFAIL, NERRS, NOUT )
360 GO TO 70
361 ELSE IF( INFO.NE.0 ) THEN
362 GO TO 70
363 END IF
364 *
365 * Reconstruct matrix from factors and compute
366 * residual.
367 *
368 CALL DPOT01( UPLO, N, A, LDA, AFAC, LDA, RWORK,
369 $ RESULT( 1 ) )
370 *
371 * Compute residual of the computed solution.
372 *
373 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK,
374 $ LDA )
375 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA,
376 $ WORK, 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
385 * pass the threshold.
386 *
387 DO 60 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 )'DPOSV ', UPLO,
392 $ N, IMAT, K, RESULT( K )
393 NFAIL = NFAIL + 1
394 END IF
395 60 CONTINUE
396 NRUN = NRUN + NT
397 70 CONTINUE
398 END IF
399 *
400 * --- Test DPOSVX ---
401 *
402 IF( .NOT.PREFAC )
403 $ CALL DLASET( UPLO, N, N, ZERO, ZERO, AFAC, LDA )
404 CALL DLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDA )
405 IF( IEQUED.GT.1 .AND. N.GT.0 ) THEN
406 *
407 * Equilibrate the matrix if FACT='F' and
408 * EQUED='Y'.
409 *
410 CALL DLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX,
411 $ EQUED )
412 END IF
413 *
414 * Solve the system and compute the condition number
415 * and error bounds using DPOSVX.
416 *
417 SRNAMT = 'DPOSVX'
418 CALL DPOSVX( FACT, UPLO, N, NRHS, A, LDA, AFAC,
419 $ LDA, EQUED, S, B, LDA, X, LDA, RCOND,
420 $ RWORK, RWORK( NRHS+1 ), WORK, IWORK,
421 $ INFO )
422 *
423 * Check the error code from DPOSVX.
424 *
425 IF( INFO.NE.IZERO ) THEN
426 CALL ALAERH( PATH, 'DPOSVX', INFO, IZERO,
427 $ FACT // UPLO, N, N, -1, -1, NRHS,
428 $ IMAT, NFAIL, NERRS, NOUT )
429 GO TO 90
430 END IF
431 *
432 IF( INFO.EQ.0 ) THEN
433 IF( .NOT.PREFAC ) THEN
434 *
435 * Reconstruct matrix from factors and compute
436 * residual.
437 *
438 CALL DPOT01( UPLO, N, A, LDA, AFAC, LDA,
439 $ RWORK( 2*NRHS+1 ), RESULT( 1 ) )
440 K1 = 1
441 ELSE
442 K1 = 2
443 END IF
444 *
445 * Compute residual of the computed solution.
446 *
447 CALL DLACPY( 'Full', N, NRHS, BSAV, LDA, WORK,
448 $ LDA )
449 CALL DPOT02( UPLO, N, NRHS, ASAV, LDA, X, LDA,
450 $ WORK, LDA, RWORK( 2*NRHS+1 ),
451 $ RESULT( 2 ) )
452 *
453 * Check solution from generated exact solution.
454 *
455 IF( NOFACT .OR. ( PREFAC .AND. LSAME( EQUED,
456 $ 'N' ) ) ) THEN
457 CALL DGET04( N, NRHS, X, LDA, XACT, LDA,
458 $ RCONDC, RESULT( 3 ) )
459 ELSE
460 CALL DGET04( N, NRHS, X, LDA, XACT, LDA,
461 $ ROLDC, RESULT( 3 ) )
462 END IF
463 *
464 * Check the error bounds from iterative
465 * refinement.
466 *
467 CALL DPOT05( UPLO, N, NRHS, ASAV, LDA, B, LDA,
468 $ X, LDA, XACT, LDA, RWORK,
469 $ RWORK( NRHS+1 ), RESULT( 4 ) )
470 ELSE
471 K1 = 6
472 END IF
473 *
474 * Compare RCOND from DPOSVX with the computed value
475 * in RCONDC.
476 *
477 RESULT( 6 ) = DGET06( RCOND, RCONDC )
478 *
479 * Print information about the tests that did not pass
480 * the threshold.
481 *
482 DO 80 K = K1, 6
483 IF( RESULT( K ).GE.THRESH ) THEN
484 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
485 $ CALL ALADHD( NOUT, PATH )
486 IF( PREFAC ) THEN
487 WRITE( NOUT, FMT = 9997 )'DPOSVX', FACT,
488 $ UPLO, N, EQUED, IMAT, K, RESULT( K )
489 ELSE
490 WRITE( NOUT, FMT = 9998 )'DPOSVX', FACT,
491 $ UPLO, N, IMAT, K, RESULT( K )
492 END IF
493 NFAIL = NFAIL + 1
494 END IF
495 80 CONTINUE
496 NRUN = NRUN + 7 - K1
497 90 CONTINUE
498 100 CONTINUE
499 110 CONTINUE
500 120 CONTINUE
501 130 CONTINUE
502 *
503 * Print a summary of the results.
504 *
505 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
506 *
507 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I1,
508 $ ', test(', I1, ')=', G12.5 )
509 9998 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5,
510 $ ', type ', I1, ', test(', I1, ')=', G12.5 )
511 9997 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5,
512 $ ', EQUED=''', A1, ''', type ', I1, ', test(', I1, ') =',
513 $ G12.5 )
514 RETURN
515 *
516 * End of DDRVPO
517 *
518 END
2 $ A, AFAC, ASAV, B, BSAV, X, XACT, S, WORK,
3 $ RWORK, IWORK, NOUT )
4 *
5 * -- LAPACK test routine (version 3.1) --
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( * ), ASAV( * ), B( * ),
18 $ BSAV( * ), RWORK( * ), S( * ), WORK( * ),
19 $ X( * ), XACT( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DDRVPO tests the driver routines DPOSV and -SVX.
26 *
27 * Arguments
28 * =========
29 *
30 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
31 * The matrix types to be used for testing. Matrices of type j
32 * (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
33 * .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
34 *
35 * NN (input) INTEGER
36 * The number of values of N contained in the vector NVAL.
37 *
38 * NVAL (input) INTEGER array, dimension (NN)
39 * The values of the matrix dimension N.
40 *
41 * NRHS (input) INTEGER
42 * The number of right hand side vectors to be generated for
43 * each linear system.
44 *
45 * THRESH (input) DOUBLE PRECISION
46 * The threshold value for the test ratios. A result is
47 * included in the output file if RESULT >= THRESH. To have
48 * every test ratio printed, use THRESH = 0.
49 *
50 * TSTERR (input) LOGICAL
51 * Flag that indicates whether error exits are to be tested.
52 *
53 * NMAX (input) INTEGER
54 * The maximum value permitted for N, used in dimensioning the
55 * work arrays.
56 *
57 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
58 *
59 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
60 *
61 * ASAV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
62 *
63 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NRHS)
64 *
65 * BSAV (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 * S (workspace) DOUBLE PRECISION array, dimension (NMAX)
72 *
73 * WORK (workspace) DOUBLE PRECISION array, dimension
74 * (NMAX*max(3,NRHS))
75 *
76 * RWORK (workspace) DOUBLE PRECISION array, dimension (NMAX+2*NRHS)
77 *
78 * IWORK (workspace) INTEGER array, dimension (NMAX)
79 *
80 * NOUT (input) INTEGER
81 * The unit number for output.
82 *
83 * =====================================================================
84 *
85 * .. Parameters ..
86 DOUBLE PRECISION ONE, ZERO
87 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
88 INTEGER NTYPES
89 PARAMETER ( NTYPES = 9 )
90 INTEGER NTESTS
91 PARAMETER ( NTESTS = 6 )
92 * ..
93 * .. Local Scalars ..
94 LOGICAL EQUIL, NOFACT, PREFAC, ZEROT
95 CHARACTER DIST, EQUED, FACT, TYPE, UPLO, XTYPE
96 CHARACTER*3 PATH
97 INTEGER I, IEQUED, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
98 $ IZERO, K, K1, KL, KU, LDA, MODE, N, NB, NBMIN,
99 $ NERRS, NFACT, NFAIL, NIMAT, NRUN, NT
100 DOUBLE PRECISION AINVNM, AMAX, ANORM, CNDNUM, RCOND, RCONDC,
101 $ ROLDC, SCOND
102 * ..
103 * .. Local Arrays ..
104 CHARACTER EQUEDS( 2 ), FACTS( 3 ), UPLOS( 2 )
105 INTEGER ISEED( 4 ), ISEEDY( 4 )
106 DOUBLE PRECISION RESULT( NTESTS )
107 * ..
108 * .. External Functions ..
109 LOGICAL LSAME
110 DOUBLE PRECISION DGET06, DLANSY
111 EXTERNAL LSAME, DGET06, DLANSY
112 * ..
113 * .. External Subroutines ..
114 EXTERNAL ALADHD, ALAERH, ALASVM, DERRVX, DGET04, DLACPY,
115 $ DLAQSY, DLARHS, DLASET, DLATB4, DLATMS, DPOEQU,
116 $ DPOSV, DPOSVX, DPOT01, DPOT02, DPOT05, DPOTRF,
117 $ DPOTRI, XLAENV
118 * ..
119 * .. Intrinsic Functions ..
120 INTRINSIC MAX
121 * ..
122 * .. Scalars in Common ..
123 LOGICAL LERR, OK
124 CHARACTER*32 SRNAMT
125 INTEGER INFOT, NUNIT
126 * ..
127 * .. Common blocks ..
128 COMMON / INFOC / INFOT, NUNIT, OK, LERR
129 COMMON / SRNAMC / SRNAMT
130 * ..
131 * .. Data statements ..
132 DATA ISEEDY / 1988, 1989, 1990, 1991 /
133 DATA UPLOS / 'U', 'L' /
134 DATA FACTS / 'F', 'N', 'E' /
135 DATA EQUEDS / 'N', 'Y' /
136 * ..
137 * .. Executable Statements ..
138 *
139 * Initialize constants and the random number seed.
140 *
141 PATH( 1: 1 ) = 'Double precision'
142 PATH( 2: 3 ) = 'PO'
143 NRUN = 0
144 NFAIL = 0
145 NERRS = 0
146 DO 10 I = 1, 4
147 ISEED( I ) = ISEEDY( I )
148 10 CONTINUE
149 *
150 * Test the error exits
151 *
152 IF( TSTERR )
153 $ CALL DERRVX( PATH, NOUT )
154 INFOT = 0
155 *
156 * Set the block size and minimum block size for testing.
157 *
158 NB = 1
159 NBMIN = 2
160 CALL XLAENV( 1, NB )
161 CALL XLAENV( 2, NBMIN )
162 *
163 * Do for each value of N in NVAL
164 *
165 DO 130 IN = 1, NN
166 N = NVAL( IN )
167 LDA = MAX( N, 1 )
168 XTYPE = 'N'
169 NIMAT = NTYPES
170 IF( N.LE.0 )
171 $ NIMAT = 1
172 *
173 DO 120 IMAT = 1, NIMAT
174 *
175 * Do the tests only if DOTYPE( IMAT ) is true.
176 *
177 IF( .NOT.DOTYPE( IMAT ) )
178 $ GO TO 120
179 *
180 * Skip types 3, 4, or 5 if the matrix size is too small.
181 *
182 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5
183 IF( ZEROT .AND. N.LT.IMAT-2 )
184 $ GO TO 120
185 *
186 * Do first for UPLO = 'U', then for UPLO = 'L'
187 *
188 DO 110 IUPLO = 1, 2
189 UPLO = UPLOS( IUPLO )
190 *
191 * Set up parameters with DLATB4 and generate a test matrix
192 * with DLATMS.
193 *
194 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
195 $ CNDNUM, DIST )
196 *
197 SRNAMT = 'DLATMS'
198 CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
199 $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
200 $ INFO )
201 *
202 * Check error code from DLATMS.
203 *
204 IF( INFO.NE.0 ) THEN
205 CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1,
206 $ -1, -1, IMAT, NFAIL, NERRS, NOUT )
207 GO TO 110
208 END IF
209 *
210 * For types 3-5, zero one row and column of the matrix to
211 * test that INFO is returned correctly.
212 *
213 IF( ZEROT ) THEN
214 IF( IMAT.EQ.3 ) THEN
215 IZERO = 1
216 ELSE IF( IMAT.EQ.4 ) THEN
217 IZERO = N
218 ELSE
219 IZERO = N / 2 + 1
220 END IF
221 IOFF = ( IZERO-1 )*LDA
222 *
223 * Set row and column IZERO of A to 0.
224 *
225 IF( IUPLO.EQ.1 ) THEN
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 IZERO = 0
247 END IF
248 *
249 * Save a copy of the matrix A in ASAV.
250 *
251 CALL DLACPY( UPLO, N, N, A, LDA, ASAV, LDA )
252 *
253 DO 100 IEQUED = 1, 2
254 EQUED = EQUEDS( IEQUED )
255 IF( IEQUED.EQ.1 ) THEN
256 NFACT = 3
257 ELSE
258 NFACT = 1
259 END IF
260 *
261 DO 90 IFACT = 1, NFACT
262 FACT = FACTS( IFACT )
263 PREFAC = LSAME( FACT, 'F' )
264 NOFACT = LSAME( FACT, 'N' )
265 EQUIL = LSAME( FACT, 'E' )
266 *
267 IF( ZEROT ) THEN
268 IF( PREFAC )
269 $ GO TO 90
270 RCONDC = ZERO
271 *
272 ELSE IF( .NOT.LSAME( FACT, 'N' ) ) THEN
273 *
274 * Compute the condition number for comparison with
275 * the value returned by DPOSVX (FACT = 'N' reuses
276 * the condition number from the previous iteration
277 * with FACT = 'F').
278 *
279 CALL DLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA )
280 IF( EQUIL .OR. IEQUED.GT.1 ) THEN
281 *
282 * Compute row and column scale factors to
283 * equilibrate the matrix A.
284 *
285 CALL DPOEQU( N, AFAC, LDA, S, SCOND, AMAX,
286 $ INFO )
287 IF( INFO.EQ.0 .AND. N.GT.0 ) THEN
288 IF( IEQUED.GT.1 )
289 $ SCOND = ZERO
290 *
291 * Equilibrate the matrix.
292 *
293 CALL DLAQSY( UPLO, N, AFAC, LDA, S, SCOND,
294 $ AMAX, EQUED )
295 END IF
296 END IF
297 *
298 * Save the condition number of the
299 * non-equilibrated system for use in DGET04.
300 *
301 IF( EQUIL )
302 $ ROLDC = RCONDC
303 *
304 * Compute the 1-norm of A.
305 *
306 ANORM = DLANSY( '1', UPLO, N, AFAC, LDA, RWORK )
307 *
308 * Factor the matrix A.
309 *
310 CALL DPOTRF( UPLO, N, AFAC, LDA, INFO )
311 *
312 * Form the inverse of A.
313 *
314 CALL DLACPY( UPLO, N, N, AFAC, LDA, A, LDA )
315 CALL DPOTRI( UPLO, N, A, LDA, INFO )
316 *
317 * Compute the 1-norm condition number of A.
318 *
319 AINVNM = DLANSY( '1', UPLO, N, A, LDA, RWORK )
320 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
321 RCONDC = ONE
322 ELSE
323 RCONDC = ( ONE / ANORM ) / AINVNM
324 END IF
325 END IF
326 *
327 * Restore the matrix A.
328 *
329 CALL DLACPY( UPLO, N, N, ASAV, LDA, A, LDA )
330 *
331 * Form an exact solution and set the right hand side.
332 *
333 SRNAMT = 'DLARHS'
334 CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
335 $ NRHS, A, LDA, XACT, LDA, B, LDA,
336 $ ISEED, INFO )
337 XTYPE = 'C'
338 CALL DLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA )
339 *
340 IF( NOFACT ) THEN
341 *
342 * --- Test DPOSV ---
343 *
344 * Compute the L*L' or U'*U factorization of the
345 * matrix and solve the system.
346 *
347 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
348 CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
349 *
350 SRNAMT = 'DPOSV '
351 CALL DPOSV( UPLO, N, NRHS, AFAC, LDA, X, LDA,
352 $ INFO )
353 *
354 * Check error code from DPOSV .
355 *
356 IF( INFO.NE.IZERO ) THEN
357 CALL ALAERH( PATH, 'DPOSV ', INFO, IZERO,
358 $ UPLO, N, N, -1, -1, NRHS, IMAT,
359 $ NFAIL, NERRS, NOUT )
360 GO TO 70
361 ELSE IF( INFO.NE.0 ) THEN
362 GO TO 70
363 END IF
364 *
365 * Reconstruct matrix from factors and compute
366 * residual.
367 *
368 CALL DPOT01( UPLO, N, A, LDA, AFAC, LDA, RWORK,
369 $ RESULT( 1 ) )
370 *
371 * Compute residual of the computed solution.
372 *
373 CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK,
374 $ LDA )
375 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA,
376 $ WORK, 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
385 * pass the threshold.
386 *
387 DO 60 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 )'DPOSV ', UPLO,
392 $ N, IMAT, K, RESULT( K )
393 NFAIL = NFAIL + 1
394 END IF
395 60 CONTINUE
396 NRUN = NRUN + NT
397 70 CONTINUE
398 END IF
399 *
400 * --- Test DPOSVX ---
401 *
402 IF( .NOT.PREFAC )
403 $ CALL DLASET( UPLO, N, N, ZERO, ZERO, AFAC, LDA )
404 CALL DLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDA )
405 IF( IEQUED.GT.1 .AND. N.GT.0 ) THEN
406 *
407 * Equilibrate the matrix if FACT='F' and
408 * EQUED='Y'.
409 *
410 CALL DLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX,
411 $ EQUED )
412 END IF
413 *
414 * Solve the system and compute the condition number
415 * and error bounds using DPOSVX.
416 *
417 SRNAMT = 'DPOSVX'
418 CALL DPOSVX( FACT, UPLO, N, NRHS, A, LDA, AFAC,
419 $ LDA, EQUED, S, B, LDA, X, LDA, RCOND,
420 $ RWORK, RWORK( NRHS+1 ), WORK, IWORK,
421 $ INFO )
422 *
423 * Check the error code from DPOSVX.
424 *
425 IF( INFO.NE.IZERO ) THEN
426 CALL ALAERH( PATH, 'DPOSVX', INFO, IZERO,
427 $ FACT // UPLO, N, N, -1, -1, NRHS,
428 $ IMAT, NFAIL, NERRS, NOUT )
429 GO TO 90
430 END IF
431 *
432 IF( INFO.EQ.0 ) THEN
433 IF( .NOT.PREFAC ) THEN
434 *
435 * Reconstruct matrix from factors and compute
436 * residual.
437 *
438 CALL DPOT01( UPLO, N, A, LDA, AFAC, LDA,
439 $ RWORK( 2*NRHS+1 ), RESULT( 1 ) )
440 K1 = 1
441 ELSE
442 K1 = 2
443 END IF
444 *
445 * Compute residual of the computed solution.
446 *
447 CALL DLACPY( 'Full', N, NRHS, BSAV, LDA, WORK,
448 $ LDA )
449 CALL DPOT02( UPLO, N, NRHS, ASAV, LDA, X, LDA,
450 $ WORK, LDA, RWORK( 2*NRHS+1 ),
451 $ RESULT( 2 ) )
452 *
453 * Check solution from generated exact solution.
454 *
455 IF( NOFACT .OR. ( PREFAC .AND. LSAME( EQUED,
456 $ 'N' ) ) ) THEN
457 CALL DGET04( N, NRHS, X, LDA, XACT, LDA,
458 $ RCONDC, RESULT( 3 ) )
459 ELSE
460 CALL DGET04( N, NRHS, X, LDA, XACT, LDA,
461 $ ROLDC, RESULT( 3 ) )
462 END IF
463 *
464 * Check the error bounds from iterative
465 * refinement.
466 *
467 CALL DPOT05( UPLO, N, NRHS, ASAV, LDA, B, LDA,
468 $ X, LDA, XACT, LDA, RWORK,
469 $ RWORK( NRHS+1 ), RESULT( 4 ) )
470 ELSE
471 K1 = 6
472 END IF
473 *
474 * Compare RCOND from DPOSVX with the computed value
475 * in RCONDC.
476 *
477 RESULT( 6 ) = DGET06( RCOND, RCONDC )
478 *
479 * Print information about the tests that did not pass
480 * the threshold.
481 *
482 DO 80 K = K1, 6
483 IF( RESULT( K ).GE.THRESH ) THEN
484 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
485 $ CALL ALADHD( NOUT, PATH )
486 IF( PREFAC ) THEN
487 WRITE( NOUT, FMT = 9997 )'DPOSVX', FACT,
488 $ UPLO, N, EQUED, IMAT, K, RESULT( K )
489 ELSE
490 WRITE( NOUT, FMT = 9998 )'DPOSVX', FACT,
491 $ UPLO, N, IMAT, K, RESULT( K )
492 END IF
493 NFAIL = NFAIL + 1
494 END IF
495 80 CONTINUE
496 NRUN = NRUN + 7 - K1
497 90 CONTINUE
498 100 CONTINUE
499 110 CONTINUE
500 120 CONTINUE
501 130 CONTINUE
502 *
503 * Print a summary of the results.
504 *
505 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
506 *
507 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I1,
508 $ ', test(', I1, ')=', G12.5 )
509 9998 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5,
510 $ ', type ', I1, ', test(', I1, ')=', G12.5 )
511 9997 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5,
512 $ ', EQUED=''', A1, ''', type ', I1, ', test(', I1, ') =',
513 $ G12.5 )
514 RETURN
515 *
516 * End of DDRVPO
517 *
518 END