1 SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS,
2 $ NSVAL, THRESH, TSTERR, A, LA, AFAC, LAFAC, B,
3 $ X, XACT, WORK, 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 LA, LAFAC, NM, NN, NNB, NNS, NOUT
12 REAL THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL DOTYPE( * )
16 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
17 $ NVAL( * )
18 REAL A( * ), AFAC( * ), B( * ), RWORK( * ),
19 $ WORK( * ), X( * ), XACT( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * SCHKGB tests SGBTRF, -TRS, -RFS, and -CON
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 * NM (input) INTEGER
36 * The number of values of M contained in the vector MVAL.
37 *
38 * MVAL (input) INTEGER array, dimension (NM)
39 * The values of the matrix row dimension M.
40 *
41 * NN (input) INTEGER
42 * The number of values of N contained in the vector NVAL.
43 *
44 * NVAL (input) INTEGER array, dimension (NN)
45 * The values of the matrix column dimension N.
46 *
47 * NNB (input) INTEGER
48 * The number of values of NB contained in the vector NBVAL.
49 *
50 * NBVAL (input) INTEGER array, dimension (NNB)
51 * The values of the blocksize NB.
52 *
53 * NNS (input) INTEGER
54 * The number of values of NRHS contained in the vector NSVAL.
55 *
56 * NSVAL (input) INTEGER array, dimension (NNS)
57 * The values of the number of right hand sides NRHS.
58 *
59 * THRESH (input) REAL
60 * The threshold value for the test ratios. A result is
61 * included in the output file if RESULT >= THRESH. To have
62 * every test ratio printed, use THRESH = 0.
63 *
64 * TSTERR (input) LOGICAL
65 * Flag that indicates whether error exits are to be tested.
66 *
67 * A (workspace) REAL array, dimension (LA)
68 *
69 * LA (input) INTEGER
70 * The length of the array A. LA >= (KLMAX+KUMAX+1)*NMAX
71 * where KLMAX is the largest entry in the local array KLVAL,
72 * KUMAX is the largest entry in the local array KUVAL and
73 * NMAX is the largest entry in the input array NVAL.
74 *
75 * AFAC (workspace) REAL array, dimension (LAFAC)
76 *
77 * LAFAC (input) INTEGER
78 * The length of the array AFAC. LAFAC >= (2*KLMAX+KUMAX+1)*NMAX
79 * where KLMAX is the largest entry in the local array KLVAL,
80 * KUMAX is the largest entry in the local array KUVAL and
81 * NMAX is the largest entry in the input array NVAL.
82 *
83 * B (workspace) REAL array, dimension (NMAX*NSMAX)
84 * where NSMAX is the largest entry in NSVAL.
85 *
86 * X (workspace) REAL array, dimension (NMAX*NSMAX)
87 *
88 * XACT (workspace) REAL array, dimension (NMAX*NSMAX)
89 *
90 * WORK (workspace) REAL array, dimension
91 * (NMAX*max(3,NSMAX,NMAX))
92 *
93 * RWORK (workspace) REAL array, dimension
94 * (max(NMAX,2*NSMAX))
95 *
96 * IWORK (workspace) INTEGER array, dimension (2*NMAX)
97 *
98 * NOUT (input) INTEGER
99 * The unit number for output.
100 *
101 * =====================================================================
102 *
103 * .. Parameters ..
104 REAL ONE, ZERO
105 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
106 INTEGER NTYPES, NTESTS
107 PARAMETER ( NTYPES = 8, NTESTS = 7 )
108 INTEGER NBW, NTRAN
109 PARAMETER ( NBW = 4, NTRAN = 3 )
110 * ..
111 * .. Local Scalars ..
112 LOGICAL TRFCON, ZEROT
113 CHARACTER DIST, NORM, TRANS, TYPE, XTYPE
114 CHARACTER*3 PATH
115 INTEGER I, I1, I2, IKL, IKU, IM, IMAT, IN, INB, INFO,
116 $ IOFF, IRHS, ITRAN, IZERO, J, K, KL, KOFF, KU,
117 $ LDA, LDAFAC, LDB, M, MODE, N, NB, NERRS, NFAIL,
118 $ NIMAT, NKL, NKU, NRHS, NRUN
119 REAL AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, RCOND,
120 $ RCONDC, RCONDI, RCONDO
121 * ..
122 * .. Local Arrays ..
123 CHARACTER TRANSS( NTRAN )
124 INTEGER ISEED( 4 ), ISEEDY( 4 ), KLVAL( NBW ),
125 $ KUVAL( NBW )
126 REAL RESULT( NTESTS )
127 * ..
128 * .. External Functions ..
129 REAL SGET06, SLANGB, SLANGE
130 EXTERNAL SGET06, SLANGB, SLANGE
131 * ..
132 * .. External Subroutines ..
133 EXTERNAL ALAERH, ALAHD, ALASUM, SCOPY, SERRGE, SGBCON,
134 $ SGBRFS, SGBT01, SGBT02, SGBT05, SGBTRF, SGBTRS,
135 $ SGET04, SLACPY, SLARHS, SLASET, SLATB4, SLATMS,
136 $ XLAENV
137 * ..
138 * .. Intrinsic Functions ..
139 INTRINSIC MAX, MIN
140 * ..
141 * .. Scalars in Common ..
142 LOGICAL LERR, OK
143 CHARACTER*32 SRNAMT
144 INTEGER INFOT, NUNIT
145 * ..
146 * .. Common blocks ..
147 COMMON / INFOC / INFOT, NUNIT, OK, LERR
148 COMMON / SRNAMC / SRNAMT
149 * ..
150 * .. Data statements ..
151 DATA ISEEDY / 1988, 1989, 1990, 1991 / ,
152 $ TRANSS / 'N', 'T', 'C' /
153 * ..
154 * .. Executable Statements ..
155 *
156 * Initialize constants and the random number seed.
157 *
158 PATH( 1: 1 ) = 'Single precision'
159 PATH( 2: 3 ) = 'GB'
160 NRUN = 0
161 NFAIL = 0
162 NERRS = 0
163 DO 10 I = 1, 4
164 ISEED( I ) = ISEEDY( I )
165 10 CONTINUE
166 *
167 * Test the error exits
168 *
169 IF( TSTERR )
170 $ CALL SERRGE( PATH, NOUT )
171 INFOT = 0
172 CALL XLAENV( 2, 2 )
173 *
174 * Initialize the first value for the lower and upper bandwidths.
175 *
176 KLVAL( 1 ) = 0
177 KUVAL( 1 ) = 0
178 *
179 * Do for each value of M in MVAL
180 *
181 DO 160 IM = 1, NM
182 M = MVAL( IM )
183 *
184 * Set values to use for the lower bandwidth.
185 *
186 KLVAL( 2 ) = M + ( M+1 ) / 4
187 *
188 * KLVAL( 2 ) = MAX( M-1, 0 )
189 *
190 KLVAL( 3 ) = ( 3*M-1 ) / 4
191 KLVAL( 4 ) = ( M+1 ) / 4
192 *
193 * Do for each value of N in NVAL
194 *
195 DO 150 IN = 1, NN
196 N = NVAL( IN )
197 XTYPE = 'N'
198 *
199 * Set values to use for the upper bandwidth.
200 *
201 KUVAL( 2 ) = N + ( N+1 ) / 4
202 *
203 * KUVAL( 2 ) = MAX( N-1, 0 )
204 *
205 KUVAL( 3 ) = ( 3*N-1 ) / 4
206 KUVAL( 4 ) = ( N+1 ) / 4
207 *
208 * Set limits on the number of loop iterations.
209 *
210 NKL = MIN( M+1, 4 )
211 IF( N.EQ.0 )
212 $ NKL = 2
213 NKU = MIN( N+1, 4 )
214 IF( M.EQ.0 )
215 $ NKU = 2
216 NIMAT = NTYPES
217 IF( M.LE.0 .OR. N.LE.0 )
218 $ NIMAT = 1
219 *
220 DO 140 IKL = 1, NKL
221 *
222 * Do for KL = 0, (5*M+1)/4, (3M-1)/4, and (M+1)/4. This
223 * order makes it easier to skip redundant values for small
224 * values of M.
225 *
226 KL = KLVAL( IKL )
227 DO 130 IKU = 1, NKU
228 *
229 * Do for KU = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This
230 * order makes it easier to skip redundant values for
231 * small values of N.
232 *
233 KU = KUVAL( IKU )
234 *
235 * Check that A and AFAC are big enough to generate this
236 * matrix.
237 *
238 LDA = KL + KU + 1
239 LDAFAC = 2*KL + KU + 1
240 IF( ( LDA*N ).GT.LA .OR. ( LDAFAC*N ).GT.LAFAC ) THEN
241 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
242 $ CALL ALAHD( NOUT, PATH )
243 IF( N*( KL+KU+1 ).GT.LA ) THEN
244 WRITE( NOUT, FMT = 9999 )LA, M, N, KL, KU,
245 $ N*( KL+KU+1 )
246 NERRS = NERRS + 1
247 END IF
248 IF( N*( 2*KL+KU+1 ).GT.LAFAC ) THEN
249 WRITE( NOUT, FMT = 9998 )LAFAC, M, N, KL, KU,
250 $ N*( 2*KL+KU+1 )
251 NERRS = NERRS + 1
252 END IF
253 GO TO 130
254 END IF
255 *
256 DO 120 IMAT = 1, NIMAT
257 *
258 * Do the tests only if DOTYPE( IMAT ) is true.
259 *
260 IF( .NOT.DOTYPE( IMAT ) )
261 $ GO TO 120
262 *
263 * Skip types 2, 3, or 4 if the matrix size is too
264 * small.
265 *
266 ZEROT = IMAT.GE.2 .AND. IMAT.LE.4
267 IF( ZEROT .AND. N.LT.IMAT-1 )
268 $ GO TO 120
269 *
270 IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 1 ) ) THEN
271 *
272 * Set up parameters with SLATB4 and generate a
273 * test matrix with SLATMS.
274 *
275 CALL SLATB4( PATH, IMAT, M, N, TYPE, KL, KU,
276 $ ANORM, MODE, CNDNUM, DIST )
277 *
278 KOFF = MAX( 1, KU+2-N )
279 DO 20 I = 1, KOFF - 1
280 A( I ) = ZERO
281 20 CONTINUE
282 SRNAMT = 'SLATMS'
283 CALL SLATMS( M, N, DIST, ISEED, TYPE, RWORK,
284 $ MODE, CNDNUM, ANORM, KL, KU, 'Z',
285 $ A( KOFF ), LDA, WORK, INFO )
286 *
287 * Check the error code from SLATMS.
288 *
289 IF( INFO.NE.0 ) THEN
290 CALL ALAERH( PATH, 'SLATMS', INFO, 0, ' ', M,
291 $ N, KL, KU, -1, IMAT, NFAIL,
292 $ NERRS, NOUT )
293 GO TO 120
294 END IF
295 ELSE IF( IZERO.GT.0 ) THEN
296 *
297 * Use the same matrix for types 3 and 4 as for
298 * type 2 by copying back the zeroed out column.
299 *
300 CALL SCOPY( I2-I1+1, B, 1, A( IOFF+I1 ), 1 )
301 END IF
302 *
303 * For types 2, 3, and 4, zero one or more columns of
304 * the matrix to test that INFO is returned correctly.
305 *
306 IZERO = 0
307 IF( ZEROT ) THEN
308 IF( IMAT.EQ.2 ) THEN
309 IZERO = 1
310 ELSE IF( IMAT.EQ.3 ) THEN
311 IZERO = MIN( M, N )
312 ELSE
313 IZERO = MIN( M, N ) / 2 + 1
314 END IF
315 IOFF = ( IZERO-1 )*LDA
316 IF( IMAT.LT.4 ) THEN
317 *
318 * Store the column to be zeroed out in B.
319 *
320 I1 = MAX( 1, KU+2-IZERO )
321 I2 = MIN( KL+KU+1, KU+1+( M-IZERO ) )
322 CALL SCOPY( I2-I1+1, A( IOFF+I1 ), 1, B, 1 )
323 *
324 DO 30 I = I1, I2
325 A( IOFF+I ) = ZERO
326 30 CONTINUE
327 ELSE
328 DO 50 J = IZERO, N
329 DO 40 I = MAX( 1, KU+2-J ),
330 $ MIN( KL+KU+1, KU+1+( M-J ) )
331 A( IOFF+I ) = ZERO
332 40 CONTINUE
333 IOFF = IOFF + LDA
334 50 CONTINUE
335 END IF
336 END IF
337 *
338 * These lines, if used in place of the calls in the
339 * loop over INB, cause the code to bomb on a Sun
340 * SPARCstation.
341 *
342 * ANORMO = SLANGB( 'O', N, KL, KU, A, LDA, RWORK )
343 * ANORMI = SLANGB( 'I', N, KL, KU, A, LDA, RWORK )
344 *
345 * Do for each blocksize in NBVAL
346 *
347 DO 110 INB = 1, NNB
348 NB = NBVAL( INB )
349 CALL XLAENV( 1, NB )
350 *
351 * Compute the LU factorization of the band matrix.
352 *
353 IF( M.GT.0 .AND. N.GT.0 )
354 $ CALL SLACPY( 'Full', KL+KU+1, N, A, LDA,
355 $ AFAC( KL+1 ), LDAFAC )
356 SRNAMT = 'SGBTRF'
357 CALL SGBTRF( M, N, KL, KU, AFAC, LDAFAC, IWORK,
358 $ INFO )
359 *
360 * Check error code from SGBTRF.
361 *
362 IF( INFO.NE.IZERO )
363 $ CALL ALAERH( PATH, 'SGBTRF', INFO, IZERO,
364 $ ' ', M, N, KL, KU, NB, IMAT,
365 $ NFAIL, NERRS, NOUT )
366 TRFCON = .FALSE.
367 *
368 *+ TEST 1
369 * Reconstruct matrix from factors and compute
370 * residual.
371 *
372 CALL SGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC,
373 $ IWORK, WORK, RESULT( 1 ) )
374 *
375 * Print information about the tests so far that
376 * did not pass the threshold.
377 *
378 IF( RESULT( 1 ).GE.THRESH ) THEN
379 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
380 $ CALL ALAHD( NOUT, PATH )
381 WRITE( NOUT, FMT = 9997 )M, N, KL, KU, NB,
382 $ IMAT, 1, RESULT( 1 )
383 NFAIL = NFAIL + 1
384 END IF
385 NRUN = NRUN + 1
386 *
387 * Skip the remaining tests if this is not the
388 * first block size or if M .ne. N.
389 *
390 IF( INB.GT.1 .OR. M.NE.N )
391 $ GO TO 110
392 *
393 ANORMO = SLANGB( 'O', N, KL, KU, A, LDA, RWORK )
394 ANORMI = SLANGB( 'I', N, KL, KU, A, LDA, RWORK )
395 *
396 IF( INFO.EQ.0 ) THEN
397 *
398 * Form the inverse of A so we can get a good
399 * estimate of CNDNUM = norm(A) * norm(inv(A)).
400 *
401 LDB = MAX( 1, N )
402 CALL SLASET( 'Full', N, N, ZERO, ONE, WORK,
403 $ LDB )
404 SRNAMT = 'SGBTRS'
405 CALL SGBTRS( 'No transpose', N, KL, KU, N,
406 $ AFAC, LDAFAC, IWORK, WORK, LDB,
407 $ INFO )
408 *
409 * Compute the 1-norm condition number of A.
410 *
411 AINVNM = SLANGE( 'O', N, N, WORK, LDB,
412 $ RWORK )
413 IF( ANORMO.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
414 RCONDO = ONE
415 ELSE
416 RCONDO = ( ONE / ANORMO ) / AINVNM
417 END IF
418 *
419 * Compute the infinity-norm condition number of
420 * A.
421 *
422 AINVNM = SLANGE( 'I', N, N, WORK, LDB,
423 $ RWORK )
424 IF( ANORMI.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
425 RCONDI = ONE
426 ELSE
427 RCONDI = ( ONE / ANORMI ) / AINVNM
428 END IF
429 ELSE
430 *
431 * Do only the condition estimate if INFO.NE.0.
432 *
433 TRFCON = .TRUE.
434 RCONDO = ZERO
435 RCONDI = ZERO
436 END IF
437 *
438 * Skip the solve tests if the matrix is singular.
439 *
440 IF( TRFCON )
441 $ GO TO 90
442 *
443 DO 80 IRHS = 1, NNS
444 NRHS = NSVAL( IRHS )
445 XTYPE = 'N'
446 *
447 DO 70 ITRAN = 1, NTRAN
448 TRANS = TRANSS( ITRAN )
449 IF( ITRAN.EQ.1 ) THEN
450 RCONDC = RCONDO
451 NORM = 'O'
452 ELSE
453 RCONDC = RCONDI
454 NORM = 'I'
455 END IF
456 *
457 *+ TEST 2:
458 * Solve and compute residual for A * X = B.
459 *
460 SRNAMT = 'SLARHS'
461 CALL SLARHS( PATH, XTYPE, ' ', TRANS, N,
462 $ N, KL, KU, NRHS, A, LDA,
463 $ XACT, LDB, B, LDB, ISEED,
464 $ INFO )
465 XTYPE = 'C'
466 CALL SLACPY( 'Full', N, NRHS, B, LDB, X,
467 $ LDB )
468 *
469 SRNAMT = 'SGBTRS'
470 CALL SGBTRS( TRANS, N, KL, KU, NRHS, AFAC,
471 $ LDAFAC, IWORK, X, LDB, INFO )
472 *
473 * Check error code from SGBTRS.
474 *
475 IF( INFO.NE.0 )
476 $ CALL ALAERH( PATH, 'SGBTRS', INFO, 0,
477 $ TRANS, N, N, KL, KU, -1,
478 $ IMAT, NFAIL, NERRS, NOUT )
479 *
480 CALL SLACPY( 'Full', N, NRHS, B, LDB,
481 $ WORK, LDB )
482 CALL SGBT02( TRANS, M, N, KL, KU, NRHS, A,
483 $ LDA, X, LDB, WORK, LDB,
484 $ RESULT( 2 ) )
485 *
486 *+ TEST 3:
487 * Check solution from generated exact
488 * solution.
489 *
490 CALL SGET04( N, NRHS, X, LDB, XACT, LDB,
491 $ RCONDC, RESULT( 3 ) )
492 *
493 *+ TESTS 4, 5, 6:
494 * Use iterative refinement to improve the
495 * solution.
496 *
497 SRNAMT = 'SGBRFS'
498 CALL SGBRFS( TRANS, N, KL, KU, NRHS, A,
499 $ LDA, AFAC, LDAFAC, IWORK, B,
500 $ LDB, X, LDB, RWORK,
501 $ RWORK( NRHS+1 ), WORK,
502 $ IWORK( N+1 ), INFO )
503 *
504 * Check error code from SGBRFS.
505 *
506 IF( INFO.NE.0 )
507 $ CALL ALAERH( PATH, 'SGBRFS', INFO, 0,
508 $ TRANS, N, N, KL, KU, NRHS,
509 $ IMAT, NFAIL, NERRS, NOUT )
510 *
511 CALL SGET04( N, NRHS, X, LDB, XACT, LDB,
512 $ RCONDC, RESULT( 4 ) )
513 CALL SGBT05( TRANS, N, KL, KU, NRHS, A,
514 $ LDA, B, LDB, X, LDB, XACT,
515 $ LDB, RWORK, RWORK( NRHS+1 ),
516 $ RESULT( 5 ) )
517 DO 60 K = 2, 6
518 IF( RESULT( K ).GE.THRESH ) THEN
519 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
520 $ CALL ALAHD( NOUT, PATH )
521 WRITE( NOUT, FMT = 9996 )TRANS, N,
522 $ KL, KU, NRHS, IMAT, K,
523 $ RESULT( K )
524 NFAIL = NFAIL + 1
525 END IF
526 60 CONTINUE
527 NRUN = NRUN + 5
528 70 CONTINUE
529 80 CONTINUE
530 *
531 *+ TEST 7:
532 * Get an estimate of RCOND = 1/CNDNUM.
533 *
534 90 CONTINUE
535 DO 100 ITRAN = 1, 2
536 IF( ITRAN.EQ.1 ) THEN
537 ANORM = ANORMO
538 RCONDC = RCONDO
539 NORM = 'O'
540 ELSE
541 ANORM = ANORMI
542 RCONDC = RCONDI
543 NORM = 'I'
544 END IF
545 SRNAMT = 'SGBCON'
546 CALL SGBCON( NORM, N, KL, KU, AFAC, LDAFAC,
547 $ IWORK, ANORM, RCOND, WORK,
548 $ IWORK( N+1 ), INFO )
549 *
550 * Check error code from SGBCON.
551 *
552 IF( INFO.NE.0 )
553 $ CALL ALAERH( PATH, 'SGBCON', INFO, 0,
554 $ NORM, N, N, KL, KU, -1, IMAT,
555 $ NFAIL, NERRS, NOUT )
556 *
557 RESULT( 7 ) = SGET06( RCOND, RCONDC )
558 *
559 * Print information about the tests that did
560 * not pass the threshold.
561 *
562 IF( RESULT( 7 ).GE.THRESH ) THEN
563 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
564 $ CALL ALAHD( NOUT, PATH )
565 WRITE( NOUT, FMT = 9995 )NORM, N, KL, KU,
566 $ IMAT, 7, RESULT( 7 )
567 NFAIL = NFAIL + 1
568 END IF
569 NRUN = NRUN + 1
570 100 CONTINUE
571 *
572 110 CONTINUE
573 120 CONTINUE
574 130 CONTINUE
575 140 CONTINUE
576 150 CONTINUE
577 160 CONTINUE
578 *
579 * Print a summary of the results.
580 *
581 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
582 *
583 9999 FORMAT( ' *** In SCHKGB, LA=', I5, ' is too small for M=', I5,
584 $ ', N=', I5, ', KL=', I4, ', KU=', I4,
585 $ / ' ==> Increase LA to at least ', I5 )
586 9998 FORMAT( ' *** In SCHKGB, LAFAC=', I5, ' is too small for M=', I5,
587 $ ', N=', I5, ', KL=', I4, ', KU=', I4,
588 $ / ' ==> Increase LAFAC to at least ', I5 )
589 9997 FORMAT( ' M =', I5, ', N =', I5, ', KL=', I5, ', KU=', I5,
590 $ ', NB =', I4, ', type ', I1, ', test(', I1, ')=', G12.5 )
591 9996 FORMAT( ' TRANS=''', A1, ''', N=', I5, ', KL=', I5, ', KU=', I5,
592 $ ', NRHS=', I3, ', type ', I1, ', test(', I1, ')=', G12.5 )
593 9995 FORMAT( ' NORM =''', A1, ''', N=', I5, ', KL=', I5, ', KU=', I5,
594 $ ',', 10X, ' type ', I1, ', test(', I1, ')=', G12.5 )
595 *
596 RETURN
597 *
598 * End of SCHKGB
599 *
600 END
2 $ NSVAL, THRESH, TSTERR, A, LA, AFAC, LAFAC, B,
3 $ X, XACT, WORK, 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 LA, LAFAC, NM, NN, NNB, NNS, NOUT
12 REAL THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL DOTYPE( * )
16 INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
17 $ NVAL( * )
18 REAL A( * ), AFAC( * ), B( * ), RWORK( * ),
19 $ WORK( * ), X( * ), XACT( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * SCHKGB tests SGBTRF, -TRS, -RFS, and -CON
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 * NM (input) INTEGER
36 * The number of values of M contained in the vector MVAL.
37 *
38 * MVAL (input) INTEGER array, dimension (NM)
39 * The values of the matrix row dimension M.
40 *
41 * NN (input) INTEGER
42 * The number of values of N contained in the vector NVAL.
43 *
44 * NVAL (input) INTEGER array, dimension (NN)
45 * The values of the matrix column dimension N.
46 *
47 * NNB (input) INTEGER
48 * The number of values of NB contained in the vector NBVAL.
49 *
50 * NBVAL (input) INTEGER array, dimension (NNB)
51 * The values of the blocksize NB.
52 *
53 * NNS (input) INTEGER
54 * The number of values of NRHS contained in the vector NSVAL.
55 *
56 * NSVAL (input) INTEGER array, dimension (NNS)
57 * The values of the number of right hand sides NRHS.
58 *
59 * THRESH (input) REAL
60 * The threshold value for the test ratios. A result is
61 * included in the output file if RESULT >= THRESH. To have
62 * every test ratio printed, use THRESH = 0.
63 *
64 * TSTERR (input) LOGICAL
65 * Flag that indicates whether error exits are to be tested.
66 *
67 * A (workspace) REAL array, dimension (LA)
68 *
69 * LA (input) INTEGER
70 * The length of the array A. LA >= (KLMAX+KUMAX+1)*NMAX
71 * where KLMAX is the largest entry in the local array KLVAL,
72 * KUMAX is the largest entry in the local array KUVAL and
73 * NMAX is the largest entry in the input array NVAL.
74 *
75 * AFAC (workspace) REAL array, dimension (LAFAC)
76 *
77 * LAFAC (input) INTEGER
78 * The length of the array AFAC. LAFAC >= (2*KLMAX+KUMAX+1)*NMAX
79 * where KLMAX is the largest entry in the local array KLVAL,
80 * KUMAX is the largest entry in the local array KUVAL and
81 * NMAX is the largest entry in the input array NVAL.
82 *
83 * B (workspace) REAL array, dimension (NMAX*NSMAX)
84 * where NSMAX is the largest entry in NSVAL.
85 *
86 * X (workspace) REAL array, dimension (NMAX*NSMAX)
87 *
88 * XACT (workspace) REAL array, dimension (NMAX*NSMAX)
89 *
90 * WORK (workspace) REAL array, dimension
91 * (NMAX*max(3,NSMAX,NMAX))
92 *
93 * RWORK (workspace) REAL array, dimension
94 * (max(NMAX,2*NSMAX))
95 *
96 * IWORK (workspace) INTEGER array, dimension (2*NMAX)
97 *
98 * NOUT (input) INTEGER
99 * The unit number for output.
100 *
101 * =====================================================================
102 *
103 * .. Parameters ..
104 REAL ONE, ZERO
105 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
106 INTEGER NTYPES, NTESTS
107 PARAMETER ( NTYPES = 8, NTESTS = 7 )
108 INTEGER NBW, NTRAN
109 PARAMETER ( NBW = 4, NTRAN = 3 )
110 * ..
111 * .. Local Scalars ..
112 LOGICAL TRFCON, ZEROT
113 CHARACTER DIST, NORM, TRANS, TYPE, XTYPE
114 CHARACTER*3 PATH
115 INTEGER I, I1, I2, IKL, IKU, IM, IMAT, IN, INB, INFO,
116 $ IOFF, IRHS, ITRAN, IZERO, J, K, KL, KOFF, KU,
117 $ LDA, LDAFAC, LDB, M, MODE, N, NB, NERRS, NFAIL,
118 $ NIMAT, NKL, NKU, NRHS, NRUN
119 REAL AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, RCOND,
120 $ RCONDC, RCONDI, RCONDO
121 * ..
122 * .. Local Arrays ..
123 CHARACTER TRANSS( NTRAN )
124 INTEGER ISEED( 4 ), ISEEDY( 4 ), KLVAL( NBW ),
125 $ KUVAL( NBW )
126 REAL RESULT( NTESTS )
127 * ..
128 * .. External Functions ..
129 REAL SGET06, SLANGB, SLANGE
130 EXTERNAL SGET06, SLANGB, SLANGE
131 * ..
132 * .. External Subroutines ..
133 EXTERNAL ALAERH, ALAHD, ALASUM, SCOPY, SERRGE, SGBCON,
134 $ SGBRFS, SGBT01, SGBT02, SGBT05, SGBTRF, SGBTRS,
135 $ SGET04, SLACPY, SLARHS, SLASET, SLATB4, SLATMS,
136 $ XLAENV
137 * ..
138 * .. Intrinsic Functions ..
139 INTRINSIC MAX, MIN
140 * ..
141 * .. Scalars in Common ..
142 LOGICAL LERR, OK
143 CHARACTER*32 SRNAMT
144 INTEGER INFOT, NUNIT
145 * ..
146 * .. Common blocks ..
147 COMMON / INFOC / INFOT, NUNIT, OK, LERR
148 COMMON / SRNAMC / SRNAMT
149 * ..
150 * .. Data statements ..
151 DATA ISEEDY / 1988, 1989, 1990, 1991 / ,
152 $ TRANSS / 'N', 'T', 'C' /
153 * ..
154 * .. Executable Statements ..
155 *
156 * Initialize constants and the random number seed.
157 *
158 PATH( 1: 1 ) = 'Single precision'
159 PATH( 2: 3 ) = 'GB'
160 NRUN = 0
161 NFAIL = 0
162 NERRS = 0
163 DO 10 I = 1, 4
164 ISEED( I ) = ISEEDY( I )
165 10 CONTINUE
166 *
167 * Test the error exits
168 *
169 IF( TSTERR )
170 $ CALL SERRGE( PATH, NOUT )
171 INFOT = 0
172 CALL XLAENV( 2, 2 )
173 *
174 * Initialize the first value for the lower and upper bandwidths.
175 *
176 KLVAL( 1 ) = 0
177 KUVAL( 1 ) = 0
178 *
179 * Do for each value of M in MVAL
180 *
181 DO 160 IM = 1, NM
182 M = MVAL( IM )
183 *
184 * Set values to use for the lower bandwidth.
185 *
186 KLVAL( 2 ) = M + ( M+1 ) / 4
187 *
188 * KLVAL( 2 ) = MAX( M-1, 0 )
189 *
190 KLVAL( 3 ) = ( 3*M-1 ) / 4
191 KLVAL( 4 ) = ( M+1 ) / 4
192 *
193 * Do for each value of N in NVAL
194 *
195 DO 150 IN = 1, NN
196 N = NVAL( IN )
197 XTYPE = 'N'
198 *
199 * Set values to use for the upper bandwidth.
200 *
201 KUVAL( 2 ) = N + ( N+1 ) / 4
202 *
203 * KUVAL( 2 ) = MAX( N-1, 0 )
204 *
205 KUVAL( 3 ) = ( 3*N-1 ) / 4
206 KUVAL( 4 ) = ( N+1 ) / 4
207 *
208 * Set limits on the number of loop iterations.
209 *
210 NKL = MIN( M+1, 4 )
211 IF( N.EQ.0 )
212 $ NKL = 2
213 NKU = MIN( N+1, 4 )
214 IF( M.EQ.0 )
215 $ NKU = 2
216 NIMAT = NTYPES
217 IF( M.LE.0 .OR. N.LE.0 )
218 $ NIMAT = 1
219 *
220 DO 140 IKL = 1, NKL
221 *
222 * Do for KL = 0, (5*M+1)/4, (3M-1)/4, and (M+1)/4. This
223 * order makes it easier to skip redundant values for small
224 * values of M.
225 *
226 KL = KLVAL( IKL )
227 DO 130 IKU = 1, NKU
228 *
229 * Do for KU = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This
230 * order makes it easier to skip redundant values for
231 * small values of N.
232 *
233 KU = KUVAL( IKU )
234 *
235 * Check that A and AFAC are big enough to generate this
236 * matrix.
237 *
238 LDA = KL + KU + 1
239 LDAFAC = 2*KL + KU + 1
240 IF( ( LDA*N ).GT.LA .OR. ( LDAFAC*N ).GT.LAFAC ) THEN
241 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
242 $ CALL ALAHD( NOUT, PATH )
243 IF( N*( KL+KU+1 ).GT.LA ) THEN
244 WRITE( NOUT, FMT = 9999 )LA, M, N, KL, KU,
245 $ N*( KL+KU+1 )
246 NERRS = NERRS + 1
247 END IF
248 IF( N*( 2*KL+KU+1 ).GT.LAFAC ) THEN
249 WRITE( NOUT, FMT = 9998 )LAFAC, M, N, KL, KU,
250 $ N*( 2*KL+KU+1 )
251 NERRS = NERRS + 1
252 END IF
253 GO TO 130
254 END IF
255 *
256 DO 120 IMAT = 1, NIMAT
257 *
258 * Do the tests only if DOTYPE( IMAT ) is true.
259 *
260 IF( .NOT.DOTYPE( IMAT ) )
261 $ GO TO 120
262 *
263 * Skip types 2, 3, or 4 if the matrix size is too
264 * small.
265 *
266 ZEROT = IMAT.GE.2 .AND. IMAT.LE.4
267 IF( ZEROT .AND. N.LT.IMAT-1 )
268 $ GO TO 120
269 *
270 IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 1 ) ) THEN
271 *
272 * Set up parameters with SLATB4 and generate a
273 * test matrix with SLATMS.
274 *
275 CALL SLATB4( PATH, IMAT, M, N, TYPE, KL, KU,
276 $ ANORM, MODE, CNDNUM, DIST )
277 *
278 KOFF = MAX( 1, KU+2-N )
279 DO 20 I = 1, KOFF - 1
280 A( I ) = ZERO
281 20 CONTINUE
282 SRNAMT = 'SLATMS'
283 CALL SLATMS( M, N, DIST, ISEED, TYPE, RWORK,
284 $ MODE, CNDNUM, ANORM, KL, KU, 'Z',
285 $ A( KOFF ), LDA, WORK, INFO )
286 *
287 * Check the error code from SLATMS.
288 *
289 IF( INFO.NE.0 ) THEN
290 CALL ALAERH( PATH, 'SLATMS', INFO, 0, ' ', M,
291 $ N, KL, KU, -1, IMAT, NFAIL,
292 $ NERRS, NOUT )
293 GO TO 120
294 END IF
295 ELSE IF( IZERO.GT.0 ) THEN
296 *
297 * Use the same matrix for types 3 and 4 as for
298 * type 2 by copying back the zeroed out column.
299 *
300 CALL SCOPY( I2-I1+1, B, 1, A( IOFF+I1 ), 1 )
301 END IF
302 *
303 * For types 2, 3, and 4, zero one or more columns of
304 * the matrix to test that INFO is returned correctly.
305 *
306 IZERO = 0
307 IF( ZEROT ) THEN
308 IF( IMAT.EQ.2 ) THEN
309 IZERO = 1
310 ELSE IF( IMAT.EQ.3 ) THEN
311 IZERO = MIN( M, N )
312 ELSE
313 IZERO = MIN( M, N ) / 2 + 1
314 END IF
315 IOFF = ( IZERO-1 )*LDA
316 IF( IMAT.LT.4 ) THEN
317 *
318 * Store the column to be zeroed out in B.
319 *
320 I1 = MAX( 1, KU+2-IZERO )
321 I2 = MIN( KL+KU+1, KU+1+( M-IZERO ) )
322 CALL SCOPY( I2-I1+1, A( IOFF+I1 ), 1, B, 1 )
323 *
324 DO 30 I = I1, I2
325 A( IOFF+I ) = ZERO
326 30 CONTINUE
327 ELSE
328 DO 50 J = IZERO, N
329 DO 40 I = MAX( 1, KU+2-J ),
330 $ MIN( KL+KU+1, KU+1+( M-J ) )
331 A( IOFF+I ) = ZERO
332 40 CONTINUE
333 IOFF = IOFF + LDA
334 50 CONTINUE
335 END IF
336 END IF
337 *
338 * These lines, if used in place of the calls in the
339 * loop over INB, cause the code to bomb on a Sun
340 * SPARCstation.
341 *
342 * ANORMO = SLANGB( 'O', N, KL, KU, A, LDA, RWORK )
343 * ANORMI = SLANGB( 'I', N, KL, KU, A, LDA, RWORK )
344 *
345 * Do for each blocksize in NBVAL
346 *
347 DO 110 INB = 1, NNB
348 NB = NBVAL( INB )
349 CALL XLAENV( 1, NB )
350 *
351 * Compute the LU factorization of the band matrix.
352 *
353 IF( M.GT.0 .AND. N.GT.0 )
354 $ CALL SLACPY( 'Full', KL+KU+1, N, A, LDA,
355 $ AFAC( KL+1 ), LDAFAC )
356 SRNAMT = 'SGBTRF'
357 CALL SGBTRF( M, N, KL, KU, AFAC, LDAFAC, IWORK,
358 $ INFO )
359 *
360 * Check error code from SGBTRF.
361 *
362 IF( INFO.NE.IZERO )
363 $ CALL ALAERH( PATH, 'SGBTRF', INFO, IZERO,
364 $ ' ', M, N, KL, KU, NB, IMAT,
365 $ NFAIL, NERRS, NOUT )
366 TRFCON = .FALSE.
367 *
368 *+ TEST 1
369 * Reconstruct matrix from factors and compute
370 * residual.
371 *
372 CALL SGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC,
373 $ IWORK, WORK, RESULT( 1 ) )
374 *
375 * Print information about the tests so far that
376 * did not pass the threshold.
377 *
378 IF( RESULT( 1 ).GE.THRESH ) THEN
379 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
380 $ CALL ALAHD( NOUT, PATH )
381 WRITE( NOUT, FMT = 9997 )M, N, KL, KU, NB,
382 $ IMAT, 1, RESULT( 1 )
383 NFAIL = NFAIL + 1
384 END IF
385 NRUN = NRUN + 1
386 *
387 * Skip the remaining tests if this is not the
388 * first block size or if M .ne. N.
389 *
390 IF( INB.GT.1 .OR. M.NE.N )
391 $ GO TO 110
392 *
393 ANORMO = SLANGB( 'O', N, KL, KU, A, LDA, RWORK )
394 ANORMI = SLANGB( 'I', N, KL, KU, A, LDA, RWORK )
395 *
396 IF( INFO.EQ.0 ) THEN
397 *
398 * Form the inverse of A so we can get a good
399 * estimate of CNDNUM = norm(A) * norm(inv(A)).
400 *
401 LDB = MAX( 1, N )
402 CALL SLASET( 'Full', N, N, ZERO, ONE, WORK,
403 $ LDB )
404 SRNAMT = 'SGBTRS'
405 CALL SGBTRS( 'No transpose', N, KL, KU, N,
406 $ AFAC, LDAFAC, IWORK, WORK, LDB,
407 $ INFO )
408 *
409 * Compute the 1-norm condition number of A.
410 *
411 AINVNM = SLANGE( 'O', N, N, WORK, LDB,
412 $ RWORK )
413 IF( ANORMO.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
414 RCONDO = ONE
415 ELSE
416 RCONDO = ( ONE / ANORMO ) / AINVNM
417 END IF
418 *
419 * Compute the infinity-norm condition number of
420 * A.
421 *
422 AINVNM = SLANGE( 'I', N, N, WORK, LDB,
423 $ RWORK )
424 IF( ANORMI.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
425 RCONDI = ONE
426 ELSE
427 RCONDI = ( ONE / ANORMI ) / AINVNM
428 END IF
429 ELSE
430 *
431 * Do only the condition estimate if INFO.NE.0.
432 *
433 TRFCON = .TRUE.
434 RCONDO = ZERO
435 RCONDI = ZERO
436 END IF
437 *
438 * Skip the solve tests if the matrix is singular.
439 *
440 IF( TRFCON )
441 $ GO TO 90
442 *
443 DO 80 IRHS = 1, NNS
444 NRHS = NSVAL( IRHS )
445 XTYPE = 'N'
446 *
447 DO 70 ITRAN = 1, NTRAN
448 TRANS = TRANSS( ITRAN )
449 IF( ITRAN.EQ.1 ) THEN
450 RCONDC = RCONDO
451 NORM = 'O'
452 ELSE
453 RCONDC = RCONDI
454 NORM = 'I'
455 END IF
456 *
457 *+ TEST 2:
458 * Solve and compute residual for A * X = B.
459 *
460 SRNAMT = 'SLARHS'
461 CALL SLARHS( PATH, XTYPE, ' ', TRANS, N,
462 $ N, KL, KU, NRHS, A, LDA,
463 $ XACT, LDB, B, LDB, ISEED,
464 $ INFO )
465 XTYPE = 'C'
466 CALL SLACPY( 'Full', N, NRHS, B, LDB, X,
467 $ LDB )
468 *
469 SRNAMT = 'SGBTRS'
470 CALL SGBTRS( TRANS, N, KL, KU, NRHS, AFAC,
471 $ LDAFAC, IWORK, X, LDB, INFO )
472 *
473 * Check error code from SGBTRS.
474 *
475 IF( INFO.NE.0 )
476 $ CALL ALAERH( PATH, 'SGBTRS', INFO, 0,
477 $ TRANS, N, N, KL, KU, -1,
478 $ IMAT, NFAIL, NERRS, NOUT )
479 *
480 CALL SLACPY( 'Full', N, NRHS, B, LDB,
481 $ WORK, LDB )
482 CALL SGBT02( TRANS, M, N, KL, KU, NRHS, A,
483 $ LDA, X, LDB, WORK, LDB,
484 $ RESULT( 2 ) )
485 *
486 *+ TEST 3:
487 * Check solution from generated exact
488 * solution.
489 *
490 CALL SGET04( N, NRHS, X, LDB, XACT, LDB,
491 $ RCONDC, RESULT( 3 ) )
492 *
493 *+ TESTS 4, 5, 6:
494 * Use iterative refinement to improve the
495 * solution.
496 *
497 SRNAMT = 'SGBRFS'
498 CALL SGBRFS( TRANS, N, KL, KU, NRHS, A,
499 $ LDA, AFAC, LDAFAC, IWORK, B,
500 $ LDB, X, LDB, RWORK,
501 $ RWORK( NRHS+1 ), WORK,
502 $ IWORK( N+1 ), INFO )
503 *
504 * Check error code from SGBRFS.
505 *
506 IF( INFO.NE.0 )
507 $ CALL ALAERH( PATH, 'SGBRFS', INFO, 0,
508 $ TRANS, N, N, KL, KU, NRHS,
509 $ IMAT, NFAIL, NERRS, NOUT )
510 *
511 CALL SGET04( N, NRHS, X, LDB, XACT, LDB,
512 $ RCONDC, RESULT( 4 ) )
513 CALL SGBT05( TRANS, N, KL, KU, NRHS, A,
514 $ LDA, B, LDB, X, LDB, XACT,
515 $ LDB, RWORK, RWORK( NRHS+1 ),
516 $ RESULT( 5 ) )
517 DO 60 K = 2, 6
518 IF( RESULT( K ).GE.THRESH ) THEN
519 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
520 $ CALL ALAHD( NOUT, PATH )
521 WRITE( NOUT, FMT = 9996 )TRANS, N,
522 $ KL, KU, NRHS, IMAT, K,
523 $ RESULT( K )
524 NFAIL = NFAIL + 1
525 END IF
526 60 CONTINUE
527 NRUN = NRUN + 5
528 70 CONTINUE
529 80 CONTINUE
530 *
531 *+ TEST 7:
532 * Get an estimate of RCOND = 1/CNDNUM.
533 *
534 90 CONTINUE
535 DO 100 ITRAN = 1, 2
536 IF( ITRAN.EQ.1 ) THEN
537 ANORM = ANORMO
538 RCONDC = RCONDO
539 NORM = 'O'
540 ELSE
541 ANORM = ANORMI
542 RCONDC = RCONDI
543 NORM = 'I'
544 END IF
545 SRNAMT = 'SGBCON'
546 CALL SGBCON( NORM, N, KL, KU, AFAC, LDAFAC,
547 $ IWORK, ANORM, RCOND, WORK,
548 $ IWORK( N+1 ), INFO )
549 *
550 * Check error code from SGBCON.
551 *
552 IF( INFO.NE.0 )
553 $ CALL ALAERH( PATH, 'SGBCON', INFO, 0,
554 $ NORM, N, N, KL, KU, -1, IMAT,
555 $ NFAIL, NERRS, NOUT )
556 *
557 RESULT( 7 ) = SGET06( RCOND, RCONDC )
558 *
559 * Print information about the tests that did
560 * not pass the threshold.
561 *
562 IF( RESULT( 7 ).GE.THRESH ) THEN
563 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
564 $ CALL ALAHD( NOUT, PATH )
565 WRITE( NOUT, FMT = 9995 )NORM, N, KL, KU,
566 $ IMAT, 7, RESULT( 7 )
567 NFAIL = NFAIL + 1
568 END IF
569 NRUN = NRUN + 1
570 100 CONTINUE
571 *
572 110 CONTINUE
573 120 CONTINUE
574 130 CONTINUE
575 140 CONTINUE
576 150 CONTINUE
577 160 CONTINUE
578 *
579 * Print a summary of the results.
580 *
581 CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
582 *
583 9999 FORMAT( ' *** In SCHKGB, LA=', I5, ' is too small for M=', I5,
584 $ ', N=', I5, ', KL=', I4, ', KU=', I4,
585 $ / ' ==> Increase LA to at least ', I5 )
586 9998 FORMAT( ' *** In SCHKGB, LAFAC=', I5, ' is too small for M=', I5,
587 $ ', N=', I5, ', KL=', I4, ', KU=', I4,
588 $ / ' ==> Increase LAFAC to at least ', I5 )
589 9997 FORMAT( ' M =', I5, ', N =', I5, ', KL=', I5, ', KU=', I5,
590 $ ', NB =', I4, ', type ', I1, ', test(', I1, ')=', G12.5 )
591 9996 FORMAT( ' TRANS=''', A1, ''', N=', I5, ', KL=', I5, ', KU=', I5,
592 $ ', NRHS=', I3, ', type ', I1, ', test(', I1, ')=', G12.5 )
593 9995 FORMAT( ' NORM =''', A1, ''', N=', I5, ', KL=', I5, ', KU=', I5,
594 $ ',', 10X, ' type ', I1, ', test(', I1, ')=', G12.5 )
595 *
596 RETURN
597 *
598 * End of SCHKGB
599 *
600 END