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