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