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 (NBVAL)
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*> \date December 2016
167*
168*> \ingroup single_lin
169*
170*  =====================================================================
171      SUBROUTINE SCHKPB( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
172     $                   THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
173     $                   XACT, WORK, RWORK, IWORK, NOUT )
174*
175*  -- LAPACK test routine (version 3.7.0) --
176*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
177*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178*     December 2016
179*
180*     .. Scalar Arguments ..
181      LOGICAL            TSTERR
182      INTEGER            NMAX, NN, NNB, NNS, NOUT
183      REAL               THRESH
184*     ..
185*     .. Array Arguments ..
186      LOGICAL            DOTYPE( * )
187      INTEGER            IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
188      REAL               A( * ), AFAC( * ), AINV( * ), B( * ),
189     $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
190*     ..
191*
192*  =====================================================================
193*
194*     .. Parameters ..
195      REAL               ONE, ZERO
196      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
197      INTEGER            NTYPES, NTESTS
198      PARAMETER          ( NTYPES = 8, NTESTS = 7 )
199      INTEGER            NBW
200      PARAMETER          ( NBW = 4 )
201*     ..
202*     .. Local Scalars ..
203      LOGICAL            ZEROT
204      CHARACTER          DIST, PACKIT, TYPE, UPLO, XTYPE
205      CHARACTER*3        PATH
206      INTEGER            I, I1, I2, IKD, IMAT, IN, INB, INFO, IOFF,
207     $                   IRHS, IUPLO, IW, IZERO, K, KD, KL, KOFF, KU,
208     $                   LDA, LDAB, MODE, N, NB, NERRS, NFAIL, NIMAT,
209     $                   NKD, NRHS, NRUN
210      REAL               AINVNM, ANORM, CNDNUM, RCOND, RCONDC
211*     ..
212*     .. Local Arrays ..
213      INTEGER            ISEED( 4 ), ISEEDY( 4 ), KDVAL( NBW )
214      REAL               RESULT( NTESTS )
215*     ..
216*     .. External Functions ..
217      REAL               SGET06, SLANGE, SLANSB
218      EXTERNAL           SGET06, SLANGE, SLANSB
219*     ..
220*     .. External Subroutines ..
221      EXTERNAL           ALAERH, ALAHD, ALASUM, SCOPY, SERRPO, SGET04,
222     $                   SLACPY, SLARHS, SLASET, SLATB4, SLATMS, SPBCON,
223     $                   SPBRFS, SPBT01, SPBT02, SPBT05, SPBTRF, SPBTRS,
224     $                   SSWAP, XLAENV
225*     ..
226*     .. Intrinsic Functions ..
227      INTRINSIC          MAX, MIN
228*     ..
229*     .. Scalars in Common ..
230      LOGICAL            LERR, OK
231      CHARACTER*32       SRNAMT
232      INTEGER            INFOT, NUNIT
233*     ..
234*     .. Common blocks ..
235      COMMON             / INFOC / INFOT, NUNIT, OK, LERR
236      COMMON             / SRNAMC / SRNAMT
237*     ..
238*     .. Data statements ..
239      DATA               ISEEDY / 1988, 1989, 1990, 1991 /
240*     ..
241*     .. Executable Statements ..
242*
243*     Initialize constants and the random number seed.
244*
245      PATH( 1: 1 ) = 'Single precision'
246      PATH( 2: 3 ) = 'PB'
247      NRUN = 0
248      NFAIL = 0
249      NERRS = 0
250      DO 10 I = 1, 4
251         ISEED( I ) = ISEEDY( I )
252   10 CONTINUE
253*
254*     Test the error exits
255*
256      IF( TSTERR )
257     $   CALL SERRPO( PATH, NOUT )
258      INFOT = 0
259      CALL XLAENV( 2, 2 )
260      KDVAL( 1 ) = 0
261*
262*     Do for each value of N in NVAL
263*
264      DO 90 IN = 1, NN
265         N = NVAL( IN )
266         LDA = MAX( N, 1 )
267         XTYPE = 'N'
268*
269*        Set limits on the number of loop iterations.
270*
271         NKD = MAX( 1, MIN( N, 4 ) )
272         NIMAT = NTYPES
273         IF( N.EQ.0 )
274     $      NIMAT = 1
275*
276         KDVAL( 2 ) = N + ( N+1 ) / 4
277         KDVAL( 3 ) = ( 3*N-1 ) / 4
278         KDVAL( 4 ) = ( N+1 ) / 4
279*
280         DO 80 IKD = 1, NKD
281*
282*           Do for KD = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This order
283*           makes it easier to skip redundant values for small values
284*           of N.
285*
286            KD = KDVAL( IKD )
287            LDAB = KD + 1
288*
289*           Do first for UPLO = 'U', then for UPLO = 'L'
290*
291            DO 70 IUPLO = 1, 2
292               KOFF = 1
293               IF( IUPLO.EQ.1 ) THEN
294                  UPLO = 'U'
295                  KOFF = MAX( 1, KD+2-N )
296                  PACKIT = 'Q'
297               ELSE
298                  UPLO = 'L'
299                  PACKIT = 'B'
300               END IF
301*
302               DO 60 IMAT = 1, NIMAT
303*
304*                 Do the tests only if DOTYPE( IMAT ) is true.
305*
306                  IF( .NOT.DOTYPE( IMAT ) )
307     $               GO TO 60
308*
309*                 Skip types 2, 3, or 4 if the matrix size is too small.
310*
311                  ZEROT = IMAT.GE.2 .AND. IMAT.LE.4
312                  IF( ZEROT .AND. N.LT.IMAT-1 )
313     $               GO TO 60
314*
315                  IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 1 ) ) THEN
316*
317*                    Set up parameters with SLATB4 and generate a test
318*                    matrix with SLATMS.
319*
320                     CALL SLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM,
321     $                            MODE, CNDNUM, DIST )
322*
323                     SRNAMT = 'SLATMS'
324                     CALL SLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
325     $                            CNDNUM, ANORM, KD, KD, PACKIT,
326     $                            A( KOFF ), LDAB, WORK, INFO )
327*
328*                    Check error code from SLATMS.
329*
330                     IF( INFO.NE.0 ) THEN
331                        CALL ALAERH( PATH, 'SLATMS', INFO, 0, UPLO, N,
332     $                               N, KD, KD, -1, IMAT, NFAIL, NERRS,
333     $                               NOUT )
334                        GO TO 60
335                     END IF
336                  ELSE IF( IZERO.GT.0 ) THEN
337*
338*                    Use the same matrix for types 3 and 4 as for type
339*                    2 by copying back the zeroed out column,
340*
341                     IW = 2*LDA + 1
342                     IF( IUPLO.EQ.1 ) THEN
343                        IOFF = ( IZERO-1 )*LDAB + KD + 1
344                        CALL SCOPY( IZERO-I1, WORK( IW ), 1,
345     $                              A( IOFF-IZERO+I1 ), 1 )
346                        IW = IW + IZERO - I1
347                        CALL SCOPY( I2-IZERO+1, WORK( IW ), 1,
348     $                              A( IOFF ), MAX( LDAB-1, 1 ) )
349                     ELSE
350                        IOFF = ( I1-1 )*LDAB + 1
351                        CALL SCOPY( IZERO-I1, WORK( IW ), 1,
352     $                              A( IOFF+IZERO-I1 ),
353     $                              MAX( LDAB-1, 1 ) )
354                        IOFF = ( IZERO-1 )*LDAB + 1
355                        IW = IW + IZERO - I1
356                        CALL SCOPY( I2-IZERO+1, WORK( IW ), 1,
357     $                              A( IOFF ), 1 )
358                     END IF
359                  END IF
360*
361*                 For types 2-4, zero one row and column of the matrix
362*                 to test that INFO is returned correctly.
363*
364                  IZERO = 0
365                  IF( ZEROT ) THEN
366                     IF( IMAT.EQ.2 ) THEN
367                        IZERO = 1
368                     ELSE IF( IMAT.EQ.3 ) THEN
369                        IZERO = N
370                     ELSE
371                        IZERO = N / 2 + 1
372                     END IF
373*
374*                    Save the zeroed out row and column in WORK(*,3)
375*
376                     IW = 2*LDA
377                     DO 20 I = 1, MIN( 2*KD+1, N )
378                        WORK( IW+I ) = ZERO
379   20                CONTINUE
380                     IW = IW + 1
381                     I1 = MAX( IZERO-KD, 1 )
382                     I2 = MIN( IZERO+KD, N )
383*
384                     IF( IUPLO.EQ.1 ) THEN
385                        IOFF = ( IZERO-1 )*LDAB + KD + 1
386                        CALL SSWAP( IZERO-I1, A( IOFF-IZERO+I1 ), 1,
387     $                              WORK( IW ), 1 )
388                        IW = IW + IZERO - I1
389                        CALL SSWAP( I2-IZERO+1, A( IOFF ),
390     $                              MAX( LDAB-1, 1 ), WORK( IW ), 1 )
391                     ELSE
392                        IOFF = ( I1-1 )*LDAB + 1
393                        CALL SSWAP( IZERO-I1, A( IOFF+IZERO-I1 ),
394     $                              MAX( LDAB-1, 1 ), WORK( IW ), 1 )
395                        IOFF = ( IZERO-1 )*LDAB + 1
396                        IW = IW + IZERO - I1
397                        CALL SSWAP( I2-IZERO+1, A( IOFF ), 1,
398     $                              WORK( IW ), 1 )
399                     END IF
400                  END IF
401*
402*                 Do for each value of NB in NBVAL
403*
404                  DO 50 INB = 1, NNB
405                     NB = NBVAL( INB )
406                     CALL XLAENV( 1, NB )
407*
408*                    Compute the L*L' or U'*U factorization of the band
409*                    matrix.
410*
411                     CALL SLACPY( 'Full', KD+1, N, A, LDAB, AFAC, LDAB )
412                     SRNAMT = 'SPBTRF'
413                     CALL SPBTRF( UPLO, N, KD, AFAC, LDAB, INFO )
414*
415*                    Check error code from SPBTRF.
416*
417                     IF( INFO.NE.IZERO ) THEN
418                        CALL ALAERH( PATH, 'SPBTRF', INFO, IZERO, UPLO,
419     $                               N, N, KD, KD, NB, IMAT, NFAIL,
420     $                               NERRS, NOUT )
421                        GO TO 50
422                     END IF
423*
424*                    Skip the tests if INFO is not 0.
425*
426                     IF( INFO.NE.0 )
427     $                  GO TO 50
428*
429*+    TEST 1
430*                    Reconstruct matrix from factors and compute
431*                    residual.
432*
433                     CALL SLACPY( 'Full', KD+1, N, AFAC, LDAB, AINV,
434     $                            LDAB )
435                     CALL SPBT01( UPLO, N, KD, A, LDAB, AINV, LDAB,
436     $                            RWORK, RESULT( 1 ) )
437*
438*                    Print the test ratio if it is .GE. THRESH.
439*
440                     IF( RESULT( 1 ).GE.THRESH ) THEN
441                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
442     $                     CALL ALAHD( NOUT, PATH )
443                        WRITE( NOUT, FMT = 9999 )UPLO, N, KD, NB, IMAT,
444     $                     1, RESULT( 1 )
445                        NFAIL = NFAIL + 1
446                     END IF
447                     NRUN = NRUN + 1
448*
449*                    Only do other tests if this is the first blocksize.
450*
451                     IF( INB.GT.1 )
452     $                  GO TO 50
453*
454*                    Form the inverse of A so we can get a good estimate
455*                    of RCONDC = 1/(norm(A) * norm(inv(A))).
456*
457                     CALL SLASET( 'Full', N, N, ZERO, ONE, AINV, LDA )
458                     SRNAMT = 'SPBTRS'
459                     CALL SPBTRS( UPLO, N, KD, N, AFAC, LDAB, AINV, LDA,
460     $                            INFO )
461*
462*                    Compute RCONDC = 1/(norm(A) * norm(inv(A))).
463*
464                     ANORM = SLANSB( '1', UPLO, N, KD, A, LDAB, RWORK )
465                     AINVNM = SLANGE( '1', N, N, AINV, LDA, RWORK )
466                     IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
467                        RCONDC = ONE
468                     ELSE
469                        RCONDC = ( ONE / ANORM ) / AINVNM
470                     END IF
471*
472                     DO 40 IRHS = 1, NNS
473                        NRHS = NSVAL( IRHS )
474*
475*+    TEST 2
476*                    Solve and compute residual for A * X = B.
477*
478                        SRNAMT = 'SLARHS'
479                        CALL SLARHS( PATH, XTYPE, UPLO, ' ', N, N, KD,
480     $                               KD, NRHS, A, LDAB, XACT, LDA, B,
481     $                               LDA, ISEED, INFO )
482                        CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
483*
484                        SRNAMT = 'SPBTRS'
485                        CALL SPBTRS( UPLO, N, KD, NRHS, AFAC, LDAB, X,
486     $                               LDA, INFO )
487*
488*                    Check error code from SPBTRS.
489*
490                        IF( INFO.NE.0 )
491     $                     CALL ALAERH( PATH, 'SPBTRS', INFO, 0, UPLO,
492     $                                  N, N, KD, KD, NRHS, IMAT, NFAIL,
493     $                                  NERRS, NOUT )
494*
495                        CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK,
496     $                               LDA )
497                        CALL SPBT02( UPLO, N, KD, NRHS, A, LDAB, X, LDA,
498     $                               WORK, LDA, RWORK, RESULT( 2 ) )
499*
500*+    TEST 3
501*                    Check solution from generated exact solution.
502*
503                        CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
504     $                               RESULT( 3 ) )
505*
506*+    TESTS 4, 5, and 6
507*                    Use iterative refinement to improve the solution.
508*
509                        SRNAMT = 'SPBRFS'
510                        CALL SPBRFS( UPLO, N, KD, NRHS, A, LDAB, AFAC,
511     $                               LDAB, B, LDA, X, LDA, RWORK,
512     $                               RWORK( NRHS+1 ), WORK, IWORK,
513     $                               INFO )
514*
515*                    Check error code from SPBRFS.
516*
517                        IF( INFO.NE.0 )
518     $                     CALL ALAERH( PATH, 'SPBRFS', INFO, 0, UPLO,
519     $                                  N, N, KD, KD, NRHS, IMAT, NFAIL,
520     $                                  NERRS, NOUT )
521*
522                        CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
523     $                               RESULT( 4 ) )
524                        CALL SPBT05( UPLO, N, KD, NRHS, A, LDAB, B, LDA,
525     $                               X, LDA, XACT, LDA, RWORK,
526     $                               RWORK( NRHS+1 ), RESULT( 5 ) )
527*
528*                       Print information about the tests that did not
529*                       pass the threshold.
530*
531                        DO 30 K = 2, 6
532                           IF( RESULT( K ).GE.THRESH ) THEN
533                              IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
534     $                           CALL ALAHD( NOUT, PATH )
535                              WRITE( NOUT, FMT = 9998 )UPLO, N, KD,
536     $                           NRHS, IMAT, K, RESULT( K )
537                              NFAIL = NFAIL + 1
538                           END IF
539   30                   CONTINUE
540                        NRUN = NRUN + 5
541   40                CONTINUE
542*
543*+    TEST 7
544*                    Get an estimate of RCOND = 1/CNDNUM.
545*
546                     SRNAMT = 'SPBCON'
547                     CALL SPBCON( UPLO, N, KD, AFAC, LDAB, ANORM, RCOND,
548     $                            WORK, IWORK, INFO )
549*
550*                    Check error code from SPBCON.
551*
552                     IF( INFO.NE.0 )
553     $                  CALL ALAERH( PATH, 'SPBCON', INFO, 0, UPLO, N,
554     $                               N, KD, KD, -1, IMAT, NFAIL, NERRS,
555     $                               NOUT )
556*
557                     RESULT( 7 ) = SGET06( RCOND, RCONDC )
558*
559*                    Print the test ratio if it is .GE. THRESH.
560*
561                     IF( RESULT( 7 ).GE.THRESH ) THEN
562                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
563     $                     CALL ALAHD( NOUT, PATH )
564                        WRITE( NOUT, FMT = 9997 )UPLO, N, KD, IMAT, 7,
565     $                     RESULT( 7 )
566                        NFAIL = NFAIL + 1
567                     END IF
568                     NRUN = NRUN + 1
569   50             CONTINUE
570   60          CONTINUE
571   70       CONTINUE
572   80    CONTINUE
573   90 CONTINUE
574*
575*     Print a summary of the results.
576*
577      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
578*
579 9999 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NB=', I4,
580     $      ', type ', I2, ', test ', I2, ', ratio= ', G12.5 )
581 9998 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ', NRHS=', I3,
582     $      ', type ', I2, ', test(', I2, ') = ', G12.5 )
583 9997 FORMAT( ' UPLO=''', A1, ''', N=', I5, ', KD=', I5, ',', 10X,
584     $      ' type ', I2, ', test(', I2, ') = ', G12.5 )
585      RETURN
586*
587*     End of SCHKPB
588*
589      END
590