1       SUBROUTINE SCHKPB( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
  2      $                   THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
  3      $                   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            NMAX, NN, NNB, NNS, NOUT
 12       REAL               THRESH
 13 *     ..
 14 *     .. Array Arguments ..
 15       LOGICAL            DOTYPE( * )
 16       INTEGER            IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
 17       REAL               A( * ), AFAC( * ), AINV( * ), B( * ),
 18      $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
 19 *     ..
 20 *
 21 *  Purpose
 22 *  =======
 23 *
 24 *  SCHKPB tests SPBTRF, -TRS, -RFS, and -CON.
 25 *
 26 *  Arguments
 27 *  =========
 28 *
 29 *  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
 30 *          The matrix types to be used for testing.  Matrices of type j
 31 *          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
 32 *          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
 33 *
 34 *  NN      (input) INTEGER
 35 *          The number of values of N contained in the vector NVAL.
 36 *
 37 *  NVAL    (input) INTEGER array, dimension (NN)
 38 *          The values of the matrix dimension N.
 39 *
 40 *  NNB     (input) INTEGER
 41 *          The number of values of NB contained in the vector NBVAL.
 42 *
 43 *  NBVAL   (input) INTEGER array, dimension (NBVAL)
 44 *          The values of the blocksize NB.
 45 *
 46 *  NNS     (input) INTEGER
 47 *          The number of values of NRHS contained in the vector NSVAL.
 48 *
 49 *  NSVAL   (input) INTEGER array, dimension (NNS)
 50 *          The values of the number of right hand sides NRHS.
 51 *
 52 *  THRESH  (input) REAL
 53 *          The threshold value for the test ratios.  A result is
 54 *          included in the output file if RESULT >= THRESH.  To have
 55 *          every test ratio printed, use THRESH = 0.
 56 *
 57 *  TSTERR  (input) LOGICAL
 58 *          Flag that indicates whether error exits are to be tested.
 59 *
 60 *  NMAX    (input) INTEGER
 61 *          The maximum value permitted for N, used in dimensioning the
 62 *          work arrays.
 63 *
 64 *  A       (workspace) REAL array, dimension (NMAX*NMAX)
 65 *
 66 *  AFAC    (workspace) REAL array, dimension (NMAX*NMAX)
 67 *
 68 *  AINV    (workspace) REAL array, dimension (NMAX*NMAX)
 69 *
 70 *  B       (workspace) REAL array, dimension (NMAX*NSMAX)
 71 *          where NSMAX is the largest entry in NSVAL.
 72 *
 73 *  X       (workspace) REAL array, dimension (NMAX*NSMAX)
 74 *
 75 *  XACT    (workspace) REAL array, dimension (NMAX*NSMAX)
 76 *
 77 *  WORK    (workspace) REAL array, dimension
 78 *                      (NMAX*max(3,NSMAX))
 79 *
 80 *  RWORK   (workspace) REAL array, dimension
 81 *                      (max(NMAX,2*NSMAX))
 82 *
 83 *  IWORK   (workspace) INTEGER array, dimension (NMAX)
 84 *
 85 *  NOUT    (input) INTEGER
 86 *          The unit number for output.
 87 *
 88 *  =====================================================================
 89 *
 90 *     .. Parameters ..
 91       REAL               ONE, ZERO
 92       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 93       INTEGER            NTYPES, NTESTS
 94       PARAMETER          ( NTYPES = 8, NTESTS = 7 )
 95       INTEGER            NBW
 96       PARAMETER          ( NBW = 4 )
 97 *     ..
 98 *     .. Local Scalars ..
 99       LOGICAL            ZEROT
100       CHARACTER          DIST, PACKIT, TYPE, UPLO, XTYPE
101       CHARACTER*3        PATH
102       INTEGER            I, I1, I2, IKD, IMAT, IN, INB, INFO, IOFF,
103      $                   IRHS, IUPLO, IW, IZERO, K, KD, KL, KOFF, KU,
104      $                   LDA, LDAB, MODE, N, NB, NERRS, NFAIL, NIMAT,
105      $                   NKD, NRHS, NRUN
106       REAL               AINVNM, ANORM, CNDNUM, RCOND, RCONDC
107 *     ..
108 *     .. Local Arrays ..
109       INTEGER            ISEED( 4 ), ISEEDY( 4 ), KDVAL( NBW )
110       REAL               RESULT( NTESTS )
111 *     ..
112 *     .. External Functions ..
113       REAL               SGET06, SLANGE, SLANSB
114       EXTERNAL           SGET06, SLANGE, SLANSB
115 *     ..
116 *     .. External Subroutines ..
117       EXTERNAL           ALAERH, ALAHD, ALASUM, SCOPY, SERRPO, SGET04,
118      $                   SLACPY, SLARHS, SLASET, SLATB4, SLATMS, SPBCON,
119      $                   SPBRFS, SPBT01, SPBT02, SPBT05, SPBTRF, SPBTRS,
120      $                   SSWAP, XLAENV
121 *     ..
122 *     .. Intrinsic Functions ..
123       INTRINSIC          MAXMIN
124 *     ..
125 *     .. Scalars in Common ..
126       LOGICAL            LERR, OK
127       CHARACTER*32       SRNAMT
128       INTEGER            INFOT, NUNIT
129 *     ..
130 *     .. Common blocks ..
131       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
132       COMMON             / SRNAMC / SRNAMT
133 *     ..
134 *     .. Data statements ..
135       DATA               ISEEDY / 1988198919901991 /
136 *     ..
137 *     .. Executable Statements ..
138 *
139 *     Initialize constants and the random number seed.
140 *
141       PATH( 11 ) = 'Single precision'
142       PATH( 23 ) = 'PB'
143       NRUN = 0
144       NFAIL = 0
145       NERRS = 0
146       DO 10 I = 14
147          ISEED( I ) = ISEEDY( I )
148    10 CONTINUE
149 *
150 *     Test the error exits
151 *
152       IF( TSTERR )
153      $   CALL SERRPO( PATH, NOUT )
154       INFOT = 0
155       CALL XLAENV( 22 )
156       KDVAL( 1 ) = 0
157 *
158 *     Do for each value of N in NVAL
159 *
160       DO 90 IN = 1, NN
161          N = NVAL( IN )
162          LDA = MAX( N, 1 )
163          XTYPE = 'N'
164 *
165 *        Set limits on the number of loop iterations.
166 *
167          NKD = MAX1MIN( N, 4 ) )
168          NIMAT = NTYPES
169          IF( N.EQ.0 )
170      $      NIMAT = 1
171 *
172          KDVAL( 2 ) = N + ( N+1 ) / 4
173          KDVAL( 3 ) = ( 3*N-1 ) / 4
174          KDVAL( 4 ) = ( N+1 ) / 4
175 *
176          DO 80 IKD = 1, NKD
177 *
178 *           Do for KD = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This order
179 *           makes it easier to skip redundant values for small values
180 *           of N.
181 *
182             KD = KDVAL( IKD )
183             LDAB = KD + 1
184 *
185 *           Do first for UPLO = 'U', then for UPLO = 'L'
186 *
187             DO 70 IUPLO = 12
188                KOFF = 1
189                IF( IUPLO.EQ.1 ) THEN
190                   UPLO = 'U'
191                   KOFF = MAX1, KD+2-N )
192                   PACKIT = 'Q'
193                ELSE
194                   UPLO = 'L'
195                   PACKIT = 'B'
196                END IF
197 *
198                DO 60 IMAT = 1, NIMAT
199 *
200 *                 Do the tests only if DOTYPE( IMAT ) is true.
201 *
202                   IF.NOT.DOTYPE( IMAT ) )
203      $               GO TO 60
204 *
205 *                 Skip types 2, 3, or 4 if the matrix size is too small.
206 *
207                   ZEROT = IMAT.GE.2 .AND. IMAT.LE.4
208                   IF( ZEROT .AND. N.LT.IMAT-1 )
209      $               GO TO 60
210 *
211                   IF.NOT.ZEROT .OR. .NOT.DOTYPE( 1 ) ) THEN
212 *
213 *                    Set up parameters with SLATB4 and generate a test
214 *                    matrix with SLATMS.
215 *
216                      CALL SLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM,
217      $                            MODE, CNDNUM, DIST )
218 *
219                      SRNAMT = 'SLATMS'
220                      CALL SLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
221      $                            CNDNUM, ANORM, KD, KD, PACKIT,
222      $                            A( KOFF ), LDAB, WORK, INFO )
223 *
224 *                    Check error code from SLATMS.
225 *
226                      IF( INFO.NE.0 ) THEN
227                         CALL ALAERH( PATH, 'SLATMS', INFO, 0, UPLO, N,
228      $                               N, KD, KD, -1, IMAT, NFAIL, NERRS,
229      $                               NOUT )
230                         GO TO 60
231                      END IF
232                   ELSE IF( IZERO.GT.0 ) THEN
233 *
234 *                    Use the same matrix for types 3 and 4 as for type
235 *                    2 by copying back the zeroed out column,
236 *
237                      IW = 2*LDA + 1
238                      IF( IUPLO.EQ.1 ) THEN
239                         IOFF = ( IZERO-1 )*LDAB + KD + 1
240                         CALL SCOPY( IZERO-I1, WORK( IW ), 1,
241      $                              A( IOFF-IZERO+I1 ), 1 )
242                         IW = IW + IZERO - I1
243                         CALL SCOPY( I2-IZERO+1, WORK( IW ), 1,
244      $                              A( IOFF ), MAX( LDAB-11 ) )
245                      ELSE
246                         IOFF = ( I1-1 )*LDAB + 1
247                         CALL SCOPY( IZERO-I1, WORK( IW ), 1,
248      $                              A( IOFF+IZERO-I1 ),
249      $                              MAX( LDAB-11 ) )
250                         IOFF = ( IZERO-1 )*LDAB + 1
251                         IW = IW + IZERO - I1
252                         CALL SCOPY( I2-IZERO+1, WORK( IW ), 1,
253      $                              A( IOFF ), 1 )
254                      END IF
255                   END IF
256 *
257 *                 For types 2-4, zero one row and column of the matrix
258 *                 to test that INFO is returned correctly.
259 *
260                   IZERO = 0
261                   IF( ZEROT ) THEN
262                      IF( IMAT.EQ.2 ) THEN
263                         IZERO = 1
264                      ELSE IF( IMAT.EQ.3 ) THEN
265                         IZERO = N
266                      ELSE
267                         IZERO = N / 2 + 1
268                      END IF
269 *
270 *                    Save the zeroed out row and column in WORK(*,3)
271 *
272                      IW = 2*LDA
273                      DO 20 I = 1MIN2*KD+1, N )
274                         WORK( IW+I ) = ZERO
275    20                CONTINUE
276                      IW = IW + 1
277                      I1 = MAX( IZERO-KD, 1 )
278                      I2 = MIN( IZERO+KD, N )
279 *
280                      IF( IUPLO.EQ.1 ) THEN
281                         IOFF = ( IZERO-1 )*LDAB + KD + 1
282                         CALL SSWAP( IZERO-I1, A( IOFF-IZERO+I1 ), 1,
283      $                              WORK( IW ), 1 )
284                         IW = IW + IZERO - I1
285                         CALL SSWAP( I2-IZERO+1, A( IOFF ),
286      $                              MAX( LDAB-11 ), WORK( IW ), 1 )
287                      ELSE
288                         IOFF = ( I1-1 )*LDAB + 1
289                         CALL SSWAP( IZERO-I1, A( IOFF+IZERO-I1 ),
290      $                              MAX( LDAB-11 ), WORK( IW ), 1 )
291                         IOFF = ( IZERO-1 )*LDAB + 1
292                         IW = IW + IZERO - I1
293                         CALL SSWAP( I2-IZERO+1, A( IOFF ), 1,
294      $                              WORK( IW ), 1 )
295                      END IF
296                   END IF
297 *
298 *                 Do for each value of NB in NBVAL
299 *
300                   DO 50 INB = 1, NNB
301                      NB = NBVAL( INB )
302                      CALL XLAENV( 1, NB )
303 *
304 *                    Compute the L*L' or U'*U factorization of the band
305 *                    matrix.
306 *
307                      CALL SLACPY( 'Full', KD+1, N, A, LDAB, AFAC, LDAB )
308                      SRNAMT = 'SPBTRF'
309                      CALL SPBTRF( UPLO, N, KD, AFAC, LDAB, INFO )
310 *
311 *                    Check error code from SPBTRF.
312 *
313                      IF( INFO.NE.IZERO ) THEN
314                         CALL ALAERH( PATH, 'SPBTRF', INFO, IZERO, UPLO,
315      $                               N, N, KD, KD, NB, IMAT, NFAIL,
316      $                               NERRS, NOUT )
317                         GO TO 50
318                      END IF
319 *
320 *                    Skip the tests if INFO is not 0.
321 *
322                      IF( INFO.NE.0 )
323      $                  GO TO 50
324 *
325 *+    TEST 1
326 *                    Reconstruct matrix from factors and compute
327 *                    residual.
328 *
329                      CALL SLACPY( 'Full', KD+1, N, AFAC, LDAB, AINV,
330      $                            LDAB )
331                      CALL SPBT01( UPLO, N, KD, A, LDAB, AINV, LDAB,
332      $                            RWORK, RESULT1 ) )
333 *
334 *                    Print the test ratio if it is .GE. THRESH.
335 *
336                      IFRESULT1 ).GE.THRESH ) THEN
337                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
338      $                     CALL ALAHD( NOUT, PATH )
339                         WRITE( NOUT, FMT = 9999 )UPLO, N, KD, NB, IMAT,
340      $                     1RESULT1 )
341                         NFAIL = NFAIL + 1
342                      END IF
343                      NRUN = NRUN + 1
344 *
345 *                    Only do other tests if this is the first blocksize.
346 *
347                      IF( INB.GT.1 )
348      $                  GO TO 50
349 *
350 *                    Form the inverse of A so we can get a good estimate
351 *                    of RCONDC = 1/(norm(A) * norm(inv(A))).
352 *
353                      CALL SLASET( 'Full', N, N, ZERO, ONE, AINV, LDA )
354                      SRNAMT = 'SPBTRS'
355                      CALL SPBTRS( UPLO, N, KD, N, AFAC, LDAB, AINV, LDA,
356      $                            INFO )
357 *
358 *                    Compute RCONDC = 1/(norm(A) * norm(inv(A))).
359 *
360                      ANORM = SLANSB( '1', UPLO, N, KD, A, LDAB, RWORK )
361                      AINVNM = SLANGE( '1', N, N, AINV, LDA, RWORK )
362                      IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
363                         RCONDC = ONE
364                      ELSE
365                         RCONDC = ( ONE / ANORM ) / AINVNM
366                      END IF
367 *
368                      DO 40 IRHS = 1, NNS
369                         NRHS = NSVAL( IRHS )
370 *
371 *+    TEST 2
372 *                    Solve and compute residual for A * X = B.
373 *
374                         SRNAMT = 'SLARHS'
375                         CALL SLARHS( PATH, XTYPE, UPLO, ' ', N, N, KD,
376      $                               KD, NRHS, A, LDAB, XACT, LDA, B,
377      $                               LDA, ISEED, INFO )
378                         CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
379 *
380                         SRNAMT = 'SPBTRS'
381                         CALL SPBTRS( UPLO, N, KD, NRHS, AFAC, LDAB, X,
382      $                               LDA, INFO )
383 *
384 *                    Check error code from SPBTRS.
385 *
386                         IF( INFO.NE.0 )
387      $                     CALL ALAERH( PATH, 'SPBTRS', INFO, 0, UPLO,
388      $                                  N, N, KD, KD, NRHS, IMAT, NFAIL,
389      $                                  NERRS, NOUT )
390 *
391                         CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK,
392      $                               LDA )
393                         CALL SPBT02( UPLO, N, KD, NRHS, A, LDAB, X, LDA,
394      $                               WORK, LDA, RWORK, RESULT2 ) )
395 *
396 *+    TEST 3
397 *                    Check solution from generated exact solution.
398 *
399                         CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
400      $                               RESULT3 ) )
401 *
402 *+    TESTS 4, 5, and 6
403 *                    Use iterative refinement to improve the solution.
404 *
405                         SRNAMT = 'SPBRFS'
406                         CALL SPBRFS( UPLO, N, KD, NRHS, A, LDAB, AFAC,
407      $                               LDAB, B, LDA, X, LDA, RWORK,
408      $                               RWORK( NRHS+1 ), WORK, IWORK,
409      $                               INFO )
410 *
411 *                    Check error code from SPBRFS.
412 *
413                         IF( INFO.NE.0 )
414      $                     CALL ALAERH( PATH, 'SPBRFS', INFO, 0, UPLO,
415      $                                  N, N, KD, KD, NRHS, IMAT, NFAIL,
416      $                                  NERRS, NOUT )
417 *
418                         CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
419      $                               RESULT4 ) )
420                         CALL SPBT05( UPLO, N, KD, NRHS, A, LDAB, B, LDA,
421      $                               X, LDA, XACT, LDA, RWORK,
422      $                               RWORK( NRHS+1 ), RESULT5 ) )
423 *
424 *                       Print information about the tests that did not
425 *                       pass the threshold.
426 *
427                         DO 30 K = 26
428                            IFRESULT( K ).GE.THRESH ) THEN
429                               IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
430      $                           CALL ALAHD( NOUT, PATH )
431                               WRITE( NOUT, FMT = 9998 )UPLO, N, KD,
432      $                           NRHS, IMAT, K, RESULT( K )
433                               NFAIL = NFAIL + 1
434                            END IF
435    30                   CONTINUE
436                         NRUN = NRUN + 5
437    40                CONTINUE
438 *
439 *+    TEST 7
440 *                    Get an estimate of RCOND = 1/CNDNUM.
441 *
442                      SRNAMT = 'SPBCON'
443                      CALL SPBCON( UPLO, N, KD, AFAC, LDAB, ANORM, RCOND,
444      $                            WORK, IWORK, INFO )
445 *
446 *                    Check error code from SPBCON.
447 *
448                      IF( INFO.NE.0 )
449      $                  CALL ALAERH( PATH, 'SPBCON', INFO, 0, UPLO, N,
450      $                               N, KD, KD, -1, IMAT, NFAIL, NERRS,
451      $                               NOUT )
452 *
453                      RESULT7 ) = SGET06( RCOND, RCONDC )
454 *
455 *                    Print the test ratio if it is .GE. THRESH.
456 *
457                      IFRESULT7 ).GE.THRESH ) THEN
458                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
459      $                     CALL ALAHD( NOUT, PATH )
460                         WRITE( NOUT, FMT = 9997 )UPLO, N, KD, IMAT, 7,
461      $                     RESULT7 )
462                         NFAIL = NFAIL + 1
463                      END IF
464                      NRUN = NRUN + 1
465    50             CONTINUE
466    60          CONTINUE
467    70       CONTINUE
468    80    CONTINUE
469    90 CONTINUE
470 *
471 *     Print a summary of the results.
472 *
473       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
474 *
475  9999 FORMAT' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NB=', I4,
476      $      ', type ', I2, ', test ', I2, ', ratio= 'G12.5 )
477  9998 FORMAT' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NRHS=', I3,
478      $      ', type ', I2, ', test(', I2, ') = 'G12.5 )
479  9997 FORMAT' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ','10X,
480      $      ' type ', I2, ', test(', I2, ') = 'G12.5 )
481       RETURN
482 *
483 *     End of SCHKPB
484 *
485       END