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