1 SUBROUTINE ZDRVPB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
2 $ A, AFAC, ASAV, B, BSAV, X, XACT, S, WORK,
3 $ RWORK, 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 NVAL( * )
17 DOUBLE PRECISION RWORK( * ), S( * )
18 COMPLEX*16 A( * ), AFAC( * ), ASAV( * ), B( * ),
19 $ BSAV( * ), WORK( * ), X( * ), XACT( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * ZDRVPB tests the driver routines ZPBSV 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) COMPLEX*16 array, dimension (NMAX*NMAX)
58 *
59 * AFAC (workspace) COMPLEX*16 array, dimension (NMAX*NMAX)
60 *
61 * ASAV (workspace) COMPLEX*16 array, dimension (NMAX*NMAX)
62 *
63 * B (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
64 *
65 * BSAV (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
66 *
67 * X (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
68 *
69 * XACT (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
70 *
71 * S (workspace) DOUBLE PRECISION array, dimension (NMAX)
72 *
73 * WORK (workspace) COMPLEX*16 array, dimension
74 * (NMAX*max(3,NRHS))
75 *
76 * RWORK (workspace) DOUBLE PRECISION array, dimension (NMAX+2*NRHS)
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 = 8, NTESTS = 6 )
88 INTEGER NBW
89 PARAMETER ( NBW = 4 )
90 * ..
91 * .. Local Scalars ..
92 LOGICAL EQUIL, NOFACT, PREFAC, ZEROT
93 CHARACTER DIST, EQUED, FACT, PACKIT, TYPE, UPLO, XTYPE
94 CHARACTER*3 PATH
95 INTEGER I, I1, I2, IEQUED, IFACT, IKD, IMAT, IN, INFO,
96 $ IOFF, IUPLO, IW, IZERO, K, K1, KD, KL, KOFF,
97 $ KU, LDA, LDAB, MODE, N, NB, NBMIN, NERRS,
98 $ NFACT, NFAIL, NIMAT, NKD, NRUN, NT
99 DOUBLE PRECISION AINVNM, AMAX, ANORM, CNDNUM, RCOND, RCONDC,
100 $ ROLDC, SCOND
101 * ..
102 * .. Local Arrays ..
103 CHARACTER EQUEDS( 2 ), FACTS( 3 )
104 INTEGER ISEED( 4 ), ISEEDY( 4 ), KDVAL( NBW )
105 DOUBLE PRECISION RESULT( NTESTS )
106 * ..
107 * .. External Functions ..
108 LOGICAL LSAME
109 DOUBLE PRECISION DGET06, ZLANGE, ZLANHB
110 EXTERNAL LSAME, DGET06, ZLANGE, ZLANHB
111 * ..
112 * .. External Subroutines ..
113 EXTERNAL ALADHD, ALAERH, ALASVM, XLAENV, ZCOPY, ZERRVX,
114 $ ZGET04, ZLACPY, ZLAIPD, ZLAQHB, ZLARHS, ZLASET,
115 $ ZLATB4, ZLATMS, ZPBEQU, ZPBSV, ZPBSVX, ZPBT01,
116 $ ZPBT02, ZPBT05, ZPBTRF, ZPBTRS, ZSWAP
117 * ..
118 * .. Intrinsic Functions ..
119 INTRINSIC DCMPLX, MAX, MIN
120 * ..
121 * .. Scalars in Common ..
122 LOGICAL LERR, OK
123 CHARACTER*32 SRNAMT
124 INTEGER INFOT, NUNIT
125 * ..
126 * .. Common blocks ..
127 COMMON / INFOC / INFOT, NUNIT, OK, LERR
128 COMMON / SRNAMC / SRNAMT
129 * ..
130 * .. Data statements ..
131 DATA ISEEDY / 1988, 1989, 1990, 1991 /
132 DATA FACTS / 'F', 'N', 'E' / , EQUEDS / 'N', 'Y' /
133 * ..
134 * .. Executable Statements ..
135 *
136 * Initialize constants and the random number seed.
137 *
138 PATH( 1: 1 ) = 'Zomplex precision'
139 PATH( 2: 3 ) = 'PB'
140 NRUN = 0
141 NFAIL = 0
142 NERRS = 0
143 DO 10 I = 1, 4
144 ISEED( I ) = ISEEDY( I )
145 10 CONTINUE
146 *
147 * Test the error exits
148 *
149 IF( TSTERR )
150 $ CALL ZERRVX( PATH, NOUT )
151 INFOT = 0
152 KDVAL( 1 ) = 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 110 IN = 1, NN
164 N = NVAL( IN )
165 LDA = MAX( N, 1 )
166 XTYPE = 'N'
167 *
168 * Set limits on the number of loop iterations.
169 *
170 NKD = MAX( 1, MIN( N, 4 ) )
171 NIMAT = NTYPES
172 IF( N.EQ.0 )
173 $ NIMAT = 1
174 *
175 KDVAL( 2 ) = N + ( N+1 ) / 4
176 KDVAL( 3 ) = ( 3*N-1 ) / 4
177 KDVAL( 4 ) = ( N+1 ) / 4
178 *
179 DO 100 IKD = 1, NKD
180 *
181 * Do for KD = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This order
182 * makes it easier to skip redundant values for small values
183 * of N.
184 *
185 KD = KDVAL( IKD )
186 LDAB = KD + 1
187 *
188 * Do first for UPLO = 'U', then for UPLO = 'L'
189 *
190 DO 90 IUPLO = 1, 2
191 KOFF = 1
192 IF( IUPLO.EQ.1 ) THEN
193 UPLO = 'U'
194 PACKIT = 'Q'
195 KOFF = MAX( 1, KD+2-N )
196 ELSE
197 UPLO = 'L'
198 PACKIT = 'B'
199 END IF
200 *
201 DO 80 IMAT = 1, NIMAT
202 *
203 * Do the tests only if DOTYPE( IMAT ) is true.
204 *
205 IF( .NOT.DOTYPE( IMAT ) )
206 $ GO TO 80
207 *
208 * Skip types 2, 3, or 4 if the matrix size is too small.
209 *
210 ZEROT = IMAT.GE.2 .AND. IMAT.LE.4
211 IF( ZEROT .AND. N.LT.IMAT-1 )
212 $ GO TO 80
213 *
214 IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 1 ) ) THEN
215 *
216 * Set up parameters with ZLATB4 and generate a test
217 * matrix with ZLATMS.
218 *
219 CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM,
220 $ MODE, CNDNUM, DIST )
221 *
222 SRNAMT = 'ZLATMS'
223 CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
224 $ CNDNUM, ANORM, KD, KD, PACKIT,
225 $ A( KOFF ), LDAB, WORK, INFO )
226 *
227 * Check error code from ZLATMS.
228 *
229 IF( INFO.NE.0 ) THEN
230 CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N,
231 $ N, -1, -1, -1, IMAT, NFAIL, NERRS,
232 $ NOUT )
233 GO TO 80
234 END IF
235 ELSE IF( IZERO.GT.0 ) THEN
236 *
237 * Use the same matrix for types 3 and 4 as for type
238 * 2 by copying back the zeroed out column,
239 *
240 IW = 2*LDA + 1
241 IF( IUPLO.EQ.1 ) THEN
242 IOFF = ( IZERO-1 )*LDAB + KD + 1
243 CALL ZCOPY( IZERO-I1, WORK( IW ), 1,
244 $ A( IOFF-IZERO+I1 ), 1 )
245 IW = IW + IZERO - I1
246 CALL ZCOPY( I2-IZERO+1, WORK( IW ), 1,
247 $ A( IOFF ), MAX( LDAB-1, 1 ) )
248 ELSE
249 IOFF = ( I1-1 )*LDAB + 1
250 CALL ZCOPY( IZERO-I1, WORK( IW ), 1,
251 $ A( IOFF+IZERO-I1 ),
252 $ MAX( LDAB-1, 1 ) )
253 IOFF = ( IZERO-1 )*LDAB + 1
254 IW = IW + IZERO - I1
255 CALL ZCOPY( I2-IZERO+1, WORK( IW ), 1,
256 $ A( IOFF ), 1 )
257 END IF
258 END IF
259 *
260 * For types 2-4, zero one row and column of the matrix
261 * to test that INFO is returned correctly.
262 *
263 IZERO = 0
264 IF( ZEROT ) THEN
265 IF( IMAT.EQ.2 ) THEN
266 IZERO = 1
267 ELSE IF( IMAT.EQ.3 ) THEN
268 IZERO = N
269 ELSE
270 IZERO = N / 2 + 1
271 END IF
272 *
273 * Save the zeroed out row and column in WORK(*,3)
274 *
275 IW = 2*LDA
276 DO 20 I = 1, MIN( 2*KD+1, N )
277 WORK( IW+I ) = ZERO
278 20 CONTINUE
279 IW = IW + 1
280 I1 = MAX( IZERO-KD, 1 )
281 I2 = MIN( IZERO+KD, N )
282 *
283 IF( IUPLO.EQ.1 ) THEN
284 IOFF = ( IZERO-1 )*LDAB + KD + 1
285 CALL ZSWAP( IZERO-I1, A( IOFF-IZERO+I1 ), 1,
286 $ WORK( IW ), 1 )
287 IW = IW + IZERO - I1
288 CALL ZSWAP( I2-IZERO+1, A( IOFF ),
289 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 )
290 ELSE
291 IOFF = ( I1-1 )*LDAB + 1
292 CALL ZSWAP( IZERO-I1, A( IOFF+IZERO-I1 ),
293 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 )
294 IOFF = ( IZERO-1 )*LDAB + 1
295 IW = IW + IZERO - I1
296 CALL ZSWAP( I2-IZERO+1, A( IOFF ), 1,
297 $ WORK( IW ), 1 )
298 END IF
299 END IF
300 *
301 * Set the imaginary part of the diagonals.
302 *
303 IF( IUPLO.EQ.1 ) THEN
304 CALL ZLAIPD( N, A( KD+1 ), LDAB, 0 )
305 ELSE
306 CALL ZLAIPD( N, A( 1 ), LDAB, 0 )
307 END IF
308 *
309 * Save a copy of the matrix A in ASAV.
310 *
311 CALL ZLACPY( 'Full', KD+1, N, A, LDAB, ASAV, LDAB )
312 *
313 DO 70 IEQUED = 1, 2
314 EQUED = EQUEDS( IEQUED )
315 IF( IEQUED.EQ.1 ) THEN
316 NFACT = 3
317 ELSE
318 NFACT = 1
319 END IF
320 *
321 DO 60 IFACT = 1, NFACT
322 FACT = FACTS( IFACT )
323 PREFAC = LSAME( FACT, 'F' )
324 NOFACT = LSAME( FACT, 'N' )
325 EQUIL = LSAME( FACT, 'E' )
326 *
327 IF( ZEROT ) THEN
328 IF( PREFAC )
329 $ GO TO 60
330 RCONDC = ZERO
331 *
332 ELSE IF( .NOT.LSAME( FACT, 'N' ) ) THEN
333 *
334 * Compute the condition number for comparison
335 * with the value returned by ZPBSVX (FACT =
336 * 'N' reuses the condition number from the
337 * previous iteration with FACT = 'F').
338 *
339 CALL ZLACPY( 'Full', KD+1, N, ASAV, LDAB,
340 $ AFAC, LDAB )
341 IF( EQUIL .OR. IEQUED.GT.1 ) THEN
342 *
343 * Compute row and column scale factors to
344 * equilibrate the matrix A.
345 *
346 CALL ZPBEQU( UPLO, N, KD, AFAC, LDAB, S,
347 $ SCOND, AMAX, INFO )
348 IF( INFO.EQ.0 .AND. N.GT.0 ) THEN
349 IF( IEQUED.GT.1 )
350 $ SCOND = ZERO
351 *
352 * Equilibrate the matrix.
353 *
354 CALL ZLAQHB( UPLO, N, KD, AFAC, LDAB,
355 $ S, SCOND, AMAX, EQUED )
356 END IF
357 END IF
358 *
359 * Save the condition number of the
360 * non-equilibrated system for use in ZGET04.
361 *
362 IF( EQUIL )
363 $ ROLDC = RCONDC
364 *
365 * Compute the 1-norm of A.
366 *
367 ANORM = ZLANHB( '1', UPLO, N, KD, AFAC, LDAB,
368 $ RWORK )
369 *
370 * Factor the matrix A.
371 *
372 CALL ZPBTRF( UPLO, N, KD, AFAC, LDAB, INFO )
373 *
374 * Form the inverse of A.
375 *
376 CALL ZLASET( 'Full', N, N, DCMPLX( ZERO ),
377 $ DCMPLX( ONE ), A, LDA )
378 SRNAMT = 'ZPBTRS'
379 CALL ZPBTRS( UPLO, N, KD, N, AFAC, LDAB, A,
380 $ LDA, INFO )
381 *
382 * Compute the 1-norm condition number of A.
383 *
384 AINVNM = ZLANGE( '1', N, N, A, LDA, RWORK )
385 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
386 RCONDC = ONE
387 ELSE
388 RCONDC = ( ONE / ANORM ) / AINVNM
389 END IF
390 END IF
391 *
392 * Restore the matrix A.
393 *
394 CALL ZLACPY( 'Full', KD+1, N, ASAV, LDAB, A,
395 $ LDAB )
396 *
397 * Form an exact solution and set the right hand
398 * side.
399 *
400 SRNAMT = 'ZLARHS'
401 CALL ZLARHS( PATH, XTYPE, UPLO, ' ', N, N, KD,
402 $ KD, NRHS, A, LDAB, XACT, LDA, B,
403 $ LDA, ISEED, INFO )
404 XTYPE = 'C'
405 CALL ZLACPY( 'Full', N, NRHS, B, LDA, BSAV,
406 $ LDA )
407 *
408 IF( NOFACT ) THEN
409 *
410 * --- Test ZPBSV ---
411 *
412 * Compute the L*L' or U'*U factorization of the
413 * matrix and solve the system.
414 *
415 CALL ZLACPY( 'Full', KD+1, N, A, LDAB, AFAC,
416 $ LDAB )
417 CALL ZLACPY( 'Full', N, NRHS, B, LDA, X,
418 $ LDA )
419 *
420 SRNAMT = 'ZPBSV '
421 CALL ZPBSV( UPLO, N, KD, NRHS, AFAC, LDAB, X,
422 $ LDA, INFO )
423 *
424 * Check error code from ZPBSV .
425 *
426 IF( INFO.NE.IZERO ) THEN
427 CALL ALAERH( PATH, 'ZPBSV ', INFO, IZERO,
428 $ UPLO, N, N, KD, KD, NRHS,
429 $ IMAT, NFAIL, NERRS, NOUT )
430 GO TO 40
431 ELSE IF( INFO.NE.0 ) THEN
432 GO TO 40
433 END IF
434 *
435 * Reconstruct matrix from factors and compute
436 * residual.
437 *
438 CALL ZPBT01( UPLO, N, KD, A, LDAB, AFAC,
439 $ LDAB, RWORK, RESULT( 1 ) )
440 *
441 * Compute residual of the computed solution.
442 *
443 CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK,
444 $ LDA )
445 CALL ZPBT02( UPLO, N, KD, NRHS, A, LDAB, X,
446 $ LDA, WORK, LDA, RWORK,
447 $ RESULT( 2 ) )
448 *
449 * Check solution from generated exact solution.
450 *
451 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA,
452 $ RCONDC, RESULT( 3 ) )
453 NT = 3
454 *
455 * Print information about the tests that did
456 * not pass the threshold.
457 *
458 DO 30 K = 1, NT
459 IF( RESULT( K ).GE.THRESH ) THEN
460 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
461 $ CALL ALADHD( NOUT, PATH )
462 WRITE( NOUT, FMT = 9999 )'ZPBSV ',
463 $ UPLO, N, KD, IMAT, K, RESULT( K )
464 NFAIL = NFAIL + 1
465 END IF
466 30 CONTINUE
467 NRUN = NRUN + NT
468 40 CONTINUE
469 END IF
470 *
471 * --- Test ZPBSVX ---
472 *
473 IF( .NOT.PREFAC )
474 $ CALL ZLASET( 'Full', KD+1, N, DCMPLX( ZERO ),
475 $ DCMPLX( ZERO ), AFAC, LDAB )
476 CALL ZLASET( 'Full', N, NRHS, DCMPLX( ZERO ),
477 $ DCMPLX( ZERO ), X, LDA )
478 IF( IEQUED.GT.1 .AND. N.GT.0 ) THEN
479 *
480 * Equilibrate the matrix if FACT='F' and
481 * EQUED='Y'
482 *
483 CALL ZLAQHB( UPLO, N, KD, A, LDAB, S, SCOND,
484 $ AMAX, EQUED )
485 END IF
486 *
487 * Solve the system and compute the condition
488 * number and error bounds using ZPBSVX.
489 *
490 SRNAMT = 'ZPBSVX'
491 CALL ZPBSVX( FACT, UPLO, N, KD, NRHS, A, LDAB,
492 $ AFAC, LDAB, EQUED, S, B, LDA, X,
493 $ LDA, RCOND, RWORK, RWORK( NRHS+1 ),
494 $ WORK, RWORK( 2*NRHS+1 ), INFO )
495 *
496 * Check the error code from ZPBSVX.
497 *
498 IF( INFO.NE.IZERO ) THEN
499 CALL ALAERH( PATH, 'ZPBSVX', INFO, IZERO,
500 $ FACT // UPLO, N, N, KD, KD,
501 $ NRHS, IMAT, NFAIL, NERRS, NOUT )
502 GO TO 60
503 END IF
504 *
505 IF( INFO.EQ.0 ) THEN
506 IF( .NOT.PREFAC ) THEN
507 *
508 * Reconstruct matrix from factors and
509 * compute residual.
510 *
511 CALL ZPBT01( UPLO, N, KD, A, LDAB, AFAC,
512 $ LDAB, RWORK( 2*NRHS+1 ),
513 $ RESULT( 1 ) )
514 K1 = 1
515 ELSE
516 K1 = 2
517 END IF
518 *
519 * Compute residual of the computed solution.
520 *
521 CALL ZLACPY( 'Full', N, NRHS, BSAV, LDA,
522 $ WORK, LDA )
523 CALL ZPBT02( UPLO, N, KD, NRHS, ASAV, LDAB,
524 $ X, LDA, WORK, LDA,
525 $ RWORK( 2*NRHS+1 ), RESULT( 2 ) )
526 *
527 * Check solution from generated exact solution.
528 *
529 IF( NOFACT .OR. ( PREFAC .AND. LSAME( EQUED,
530 $ 'N' ) ) ) THEN
531 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA,
532 $ RCONDC, RESULT( 3 ) )
533 ELSE
534 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA,
535 $ ROLDC, RESULT( 3 ) )
536 END IF
537 *
538 * Check the error bounds from iterative
539 * refinement.
540 *
541 CALL ZPBT05( UPLO, N, KD, NRHS, ASAV, LDAB,
542 $ B, LDA, X, LDA, XACT, LDA,
543 $ RWORK, RWORK( NRHS+1 ),
544 $ RESULT( 4 ) )
545 ELSE
546 K1 = 6
547 END IF
548 *
549 * Compare RCOND from ZPBSVX with the computed
550 * value in RCONDC.
551 *
552 RESULT( 6 ) = DGET06( RCOND, RCONDC )
553 *
554 * Print information about the tests that did not
555 * pass the threshold.
556 *
557 DO 50 K = K1, 6
558 IF( RESULT( K ).GE.THRESH ) THEN
559 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
560 $ CALL ALADHD( NOUT, PATH )
561 IF( PREFAC ) THEN
562 WRITE( NOUT, FMT = 9997 )'ZPBSVX',
563 $ FACT, UPLO, N, KD, EQUED, IMAT, K,
564 $ RESULT( K )
565 ELSE
566 WRITE( NOUT, FMT = 9998 )'ZPBSVX',
567 $ FACT, UPLO, N, KD, IMAT, K,
568 $ RESULT( K )
569 END IF
570 NFAIL = NFAIL + 1
571 END IF
572 50 CONTINUE
573 NRUN = NRUN + 7 - K1
574 60 CONTINUE
575 70 CONTINUE
576 80 CONTINUE
577 90 CONTINUE
578 100 CONTINUE
579 110 CONTINUE
580 *
581 * Print a summary of the results.
582 *
583 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
584 *
585 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', KD =', I5,
586 $ ', type ', I1, ', test(', I1, ')=', G12.5 )
587 9998 FORMAT( 1X, A, '( ''', A1, ''', ''', A1, ''', ', I5, ', ', I5,
588 $ ', ... ), type ', I1, ', test(', I1, ')=', G12.5 )
589 9997 FORMAT( 1X, A, '( ''', A1, ''', ''', A1, ''', ', I5, ', ', I5,
590 $ ', ... ), EQUED=''', A1, ''', type ', I1, ', test(', I1,
591 $ ')=', G12.5 )
592 RETURN
593 *
594 * End of ZDRVPB
595 *
596 END
2 $ A, AFAC, ASAV, B, BSAV, X, XACT, S, WORK,
3 $ RWORK, 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 NVAL( * )
17 DOUBLE PRECISION RWORK( * ), S( * )
18 COMPLEX*16 A( * ), AFAC( * ), ASAV( * ), B( * ),
19 $ BSAV( * ), WORK( * ), X( * ), XACT( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * ZDRVPB tests the driver routines ZPBSV 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) COMPLEX*16 array, dimension (NMAX*NMAX)
58 *
59 * AFAC (workspace) COMPLEX*16 array, dimension (NMAX*NMAX)
60 *
61 * ASAV (workspace) COMPLEX*16 array, dimension (NMAX*NMAX)
62 *
63 * B (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
64 *
65 * BSAV (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
66 *
67 * X (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
68 *
69 * XACT (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
70 *
71 * S (workspace) DOUBLE PRECISION array, dimension (NMAX)
72 *
73 * WORK (workspace) COMPLEX*16 array, dimension
74 * (NMAX*max(3,NRHS))
75 *
76 * RWORK (workspace) DOUBLE PRECISION array, dimension (NMAX+2*NRHS)
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 = 8, NTESTS = 6 )
88 INTEGER NBW
89 PARAMETER ( NBW = 4 )
90 * ..
91 * .. Local Scalars ..
92 LOGICAL EQUIL, NOFACT, PREFAC, ZEROT
93 CHARACTER DIST, EQUED, FACT, PACKIT, TYPE, UPLO, XTYPE
94 CHARACTER*3 PATH
95 INTEGER I, I1, I2, IEQUED, IFACT, IKD, IMAT, IN, INFO,
96 $ IOFF, IUPLO, IW, IZERO, K, K1, KD, KL, KOFF,
97 $ KU, LDA, LDAB, MODE, N, NB, NBMIN, NERRS,
98 $ NFACT, NFAIL, NIMAT, NKD, NRUN, NT
99 DOUBLE PRECISION AINVNM, AMAX, ANORM, CNDNUM, RCOND, RCONDC,
100 $ ROLDC, SCOND
101 * ..
102 * .. Local Arrays ..
103 CHARACTER EQUEDS( 2 ), FACTS( 3 )
104 INTEGER ISEED( 4 ), ISEEDY( 4 ), KDVAL( NBW )
105 DOUBLE PRECISION RESULT( NTESTS )
106 * ..
107 * .. External Functions ..
108 LOGICAL LSAME
109 DOUBLE PRECISION DGET06, ZLANGE, ZLANHB
110 EXTERNAL LSAME, DGET06, ZLANGE, ZLANHB
111 * ..
112 * .. External Subroutines ..
113 EXTERNAL ALADHD, ALAERH, ALASVM, XLAENV, ZCOPY, ZERRVX,
114 $ ZGET04, ZLACPY, ZLAIPD, ZLAQHB, ZLARHS, ZLASET,
115 $ ZLATB4, ZLATMS, ZPBEQU, ZPBSV, ZPBSVX, ZPBT01,
116 $ ZPBT02, ZPBT05, ZPBTRF, ZPBTRS, ZSWAP
117 * ..
118 * .. Intrinsic Functions ..
119 INTRINSIC DCMPLX, MAX, MIN
120 * ..
121 * .. Scalars in Common ..
122 LOGICAL LERR, OK
123 CHARACTER*32 SRNAMT
124 INTEGER INFOT, NUNIT
125 * ..
126 * .. Common blocks ..
127 COMMON / INFOC / INFOT, NUNIT, OK, LERR
128 COMMON / SRNAMC / SRNAMT
129 * ..
130 * .. Data statements ..
131 DATA ISEEDY / 1988, 1989, 1990, 1991 /
132 DATA FACTS / 'F', 'N', 'E' / , EQUEDS / 'N', 'Y' /
133 * ..
134 * .. Executable Statements ..
135 *
136 * Initialize constants and the random number seed.
137 *
138 PATH( 1: 1 ) = 'Zomplex precision'
139 PATH( 2: 3 ) = 'PB'
140 NRUN = 0
141 NFAIL = 0
142 NERRS = 0
143 DO 10 I = 1, 4
144 ISEED( I ) = ISEEDY( I )
145 10 CONTINUE
146 *
147 * Test the error exits
148 *
149 IF( TSTERR )
150 $ CALL ZERRVX( PATH, NOUT )
151 INFOT = 0
152 KDVAL( 1 ) = 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 110 IN = 1, NN
164 N = NVAL( IN )
165 LDA = MAX( N, 1 )
166 XTYPE = 'N'
167 *
168 * Set limits on the number of loop iterations.
169 *
170 NKD = MAX( 1, MIN( N, 4 ) )
171 NIMAT = NTYPES
172 IF( N.EQ.0 )
173 $ NIMAT = 1
174 *
175 KDVAL( 2 ) = N + ( N+1 ) / 4
176 KDVAL( 3 ) = ( 3*N-1 ) / 4
177 KDVAL( 4 ) = ( N+1 ) / 4
178 *
179 DO 100 IKD = 1, NKD
180 *
181 * Do for KD = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This order
182 * makes it easier to skip redundant values for small values
183 * of N.
184 *
185 KD = KDVAL( IKD )
186 LDAB = KD + 1
187 *
188 * Do first for UPLO = 'U', then for UPLO = 'L'
189 *
190 DO 90 IUPLO = 1, 2
191 KOFF = 1
192 IF( IUPLO.EQ.1 ) THEN
193 UPLO = 'U'
194 PACKIT = 'Q'
195 KOFF = MAX( 1, KD+2-N )
196 ELSE
197 UPLO = 'L'
198 PACKIT = 'B'
199 END IF
200 *
201 DO 80 IMAT = 1, NIMAT
202 *
203 * Do the tests only if DOTYPE( IMAT ) is true.
204 *
205 IF( .NOT.DOTYPE( IMAT ) )
206 $ GO TO 80
207 *
208 * Skip types 2, 3, or 4 if the matrix size is too small.
209 *
210 ZEROT = IMAT.GE.2 .AND. IMAT.LE.4
211 IF( ZEROT .AND. N.LT.IMAT-1 )
212 $ GO TO 80
213 *
214 IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 1 ) ) THEN
215 *
216 * Set up parameters with ZLATB4 and generate a test
217 * matrix with ZLATMS.
218 *
219 CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM,
220 $ MODE, CNDNUM, DIST )
221 *
222 SRNAMT = 'ZLATMS'
223 CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
224 $ CNDNUM, ANORM, KD, KD, PACKIT,
225 $ A( KOFF ), LDAB, WORK, INFO )
226 *
227 * Check error code from ZLATMS.
228 *
229 IF( INFO.NE.0 ) THEN
230 CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N,
231 $ N, -1, -1, -1, IMAT, NFAIL, NERRS,
232 $ NOUT )
233 GO TO 80
234 END IF
235 ELSE IF( IZERO.GT.0 ) THEN
236 *
237 * Use the same matrix for types 3 and 4 as for type
238 * 2 by copying back the zeroed out column,
239 *
240 IW = 2*LDA + 1
241 IF( IUPLO.EQ.1 ) THEN
242 IOFF = ( IZERO-1 )*LDAB + KD + 1
243 CALL ZCOPY( IZERO-I1, WORK( IW ), 1,
244 $ A( IOFF-IZERO+I1 ), 1 )
245 IW = IW + IZERO - I1
246 CALL ZCOPY( I2-IZERO+1, WORK( IW ), 1,
247 $ A( IOFF ), MAX( LDAB-1, 1 ) )
248 ELSE
249 IOFF = ( I1-1 )*LDAB + 1
250 CALL ZCOPY( IZERO-I1, WORK( IW ), 1,
251 $ A( IOFF+IZERO-I1 ),
252 $ MAX( LDAB-1, 1 ) )
253 IOFF = ( IZERO-1 )*LDAB + 1
254 IW = IW + IZERO - I1
255 CALL ZCOPY( I2-IZERO+1, WORK( IW ), 1,
256 $ A( IOFF ), 1 )
257 END IF
258 END IF
259 *
260 * For types 2-4, zero one row and column of the matrix
261 * to test that INFO is returned correctly.
262 *
263 IZERO = 0
264 IF( ZEROT ) THEN
265 IF( IMAT.EQ.2 ) THEN
266 IZERO = 1
267 ELSE IF( IMAT.EQ.3 ) THEN
268 IZERO = N
269 ELSE
270 IZERO = N / 2 + 1
271 END IF
272 *
273 * Save the zeroed out row and column in WORK(*,3)
274 *
275 IW = 2*LDA
276 DO 20 I = 1, MIN( 2*KD+1, N )
277 WORK( IW+I ) = ZERO
278 20 CONTINUE
279 IW = IW + 1
280 I1 = MAX( IZERO-KD, 1 )
281 I2 = MIN( IZERO+KD, N )
282 *
283 IF( IUPLO.EQ.1 ) THEN
284 IOFF = ( IZERO-1 )*LDAB + KD + 1
285 CALL ZSWAP( IZERO-I1, A( IOFF-IZERO+I1 ), 1,
286 $ WORK( IW ), 1 )
287 IW = IW + IZERO - I1
288 CALL ZSWAP( I2-IZERO+1, A( IOFF ),
289 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 )
290 ELSE
291 IOFF = ( I1-1 )*LDAB + 1
292 CALL ZSWAP( IZERO-I1, A( IOFF+IZERO-I1 ),
293 $ MAX( LDAB-1, 1 ), WORK( IW ), 1 )
294 IOFF = ( IZERO-1 )*LDAB + 1
295 IW = IW + IZERO - I1
296 CALL ZSWAP( I2-IZERO+1, A( IOFF ), 1,
297 $ WORK( IW ), 1 )
298 END IF
299 END IF
300 *
301 * Set the imaginary part of the diagonals.
302 *
303 IF( IUPLO.EQ.1 ) THEN
304 CALL ZLAIPD( N, A( KD+1 ), LDAB, 0 )
305 ELSE
306 CALL ZLAIPD( N, A( 1 ), LDAB, 0 )
307 END IF
308 *
309 * Save a copy of the matrix A in ASAV.
310 *
311 CALL ZLACPY( 'Full', KD+1, N, A, LDAB, ASAV, LDAB )
312 *
313 DO 70 IEQUED = 1, 2
314 EQUED = EQUEDS( IEQUED )
315 IF( IEQUED.EQ.1 ) THEN
316 NFACT = 3
317 ELSE
318 NFACT = 1
319 END IF
320 *
321 DO 60 IFACT = 1, NFACT
322 FACT = FACTS( IFACT )
323 PREFAC = LSAME( FACT, 'F' )
324 NOFACT = LSAME( FACT, 'N' )
325 EQUIL = LSAME( FACT, 'E' )
326 *
327 IF( ZEROT ) THEN
328 IF( PREFAC )
329 $ GO TO 60
330 RCONDC = ZERO
331 *
332 ELSE IF( .NOT.LSAME( FACT, 'N' ) ) THEN
333 *
334 * Compute the condition number for comparison
335 * with the value returned by ZPBSVX (FACT =
336 * 'N' reuses the condition number from the
337 * previous iteration with FACT = 'F').
338 *
339 CALL ZLACPY( 'Full', KD+1, N, ASAV, LDAB,
340 $ AFAC, LDAB )
341 IF( EQUIL .OR. IEQUED.GT.1 ) THEN
342 *
343 * Compute row and column scale factors to
344 * equilibrate the matrix A.
345 *
346 CALL ZPBEQU( UPLO, N, KD, AFAC, LDAB, S,
347 $ SCOND, AMAX, INFO )
348 IF( INFO.EQ.0 .AND. N.GT.0 ) THEN
349 IF( IEQUED.GT.1 )
350 $ SCOND = ZERO
351 *
352 * Equilibrate the matrix.
353 *
354 CALL ZLAQHB( UPLO, N, KD, AFAC, LDAB,
355 $ S, SCOND, AMAX, EQUED )
356 END IF
357 END IF
358 *
359 * Save the condition number of the
360 * non-equilibrated system for use in ZGET04.
361 *
362 IF( EQUIL )
363 $ ROLDC = RCONDC
364 *
365 * Compute the 1-norm of A.
366 *
367 ANORM = ZLANHB( '1', UPLO, N, KD, AFAC, LDAB,
368 $ RWORK )
369 *
370 * Factor the matrix A.
371 *
372 CALL ZPBTRF( UPLO, N, KD, AFAC, LDAB, INFO )
373 *
374 * Form the inverse of A.
375 *
376 CALL ZLASET( 'Full', N, N, DCMPLX( ZERO ),
377 $ DCMPLX( ONE ), A, LDA )
378 SRNAMT = 'ZPBTRS'
379 CALL ZPBTRS( UPLO, N, KD, N, AFAC, LDAB, A,
380 $ LDA, INFO )
381 *
382 * Compute the 1-norm condition number of A.
383 *
384 AINVNM = ZLANGE( '1', N, N, A, LDA, RWORK )
385 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
386 RCONDC = ONE
387 ELSE
388 RCONDC = ( ONE / ANORM ) / AINVNM
389 END IF
390 END IF
391 *
392 * Restore the matrix A.
393 *
394 CALL ZLACPY( 'Full', KD+1, N, ASAV, LDAB, A,
395 $ LDAB )
396 *
397 * Form an exact solution and set the right hand
398 * side.
399 *
400 SRNAMT = 'ZLARHS'
401 CALL ZLARHS( PATH, XTYPE, UPLO, ' ', N, N, KD,
402 $ KD, NRHS, A, LDAB, XACT, LDA, B,
403 $ LDA, ISEED, INFO )
404 XTYPE = 'C'
405 CALL ZLACPY( 'Full', N, NRHS, B, LDA, BSAV,
406 $ LDA )
407 *
408 IF( NOFACT ) THEN
409 *
410 * --- Test ZPBSV ---
411 *
412 * Compute the L*L' or U'*U factorization of the
413 * matrix and solve the system.
414 *
415 CALL ZLACPY( 'Full', KD+1, N, A, LDAB, AFAC,
416 $ LDAB )
417 CALL ZLACPY( 'Full', N, NRHS, B, LDA, X,
418 $ LDA )
419 *
420 SRNAMT = 'ZPBSV '
421 CALL ZPBSV( UPLO, N, KD, NRHS, AFAC, LDAB, X,
422 $ LDA, INFO )
423 *
424 * Check error code from ZPBSV .
425 *
426 IF( INFO.NE.IZERO ) THEN
427 CALL ALAERH( PATH, 'ZPBSV ', INFO, IZERO,
428 $ UPLO, N, N, KD, KD, NRHS,
429 $ IMAT, NFAIL, NERRS, NOUT )
430 GO TO 40
431 ELSE IF( INFO.NE.0 ) THEN
432 GO TO 40
433 END IF
434 *
435 * Reconstruct matrix from factors and compute
436 * residual.
437 *
438 CALL ZPBT01( UPLO, N, KD, A, LDAB, AFAC,
439 $ LDAB, RWORK, RESULT( 1 ) )
440 *
441 * Compute residual of the computed solution.
442 *
443 CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK,
444 $ LDA )
445 CALL ZPBT02( UPLO, N, KD, NRHS, A, LDAB, X,
446 $ LDA, WORK, LDA, RWORK,
447 $ RESULT( 2 ) )
448 *
449 * Check solution from generated exact solution.
450 *
451 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA,
452 $ RCONDC, RESULT( 3 ) )
453 NT = 3
454 *
455 * Print information about the tests that did
456 * not pass the threshold.
457 *
458 DO 30 K = 1, NT
459 IF( RESULT( K ).GE.THRESH ) THEN
460 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
461 $ CALL ALADHD( NOUT, PATH )
462 WRITE( NOUT, FMT = 9999 )'ZPBSV ',
463 $ UPLO, N, KD, IMAT, K, RESULT( K )
464 NFAIL = NFAIL + 1
465 END IF
466 30 CONTINUE
467 NRUN = NRUN + NT
468 40 CONTINUE
469 END IF
470 *
471 * --- Test ZPBSVX ---
472 *
473 IF( .NOT.PREFAC )
474 $ CALL ZLASET( 'Full', KD+1, N, DCMPLX( ZERO ),
475 $ DCMPLX( ZERO ), AFAC, LDAB )
476 CALL ZLASET( 'Full', N, NRHS, DCMPLX( ZERO ),
477 $ DCMPLX( ZERO ), X, LDA )
478 IF( IEQUED.GT.1 .AND. N.GT.0 ) THEN
479 *
480 * Equilibrate the matrix if FACT='F' and
481 * EQUED='Y'
482 *
483 CALL ZLAQHB( UPLO, N, KD, A, LDAB, S, SCOND,
484 $ AMAX, EQUED )
485 END IF
486 *
487 * Solve the system and compute the condition
488 * number and error bounds using ZPBSVX.
489 *
490 SRNAMT = 'ZPBSVX'
491 CALL ZPBSVX( FACT, UPLO, N, KD, NRHS, A, LDAB,
492 $ AFAC, LDAB, EQUED, S, B, LDA, X,
493 $ LDA, RCOND, RWORK, RWORK( NRHS+1 ),
494 $ WORK, RWORK( 2*NRHS+1 ), INFO )
495 *
496 * Check the error code from ZPBSVX.
497 *
498 IF( INFO.NE.IZERO ) THEN
499 CALL ALAERH( PATH, 'ZPBSVX', INFO, IZERO,
500 $ FACT // UPLO, N, N, KD, KD,
501 $ NRHS, IMAT, NFAIL, NERRS, NOUT )
502 GO TO 60
503 END IF
504 *
505 IF( INFO.EQ.0 ) THEN
506 IF( .NOT.PREFAC ) THEN
507 *
508 * Reconstruct matrix from factors and
509 * compute residual.
510 *
511 CALL ZPBT01( UPLO, N, KD, A, LDAB, AFAC,
512 $ LDAB, RWORK( 2*NRHS+1 ),
513 $ RESULT( 1 ) )
514 K1 = 1
515 ELSE
516 K1 = 2
517 END IF
518 *
519 * Compute residual of the computed solution.
520 *
521 CALL ZLACPY( 'Full', N, NRHS, BSAV, LDA,
522 $ WORK, LDA )
523 CALL ZPBT02( UPLO, N, KD, NRHS, ASAV, LDAB,
524 $ X, LDA, WORK, LDA,
525 $ RWORK( 2*NRHS+1 ), RESULT( 2 ) )
526 *
527 * Check solution from generated exact solution.
528 *
529 IF( NOFACT .OR. ( PREFAC .AND. LSAME( EQUED,
530 $ 'N' ) ) ) THEN
531 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA,
532 $ RCONDC, RESULT( 3 ) )
533 ELSE
534 CALL ZGET04( N, NRHS, X, LDA, XACT, LDA,
535 $ ROLDC, RESULT( 3 ) )
536 END IF
537 *
538 * Check the error bounds from iterative
539 * refinement.
540 *
541 CALL ZPBT05( UPLO, N, KD, NRHS, ASAV, LDAB,
542 $ B, LDA, X, LDA, XACT, LDA,
543 $ RWORK, RWORK( NRHS+1 ),
544 $ RESULT( 4 ) )
545 ELSE
546 K1 = 6
547 END IF
548 *
549 * Compare RCOND from ZPBSVX with the computed
550 * value in RCONDC.
551 *
552 RESULT( 6 ) = DGET06( RCOND, RCONDC )
553 *
554 * Print information about the tests that did not
555 * pass the threshold.
556 *
557 DO 50 K = K1, 6
558 IF( RESULT( K ).GE.THRESH ) THEN
559 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
560 $ CALL ALADHD( NOUT, PATH )
561 IF( PREFAC ) THEN
562 WRITE( NOUT, FMT = 9997 )'ZPBSVX',
563 $ FACT, UPLO, N, KD, EQUED, IMAT, K,
564 $ RESULT( K )
565 ELSE
566 WRITE( NOUT, FMT = 9998 )'ZPBSVX',
567 $ FACT, UPLO, N, KD, IMAT, K,
568 $ RESULT( K )
569 END IF
570 NFAIL = NFAIL + 1
571 END IF
572 50 CONTINUE
573 NRUN = NRUN + 7 - K1
574 60 CONTINUE
575 70 CONTINUE
576 80 CONTINUE
577 90 CONTINUE
578 100 CONTINUE
579 110 CONTINUE
580 *
581 * Print a summary of the results.
582 *
583 CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
584 *
585 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', KD =', I5,
586 $ ', type ', I1, ', test(', I1, ')=', G12.5 )
587 9998 FORMAT( 1X, A, '( ''', A1, ''', ''', A1, ''', ', I5, ', ', I5,
588 $ ', ... ), type ', I1, ', test(', I1, ')=', G12.5 )
589 9997 FORMAT( 1X, A, '( ''', A1, ''', ''', A1, ''', ', I5, ', ', I5,
590 $ ', ... ), EQUED=''', A1, ''', type ', I1, ', test(', I1,
591 $ ')=', G12.5 )
592 RETURN
593 *
594 * End of ZDRVPB
595 *
596 END