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