1*> \brief \b DCHKSB
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 DCHKSB( NSIZES, NN, NWDTHS, KK, NTYPES, DOTYPE, ISEED,
12*                          THRESH, NOUNIT, A, LDA, SD, SE, U, LDU, WORK,
13*                          LWORK, RESULT, INFO )
14*
15*       .. Scalar Arguments ..
16*       INTEGER            INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES,
17*      $                   NWDTHS
18*       DOUBLE PRECISION   THRESH
19*       ..
20*       .. Array Arguments ..
21*       LOGICAL            DOTYPE( * )
22*       INTEGER            ISEED( 4 ), KK( * ), NN( * )
23*       DOUBLE PRECISION   A( LDA, * ), RESULT( * ), SD( * ), SE( * ),
24*      $                   U( LDU, * ), WORK( * )
25*       ..
26*
27*
28*> \par Purpose:
29*  =============
30*>
31*> \verbatim
32*>
33*> DCHKSB tests the reduction of a symmetric band matrix to tridiagonal
34*> form, used with the symmetric eigenvalue problem.
35*>
36*> DSBTRD factors a symmetric band matrix A as  U S U' , where ' means
37*> transpose, S is symmetric tridiagonal, and U is orthogonal.
38*> DSBTRD can use either just the lower or just the upper triangle
39*> of A; DCHKSB checks both cases.
40*>
41*> When DCHKSB is called, a number of matrix "sizes" ("n's"), a number
42*> of bandwidths ("k's"), and a number of matrix "types" are
43*> specified.  For each size ("n"), each bandwidth ("k") less than or
44*> equal to "n", and each type of matrix, one matrix will be generated
45*> and used to test the symmetric banded reduction routine.  For each
46*> matrix, a number of tests will be performed:
47*>
48*> (1)     | A - V S V' | / ( |A| n ulp )  computed by DSBTRD with
49*>                                         UPLO='U'
50*>
51*> (2)     | I - UU' | / ( n ulp )
52*>
53*> (3)     | A - V S V' | / ( |A| n ulp )  computed by DSBTRD with
54*>                                         UPLO='L'
55*>
56*> (4)     | I - UU' | / ( n ulp )
57*>
58*> The "sizes" are specified by an array NN(1:NSIZES); the value of
59*> each element NN(j) specifies one size.
60*> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
61*> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
62*> Currently, the list of possible types is:
63*>
64*> (1)  The zero matrix.
65*> (2)  The identity matrix.
66*>
67*> (3)  A diagonal matrix with evenly spaced entries
68*>      1, ..., ULP  and random signs.
69*>      (ULP = (first number larger than 1) - 1 )
70*> (4)  A diagonal matrix with geometrically spaced entries
71*>      1, ..., ULP  and random signs.
72*> (5)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
73*>      and random signs.
74*>
75*> (6)  Same as (4), but multiplied by SQRT( overflow threshold )
76*> (7)  Same as (4), but multiplied by SQRT( underflow threshold )
77*>
78*> (8)  A matrix of the form  U' D U, where U is orthogonal and
79*>      D has evenly spaced entries 1, ..., ULP with random signs
80*>      on the diagonal.
81*>
82*> (9)  A matrix of the form  U' D U, where U is orthogonal and
83*>      D has geometrically spaced entries 1, ..., ULP with random
84*>      signs on the diagonal.
85*>
86*> (10) A matrix of the form  U' D U, where U is orthogonal and
87*>      D has "clustered" entries 1, ULP,..., ULP with random
88*>      signs on the diagonal.
89*>
90*> (11) Same as (8), but multiplied by SQRT( overflow threshold )
91*> (12) Same as (8), but multiplied by SQRT( underflow threshold )
92*>
93*> (13) Symmetric matrix with random entries chosen from (-1,1).
94*> (14) Same as (13), but multiplied by SQRT( overflow threshold )
95*> (15) Same as (13), but multiplied by SQRT( underflow threshold )
96*> \endverbatim
97*
98*  Arguments:
99*  ==========
100*
101*> \param[in] NSIZES
102*> \verbatim
103*>          NSIZES is INTEGER
104*>          The number of sizes of matrices to use.  If it is zero,
105*>          DCHKSB does nothing.  It must be at least zero.
106*> \endverbatim
107*>
108*> \param[in] NN
109*> \verbatim
110*>          NN is INTEGER array, dimension (NSIZES)
111*>          An array containing the sizes to be used for the matrices.
112*>          Zero values will be skipped.  The values must be at least
113*>          zero.
114*> \endverbatim
115*>
116*> \param[in] NWDTHS
117*> \verbatim
118*>          NWDTHS is INTEGER
119*>          The number of bandwidths to use.  If it is zero,
120*>          DCHKSB does nothing.  It must be at least zero.
121*> \endverbatim
122*>
123*> \param[in] KK
124*> \verbatim
125*>          KK is INTEGER array, dimension (NWDTHS)
126*>          An array containing the bandwidths to be used for the band
127*>          matrices.  The values must be at least zero.
128*> \endverbatim
129*>
130*> \param[in] NTYPES
131*> \verbatim
132*>          NTYPES is INTEGER
133*>          The number of elements in DOTYPE.   If it is zero, DCHKSB
134*>          does nothing.  It must be at least zero.  If it is MAXTYP+1
135*>          and NSIZES is 1, then an additional type, MAXTYP+1 is
136*>          defined, which is to use whatever matrix is in A.  This
137*>          is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
138*>          DOTYPE(MAXTYP+1) is .TRUE. .
139*> \endverbatim
140*>
141*> \param[in] DOTYPE
142*> \verbatim
143*>          DOTYPE is LOGICAL array, dimension (NTYPES)
144*>          If DOTYPE(j) is .TRUE., then for each size in NN a
145*>          matrix of that size and of type j will be generated.
146*>          If NTYPES is smaller than the maximum number of types
147*>          defined (PARAMETER MAXTYP), then types NTYPES+1 through
148*>          MAXTYP will not be generated.  If NTYPES is larger
149*>          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
150*>          will be ignored.
151*> \endverbatim
152*>
153*> \param[in,out] ISEED
154*> \verbatim
155*>          ISEED is INTEGER array, dimension (4)
156*>          On entry ISEED specifies the seed of the random number
157*>          generator. The array elements should be between 0 and 4095;
158*>          if not they will be reduced mod 4096.  Also, ISEED(4) must
159*>          be odd.  The random number generator uses a linear
160*>          congruential sequence limited to small integers, and so
161*>          should produce machine independent random numbers. The
162*>          values of ISEED are changed on exit, and can be used in the
163*>          next call to DCHKSB to continue the same random number
164*>          sequence.
165*> \endverbatim
166*>
167*> \param[in] THRESH
168*> \verbatim
169*>          THRESH is DOUBLE PRECISION
170*>          A test will count as "failed" if the "error", computed as
171*>          described above, exceeds THRESH.  Note that the error
172*>          is scaled to be O(1), so THRESH should be a reasonably
173*>          small multiple of 1, e.g., 10 or 100.  In particular,
174*>          it should not depend on the precision (single vs. double)
175*>          or the size of the matrix.  It must be at least zero.
176*> \endverbatim
177*>
178*> \param[in] NOUNIT
179*> \verbatim
180*>          NOUNIT is INTEGER
181*>          The FORTRAN unit number for printing out error messages
182*>          (e.g., if a routine returns IINFO not equal to 0.)
183*> \endverbatim
184*>
185*> \param[in,out] A
186*> \verbatim
187*>          A is DOUBLE PRECISION array, dimension
188*>                            (LDA, max(NN))
189*>          Used to hold the matrix whose eigenvalues are to be
190*>          computed.
191*> \endverbatim
192*>
193*> \param[in] LDA
194*> \verbatim
195*>          LDA is INTEGER
196*>          The leading dimension of A.  It must be at least 2 (not 1!)
197*>          and at least max( KK )+1.
198*> \endverbatim
199*>
200*> \param[out] SD
201*> \verbatim
202*>          SD is DOUBLE PRECISION array, dimension (max(NN))
203*>          Used to hold the diagonal of the tridiagonal matrix computed
204*>          by DSBTRD.
205*> \endverbatim
206*>
207*> \param[out] SE
208*> \verbatim
209*>          SE is DOUBLE PRECISION array, dimension (max(NN))
210*>          Used to hold the off-diagonal of the tridiagonal matrix
211*>          computed by DSBTRD.
212*> \endverbatim
213*>
214*> \param[out] U
215*> \verbatim
216*>          U is DOUBLE PRECISION array, dimension (LDU, max(NN))
217*>          Used to hold the orthogonal matrix computed by DSBTRD.
218*> \endverbatim
219*>
220*> \param[in] LDU
221*> \verbatim
222*>          LDU is INTEGER
223*>          The leading dimension of U.  It must be at least 1
224*>          and at least max( NN ).
225*> \endverbatim
226*>
227*> \param[out] WORK
228*> \verbatim
229*>          WORK is DOUBLE PRECISION array, dimension (LWORK)
230*> \endverbatim
231*>
232*> \param[in] LWORK
233*> \verbatim
234*>          LWORK is INTEGER
235*>          The number of entries in WORK.  This must be at least
236*>          max( LDA+1, max(NN)+1 )*max(NN).
237*> \endverbatim
238*>
239*> \param[out] RESULT
240*> \verbatim
241*>          RESULT is DOUBLE PRECISION array, dimension (4)
242*>          The values computed by the tests described above.
243*>          The values are currently limited to 1/ulp, to avoid
244*>          overflow.
245*> \endverbatim
246*>
247*> \param[out] INFO
248*> \verbatim
249*>          INFO is INTEGER
250*>          If 0, then everything ran OK.
251*>
252*>-----------------------------------------------------------------------
253*>
254*>       Some Local Variables and Parameters:
255*>       ---- ----- --------- --- ----------
256*>       ZERO, ONE       Real 0 and 1.
257*>       MAXTYP          The number of types defined.
258*>       NTEST           The number of tests performed, or which can
259*>                       be performed so far, for the current matrix.
260*>       NTESTT          The total number of tests performed so far.
261*>       NMAX            Largest value in NN.
262*>       NMATS           The number of matrices generated so far.
263*>       NERRS           The number of tests which have exceeded THRESH
264*>                       so far.
265*>       COND, IMODE     Values to be passed to the matrix generators.
266*>       ANORM           Norm of A; passed to matrix generators.
267*>
268*>       OVFL, UNFL      Overflow and underflow thresholds.
269*>       ULP, ULPINV     Finest relative precision and its inverse.
270*>       RTOVFL, RTUNFL  Square roots of the previous 2 values.
271*>               The following four arrays decode JTYPE:
272*>       KTYPE(j)        The general type (1-10) for type "j".
273*>       KMODE(j)        The MODE value to be passed to the matrix
274*>                       generator for type "j".
275*>       KMAGN(j)        The order of magnitude ( O(1),
276*>                       O(overflow^(1/2) ), O(underflow^(1/2) )
277*> \endverbatim
278*
279*  Authors:
280*  ========
281*
282*> \author Univ. of Tennessee
283*> \author Univ. of California Berkeley
284*> \author Univ. of Colorado Denver
285*> \author NAG Ltd.
286*
287*> \ingroup double_eig
288*
289*  =====================================================================
290      SUBROUTINE DCHKSB( NSIZES, NN, NWDTHS, KK, NTYPES, DOTYPE, ISEED,
291     $                   THRESH, NOUNIT, A, LDA, SD, SE, U, LDU, WORK,
292     $                   LWORK, RESULT, INFO )
293*
294*  -- LAPACK test routine --
295*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
296*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
297*
298*     .. Scalar Arguments ..
299      INTEGER            INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES,
300     $                   NWDTHS
301      DOUBLE PRECISION   THRESH
302*     ..
303*     .. Array Arguments ..
304      LOGICAL            DOTYPE( * )
305      INTEGER            ISEED( 4 ), KK( * ), NN( * )
306      DOUBLE PRECISION   A( LDA, * ), RESULT( * ), SD( * ), SE( * ),
307     $                   U( LDU, * ), WORK( * )
308*     ..
309*
310*  =====================================================================
311*
312*     .. Parameters ..
313      DOUBLE PRECISION   ZERO, ONE, TWO, TEN
314      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
315     $                   TEN = 10.0D0 )
316      DOUBLE PRECISION   HALF
317      PARAMETER          ( HALF = ONE / TWO )
318      INTEGER            MAXTYP
319      PARAMETER          ( MAXTYP = 15 )
320*     ..
321*     .. Local Scalars ..
322      LOGICAL            BADNN, BADNNB
323      INTEGER            I, IINFO, IMODE, ITYPE, J, JC, JCOL, JR, JSIZE,
324     $                   JTYPE, JWIDTH, K, KMAX, MTYPES, N, NERRS,
325     $                   NMATS, NMAX, NTEST, NTESTT
326      DOUBLE PRECISION   ANINV, ANORM, COND, OVFL, RTOVFL, RTUNFL,
327     $                   TEMP1, ULP, ULPINV, UNFL
328*     ..
329*     .. Local Arrays ..
330      INTEGER            IDUMMA( 1 ), IOLDSD( 4 ), KMAGN( MAXTYP ),
331     $                   KMODE( MAXTYP ), KTYPE( MAXTYP )
332*     ..
333*     .. External Functions ..
334      DOUBLE PRECISION   DLAMCH
335      EXTERNAL           DLAMCH
336*     ..
337*     .. External Subroutines ..
338      EXTERNAL           DLACPY, DLASET, DLASUM, DLATMR, DLATMS, DSBT21,
339     $                   DSBTRD, XERBLA
340*     ..
341*     .. Intrinsic Functions ..
342      INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
343*     ..
344*     .. Data statements ..
345      DATA               KTYPE / 1, 2, 5*4, 5*5, 3*8 /
346      DATA               KMAGN / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
347     $                   2, 3 /
348      DATA               KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
349     $                   0, 0 /
350*     ..
351*     .. Executable Statements ..
352*
353*     Check for errors
354*
355      NTESTT = 0
356      INFO = 0
357*
358*     Important constants
359*
360      BADNN = .FALSE.
361      NMAX = 1
362      DO 10 J = 1, NSIZES
363         NMAX = MAX( NMAX, NN( J ) )
364         IF( NN( J ).LT.0 )
365     $      BADNN = .TRUE.
366   10 CONTINUE
367*
368      BADNNB = .FALSE.
369      KMAX = 0
370      DO 20 J = 1, NSIZES
371         KMAX = MAX( KMAX, KK( J ) )
372         IF( KK( J ).LT.0 )
373     $      BADNNB = .TRUE.
374   20 CONTINUE
375      KMAX = MIN( NMAX-1, KMAX )
376*
377*     Check for errors
378*
379      IF( NSIZES.LT.0 ) THEN
380         INFO = -1
381      ELSE IF( BADNN ) THEN
382         INFO = -2
383      ELSE IF( NWDTHS.LT.0 ) THEN
384         INFO = -3
385      ELSE IF( BADNNB ) THEN
386         INFO = -4
387      ELSE IF( NTYPES.LT.0 ) THEN
388         INFO = -5
389      ELSE IF( LDA.LT.KMAX+1 ) THEN
390         INFO = -11
391      ELSE IF( LDU.LT.NMAX ) THEN
392         INFO = -15
393      ELSE IF( ( MAX( LDA, NMAX )+1 )*NMAX.GT.LWORK ) THEN
394         INFO = -17
395      END IF
396*
397      IF( INFO.NE.0 ) THEN
398         CALL XERBLA( 'DCHKSB', -INFO )
399         RETURN
400      END IF
401*
402*     Quick return if possible
403*
404      IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 .OR. NWDTHS.EQ.0 )
405     $   RETURN
406*
407*     More Important constants
408*
409      UNFL = DLAMCH( 'Safe minimum' )
410      OVFL = ONE / UNFL
411      ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
412      ULPINV = ONE / ULP
413      RTUNFL = SQRT( UNFL )
414      RTOVFL = SQRT( OVFL )
415*
416*     Loop over sizes, types
417*
418      NERRS = 0
419      NMATS = 0
420*
421      DO 190 JSIZE = 1, NSIZES
422         N = NN( JSIZE )
423         ANINV = ONE / DBLE( MAX( 1, N ) )
424*
425         DO 180 JWIDTH = 1, NWDTHS
426            K = KK( JWIDTH )
427            IF( K.GT.N )
428     $         GO TO 180
429            K = MAX( 0, MIN( N-1, K ) )
430*
431            IF( NSIZES.NE.1 ) THEN
432               MTYPES = MIN( MAXTYP, NTYPES )
433            ELSE
434               MTYPES = MIN( MAXTYP+1, NTYPES )
435            END IF
436*
437            DO 170 JTYPE = 1, MTYPES
438               IF( .NOT.DOTYPE( JTYPE ) )
439     $            GO TO 170
440               NMATS = NMATS + 1
441               NTEST = 0
442*
443               DO 30 J = 1, 4
444                  IOLDSD( J ) = ISEED( J )
445   30          CONTINUE
446*
447*              Compute "A".
448*              Store as "Upper"; later, we will copy to other format.
449*
450*              Control parameters:
451*
452*                  KMAGN  KMODE        KTYPE
453*              =1  O(1)   clustered 1  zero
454*              =2  large  clustered 2  identity
455*              =3  small  exponential  (none)
456*              =4         arithmetic   diagonal, (w/ eigenvalues)
457*              =5         random log   symmetric, w/ eigenvalues
458*              =6         random       (none)
459*              =7                      random diagonal
460*              =8                      random symmetric
461*              =9                      positive definite
462*              =10                     diagonally dominant tridiagonal
463*
464               IF( MTYPES.GT.MAXTYP )
465     $            GO TO 100
466*
467               ITYPE = KTYPE( JTYPE )
468               IMODE = KMODE( JTYPE )
469*
470*              Compute norm
471*
472               GO TO ( 40, 50, 60 )KMAGN( JTYPE )
473*
474   40          CONTINUE
475               ANORM = ONE
476               GO TO 70
477*
478   50          CONTINUE
479               ANORM = ( RTOVFL*ULP )*ANINV
480               GO TO 70
481*
482   60          CONTINUE
483               ANORM = RTUNFL*N*ULPINV
484               GO TO 70
485*
486   70          CONTINUE
487*
488               CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
489               IINFO = 0
490               IF( JTYPE.LE.15 ) THEN
491                  COND = ULPINV
492               ELSE
493                  COND = ULPINV*ANINV / TEN
494               END IF
495*
496*              Special Matrices -- Identity & Jordan block
497*
498*                 Zero
499*
500               IF( ITYPE.EQ.1 ) THEN
501                  IINFO = 0
502*
503               ELSE IF( ITYPE.EQ.2 ) THEN
504*
505*                 Identity
506*
507                  DO 80 JCOL = 1, N
508                     A( K+1, JCOL ) = ANORM
509   80             CONTINUE
510*
511               ELSE IF( ITYPE.EQ.4 ) THEN
512*
513*                 Diagonal Matrix, [Eigen]values Specified
514*
515                  CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
516     $                         ANORM, 0, 0, 'Q', A( K+1, 1 ), LDA,
517     $                         WORK( N+1 ), IINFO )
518*
519               ELSE IF( ITYPE.EQ.5 ) THEN
520*
521*                 Symmetric, eigenvalues specified
522*
523                  CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
524     $                         ANORM, K, K, 'Q', A, LDA, WORK( N+1 ),
525     $                         IINFO )
526*
527               ELSE IF( ITYPE.EQ.7 ) THEN
528*
529*                 Diagonal, random eigenvalues
530*
531                  CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
532     $                         'T', 'N', WORK( N+1 ), 1, ONE,
533     $                         WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
534     $                         ZERO, ANORM, 'Q', A( K+1, 1 ), LDA,
535     $                         IDUMMA, IINFO )
536*
537               ELSE IF( ITYPE.EQ.8 ) THEN
538*
539*                 Symmetric, random eigenvalues
540*
541                  CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
542     $                         'T', 'N', WORK( N+1 ), 1, ONE,
543     $                         WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, K, K,
544     $                         ZERO, ANORM, 'Q', A, LDA, IDUMMA, IINFO )
545*
546               ELSE IF( ITYPE.EQ.9 ) THEN
547*
548*                 Positive definite, eigenvalues specified.
549*
550                  CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
551     $                         ANORM, K, K, 'Q', A, LDA, WORK( N+1 ),
552     $                         IINFO )
553*
554               ELSE IF( ITYPE.EQ.10 ) THEN
555*
556*                 Positive definite tridiagonal, eigenvalues specified.
557*
558                  IF( N.GT.1 )
559     $               K = MAX( 1, K )
560                  CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
561     $                         ANORM, 1, 1, 'Q', A( K, 1 ), LDA,
562     $                         WORK( N+1 ), IINFO )
563                  DO 90 I = 2, N
564                     TEMP1 = ABS( A( K, I ) ) /
565     $                       SQRT( ABS( A( K+1, I-1 )*A( K+1, I ) ) )
566                     IF( TEMP1.GT.HALF ) THEN
567                        A( K, I ) = HALF*SQRT( ABS( A( K+1,
568     $                              I-1 )*A( K+1, I ) ) )
569                     END IF
570   90             CONTINUE
571*
572               ELSE
573*
574                  IINFO = 1
575               END IF
576*
577               IF( IINFO.NE.0 ) THEN
578                  WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N,
579     $               JTYPE, IOLDSD
580                  INFO = ABS( IINFO )
581                  RETURN
582               END IF
583*
584  100          CONTINUE
585*
586*              Call DSBTRD to compute S and U from upper triangle.
587*
588               CALL DLACPY( ' ', K+1, N, A, LDA, WORK, LDA )
589*
590               NTEST = 1
591               CALL DSBTRD( 'V', 'U', N, K, WORK, LDA, SD, SE, U, LDU,
592     $                      WORK( LDA*N+1 ), IINFO )
593*
594               IF( IINFO.NE.0 ) THEN
595                  WRITE( NOUNIT, FMT = 9999 )'DSBTRD(U)', IINFO, N,
596     $               JTYPE, IOLDSD
597                  INFO = ABS( IINFO )
598                  IF( IINFO.LT.0 ) THEN
599                     RETURN
600                  ELSE
601                     RESULT( 1 ) = ULPINV
602                     GO TO 150
603                  END IF
604               END IF
605*
606*              Do tests 1 and 2
607*
608               CALL DSBT21( 'Upper', N, K, 1, A, LDA, SD, SE, U, LDU,
609     $                      WORK, RESULT( 1 ) )
610*
611*              Convert A from Upper-Triangle-Only storage to
612*              Lower-Triangle-Only storage.
613*
614               DO 120 JC = 1, N
615                  DO 110 JR = 0, MIN( K, N-JC )
616                     A( JR+1, JC ) = A( K+1-JR, JC+JR )
617  110             CONTINUE
618  120          CONTINUE
619               DO 140 JC = N + 1 - K, N
620                  DO 130 JR = MIN( K, N-JC ) + 1, K
621                     A( JR+1, JC ) = ZERO
622  130             CONTINUE
623  140          CONTINUE
624*
625*              Call DSBTRD to compute S and U from lower triangle
626*
627               CALL DLACPY( ' ', K+1, N, A, LDA, WORK, LDA )
628*
629               NTEST = 3
630               CALL DSBTRD( 'V', 'L', N, K, WORK, LDA, SD, SE, U, LDU,
631     $                      WORK( LDA*N+1 ), IINFO )
632*
633               IF( IINFO.NE.0 ) THEN
634                  WRITE( NOUNIT, FMT = 9999 )'DSBTRD(L)', IINFO, N,
635     $               JTYPE, IOLDSD
636                  INFO = ABS( IINFO )
637                  IF( IINFO.LT.0 ) THEN
638                     RETURN
639                  ELSE
640                     RESULT( 3 ) = ULPINV
641                     GO TO 150
642                  END IF
643               END IF
644               NTEST = 4
645*
646*              Do tests 3 and 4
647*
648               CALL DSBT21( 'Lower', N, K, 1, A, LDA, SD, SE, U, LDU,
649     $                      WORK, RESULT( 3 ) )
650*
651*              End of Loop -- Check for RESULT(j) > THRESH
652*
653  150          CONTINUE
654               NTESTT = NTESTT + NTEST
655*
656*              Print out tests which fail.
657*
658               DO 160 JR = 1, NTEST
659                  IF( RESULT( JR ).GE.THRESH ) THEN
660*
661*                    If this is the first test to fail,
662*                    print a header to the data file.
663*
664                     IF( NERRS.EQ.0 ) THEN
665                        WRITE( NOUNIT, FMT = 9998 )'DSB'
666                        WRITE( NOUNIT, FMT = 9997 )
667                        WRITE( NOUNIT, FMT = 9996 )
668                        WRITE( NOUNIT, FMT = 9995 )'Symmetric'
669                        WRITE( NOUNIT, FMT = 9994 )'orthogonal', '''',
670     $                     'transpose', ( '''', J = 1, 4 )
671                     END IF
672                     NERRS = NERRS + 1
673                     WRITE( NOUNIT, FMT = 9993 )N, K, IOLDSD, JTYPE,
674     $                  JR, RESULT( JR )
675                  END IF
676  160          CONTINUE
677*
678  170       CONTINUE
679  180    CONTINUE
680  190 CONTINUE
681*
682*     Summary
683*
684      CALL DLASUM( 'DSB', NOUNIT, NERRS, NTESTT )
685      RETURN
686*
687 9999 FORMAT( ' DCHKSB: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
688     $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
689*
690 9998 FORMAT( / 1X, A3,
691     $      ' -- Real Symmetric Banded Tridiagonal Reduction Routines' )
692 9997 FORMAT( ' Matrix types (see DCHKSB for details): ' )
693*
694 9996 FORMAT( / ' Special Matrices:',
695     $      / '  1=Zero matrix.                        ',
696     $      '  5=Diagonal: clustered entries.',
697     $      / '  2=Identity matrix.                    ',
698     $      '  6=Diagonal: large, evenly spaced.',
699     $      / '  3=Diagonal: evenly spaced entries.    ',
700     $      '  7=Diagonal: small, evenly spaced.',
701     $      / '  4=Diagonal: geometr. spaced entries.' )
702 9995 FORMAT( ' Dense ', A, ' Banded Matrices:',
703     $      / '  8=Evenly spaced eigenvals.            ',
704     $      ' 12=Small, evenly spaced eigenvals.',
705     $      / '  9=Geometrically spaced eigenvals.     ',
706     $      ' 13=Matrix with random O(1) entries.',
707     $      / ' 10=Clustered eigenvalues.              ',
708     $      ' 14=Matrix with large random entries.',
709     $      / ' 11=Large, evenly spaced eigenvals.     ',
710     $      ' 15=Matrix with small random entries.' )
711*
712 9994 FORMAT( / ' Tests performed:   (S is Tridiag,  U is ', A, ',',
713     $      / 20X, A, ' means ', A, '.', / ' UPLO=''U'':',
714     $      / '  1= | A - U S U', A1, ' | / ( |A| n ulp )     ',
715     $      '  2= | I - U U', A1, ' | / ( n ulp )', / ' UPLO=''L'':',
716     $      / '  3= | A - U S U', A1, ' | / ( |A| n ulp )     ',
717     $      '  4= | I - U U', A1, ' | / ( n ulp )' )
718 9993 FORMAT( ' N=', I5, ', K=', I4, ', seed=', 4( I4, ',' ), ' type ',
719     $      I2, ', test(', I2, ')=', G10.3 )
720*
721*     End of DCHKSB
722*
723      END
724