LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dchksb.f
Go to the documentation of this file.
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 *> \date November 2011
288 *
289 *> \ingroup double_eig
290 *
291 * =====================================================================
292  SUBROUTINE dchksb( NSIZES, NN, NWDTHS, KK, NTYPES, DOTYPE, ISEED,
293  $ thresh, nounit, a, lda, sd, se, u, ldu, work,
294  $ lwork, result, info )
295 *
296 * -- LAPACK test routine (version 3.4.0) --
297 * -- LAPACK is a software package provided by Univ. of Tennessee, --
298 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
299 * November 2011
300 *
301 * .. Scalar Arguments ..
302  INTEGER info, lda, ldu, lwork, nounit, nsizes, ntypes,
303  $ nwdths
304  DOUBLE PRECISION thresh
305 * ..
306 * .. Array Arguments ..
307  LOGICAL dotype( * )
308  INTEGER iseed( 4 ), kk( * ), nn( * )
309  DOUBLE PRECISION a( lda, * ), result( * ), sd( * ), se( * ),
310  $ u( ldu, * ), work( * )
311 * ..
312 *
313 * =====================================================================
314 *
315 * .. Parameters ..
316  DOUBLE PRECISION zero, one, two, ten
317  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
318  $ ten = 10.0d0 )
319  DOUBLE PRECISION half
320  parameter( half = one / two )
321  INTEGER maxtyp
322  parameter( maxtyp = 15 )
323 * ..
324 * .. Local Scalars ..
325  LOGICAL badnn, badnnb
326  INTEGER i, iinfo, imode, itype, j, jc, jcol, jr, jsize,
327  $ jtype, jwidth, k, kmax, mtypes, n, nerrs,
328  $ nmats, nmax, ntest, ntestt
329  DOUBLE PRECISION aninv, anorm, cond, ovfl, rtovfl, rtunfl,
330  $ temp1, ulp, ulpinv, unfl
331 * ..
332 * .. Local Arrays ..
333  INTEGER idumma( 1 ), ioldsd( 4 ), kmagn( maxtyp ),
334  $ kmode( maxtyp ), ktype( maxtyp )
335 * ..
336 * .. External Functions ..
337  DOUBLE PRECISION dlamch
338  EXTERNAL dlamch
339 * ..
340 * .. External Subroutines ..
341  EXTERNAL dlacpy, dlaset, dlasum, dlatmr, dlatms, dsbt21,
342  $ dsbtrd, xerbla
343 * ..
344 * .. Intrinsic Functions ..
345  INTRINSIC abs, dble, max, min, sqrt
346 * ..
347 * .. Data statements ..
348  DATA ktype / 1, 2, 5*4, 5*5, 3*8 /
349  DATA kmagn / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
350  $ 2, 3 /
351  DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
352  $ 0, 0 /
353 * ..
354 * .. Executable Statements ..
355 *
356 * Check for errors
357 *
358  ntestt = 0
359  info = 0
360 *
361 * Important constants
362 *
363  badnn = .false.
364  nmax = 1
365  DO 10 j = 1, nsizes
366  nmax = max( nmax, nn( j ) )
367  IF( nn( j ).LT.0 )
368  $ badnn = .true.
369  10 continue
370 *
371  badnnb = .false.
372  kmax = 0
373  DO 20 j = 1, nsizes
374  kmax = max( kmax, kk( j ) )
375  IF( kk( j ).LT.0 )
376  $ badnnb = .true.
377  20 continue
378  kmax = min( nmax-1, kmax )
379 *
380 * Check for errors
381 *
382  IF( nsizes.LT.0 ) THEN
383  info = -1
384  ELSE IF( badnn ) THEN
385  info = -2
386  ELSE IF( nwdths.LT.0 ) THEN
387  info = -3
388  ELSE IF( badnnb ) THEN
389  info = -4
390  ELSE IF( ntypes.LT.0 ) THEN
391  info = -5
392  ELSE IF( lda.LT.kmax+1 ) THEN
393  info = -11
394  ELSE IF( ldu.LT.nmax ) THEN
395  info = -15
396  ELSE IF( ( max( lda, nmax )+1 )*nmax.GT.lwork ) THEN
397  info = -17
398  END IF
399 *
400  IF( info.NE.0 ) THEN
401  CALL xerbla( 'DCHKSB', -info )
402  return
403  END IF
404 *
405 * Quick return if possible
406 *
407  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 .OR. nwdths.EQ.0 )
408  $ return
409 *
410 * More Important constants
411 *
412  unfl = dlamch( 'Safe minimum' )
413  ovfl = one / unfl
414  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
415  ulpinv = one / ulp
416  rtunfl = sqrt( unfl )
417  rtovfl = sqrt( ovfl )
418 *
419 * Loop over sizes, types
420 *
421  nerrs = 0
422  nmats = 0
423 *
424  DO 190 jsize = 1, nsizes
425  n = nn( jsize )
426  aninv = one / dble( max( 1, n ) )
427 *
428  DO 180 jwidth = 1, nwdths
429  k = kk( jwidth )
430  IF( k.GT.n )
431  $ go to 180
432  k = max( 0, min( n-1, k ) )
433 *
434  IF( nsizes.NE.1 ) THEN
435  mtypes = min( maxtyp, ntypes )
436  ELSE
437  mtypes = min( maxtyp+1, ntypes )
438  END IF
439 *
440  DO 170 jtype = 1, mtypes
441  IF( .NOT.dotype( jtype ) )
442  $ go to 170
443  nmats = nmats + 1
444  ntest = 0
445 *
446  DO 30 j = 1, 4
447  ioldsd( j ) = iseed( j )
448  30 continue
449 *
450 * Compute "A".
451 * Store as "Upper"; later, we will copy to other format.
452 *
453 * Control parameters:
454 *
455 * KMAGN KMODE KTYPE
456 * =1 O(1) clustered 1 zero
457 * =2 large clustered 2 identity
458 * =3 small exponential (none)
459 * =4 arithmetic diagonal, (w/ eigenvalues)
460 * =5 random log symmetric, w/ eigenvalues
461 * =6 random (none)
462 * =7 random diagonal
463 * =8 random symmetric
464 * =9 positive definite
465 * =10 diagonally dominant tridiagonal
466 *
467  IF( mtypes.GT.maxtyp )
468  $ go to 100
469 *
470  itype = ktype( jtype )
471  imode = kmode( jtype )
472 *
473 * Compute norm
474 *
475  go to( 40, 50, 60 )kmagn( jtype )
476 *
477  40 continue
478  anorm = one
479  go to 70
480 *
481  50 continue
482  anorm = ( rtovfl*ulp )*aninv
483  go to 70
484 *
485  60 continue
486  anorm = rtunfl*n*ulpinv
487  go to 70
488 *
489  70 continue
490 *
491  CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
492  iinfo = 0
493  IF( jtype.LE.15 ) THEN
494  cond = ulpinv
495  ELSE
496  cond = ulpinv*aninv / ten
497  END IF
498 *
499 * Special Matrices -- Identity & Jordan block
500 *
501 * Zero
502 *
503  IF( itype.EQ.1 ) THEN
504  iinfo = 0
505 *
506  ELSE IF( itype.EQ.2 ) THEN
507 *
508 * Identity
509 *
510  DO 80 jcol = 1, n
511  a( k+1, jcol ) = anorm
512  80 continue
513 *
514  ELSE IF( itype.EQ.4 ) THEN
515 *
516 * Diagonal Matrix, [Eigen]values Specified
517 *
518  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
519  $ anorm, 0, 0, 'Q', a( k+1, 1 ), lda,
520  $ work( n+1 ), iinfo )
521 *
522  ELSE IF( itype.EQ.5 ) THEN
523 *
524 * Symmetric, eigenvalues specified
525 *
526  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
527  $ anorm, k, k, 'Q', a, lda, work( n+1 ),
528  $ iinfo )
529 *
530  ELSE IF( itype.EQ.7 ) THEN
531 *
532 * Diagonal, random eigenvalues
533 *
534  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
535  $ 'T', 'N', work( n+1 ), 1, one,
536  $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
537  $ zero, anorm, 'Q', a( k+1, 1 ), lda,
538  $ idumma, iinfo )
539 *
540  ELSE IF( itype.EQ.8 ) THEN
541 *
542 * Symmetric, random eigenvalues
543 *
544  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
545  $ 'T', 'N', work( n+1 ), 1, one,
546  $ work( 2*n+1 ), 1, one, 'N', idumma, k, k,
547  $ zero, anorm, 'Q', a, lda, idumma, iinfo )
548 *
549  ELSE IF( itype.EQ.9 ) THEN
550 *
551 * Positive definite, eigenvalues specified.
552 *
553  CALL dlatms( n, n, 'S', iseed, 'P', work, imode, cond,
554  $ anorm, k, k, 'Q', a, lda, work( n+1 ),
555  $ iinfo )
556 *
557  ELSE IF( itype.EQ.10 ) THEN
558 *
559 * Positive definite tridiagonal, eigenvalues specified.
560 *
561  IF( n.GT.1 )
562  $ k = max( 1, k )
563  CALL dlatms( n, n, 'S', iseed, 'P', work, imode, cond,
564  $ anorm, 1, 1, 'Q', a( k, 1 ), lda,
565  $ work( n+1 ), iinfo )
566  DO 90 i = 2, n
567  temp1 = abs( a( k, i ) ) /
568  $ sqrt( abs( a( k+1, i-1 )*a( k+1, i ) ) )
569  IF( temp1.GT.half ) THEN
570  a( k, i ) = half*sqrt( abs( a( k+1,
571  $ i-1 )*a( k+1, i ) ) )
572  END IF
573  90 continue
574 *
575  ELSE
576 *
577  iinfo = 1
578  END IF
579 *
580  IF( iinfo.NE.0 ) THEN
581  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n,
582  $ jtype, ioldsd
583  info = abs( iinfo )
584  return
585  END IF
586 *
587  100 continue
588 *
589 * Call DSBTRD to compute S and U from upper triangle.
590 *
591  CALL dlacpy( ' ', k+1, n, a, lda, work, lda )
592 *
593  ntest = 1
594  CALL dsbtrd( 'V', 'U', n, k, work, lda, sd, se, u, ldu,
595  $ work( lda*n+1 ), iinfo )
596 *
597  IF( iinfo.NE.0 ) THEN
598  WRITE( nounit, fmt = 9999 )'DSBTRD(U)', iinfo, n,
599  $ jtype, ioldsd
600  info = abs( iinfo )
601  IF( iinfo.LT.0 ) THEN
602  return
603  ELSE
604  result( 1 ) = ulpinv
605  go to 150
606  END IF
607  END IF
608 *
609 * Do tests 1 and 2
610 *
611  CALL dsbt21( 'Upper', n, k, 1, a, lda, sd, se, u, ldu,
612  $ work, result( 1 ) )
613 *
614 * Convert A from Upper-Triangle-Only storage to
615 * Lower-Triangle-Only storage.
616 *
617  DO 120 jc = 1, n
618  DO 110 jr = 0, min( k, n-jc )
619  a( jr+1, jc ) = a( k+1-jr, jc+jr )
620  110 continue
621  120 continue
622  DO 140 jc = n + 1 - k, n
623  DO 130 jr = min( k, n-jc ) + 1, k
624  a( jr+1, jc ) = zero
625  130 continue
626  140 continue
627 *
628 * Call DSBTRD to compute S and U from lower triangle
629 *
630  CALL dlacpy( ' ', k+1, n, a, lda, work, lda )
631 *
632  ntest = 3
633  CALL dsbtrd( 'V', 'L', n, k, work, lda, sd, se, u, ldu,
634  $ work( lda*n+1 ), iinfo )
635 *
636  IF( iinfo.NE.0 ) THEN
637  WRITE( nounit, fmt = 9999 )'DSBTRD(L)', iinfo, n,
638  $ jtype, ioldsd
639  info = abs( iinfo )
640  IF( iinfo.LT.0 ) THEN
641  return
642  ELSE
643  result( 3 ) = ulpinv
644  go to 150
645  END IF
646  END IF
647  ntest = 4
648 *
649 * Do tests 3 and 4
650 *
651  CALL dsbt21( 'Lower', n, k, 1, a, lda, sd, se, u, ldu,
652  $ work, result( 3 ) )
653 *
654 * End of Loop -- Check for RESULT(j) > THRESH
655 *
656  150 continue
657  ntestt = ntestt + ntest
658 *
659 * Print out tests which fail.
660 *
661  DO 160 jr = 1, ntest
662  IF( result( jr ).GE.thresh ) THEN
663 *
664 * If this is the first test to fail,
665 * print a header to the data file.
666 *
667  IF( nerrs.EQ.0 ) THEN
668  WRITE( nounit, fmt = 9998 )'DSB'
669  WRITE( nounit, fmt = 9997 )
670  WRITE( nounit, fmt = 9996 )
671  WRITE( nounit, fmt = 9995 )'Symmetric'
672  WRITE( nounit, fmt = 9994 )'orthogonal', '''',
673  $ 'transpose', ( '''', j = 1, 4 )
674  END IF
675  nerrs = nerrs + 1
676  WRITE( nounit, fmt = 9993 )n, k, ioldsd, jtype,
677  $ jr, result( jr )
678  END IF
679  160 continue
680 *
681  170 continue
682  180 continue
683  190 continue
684 *
685 * Summary
686 *
687  CALL dlasum( 'DSB', nounit, nerrs, ntestt )
688  return
689 *
690  9999 format( ' DCHKSB: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
691  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
692 *
693  9998 format( / 1x, a3,
694  $ ' -- Real Symmetric Banded Tridiagonal Reduction Routines' )
695  9997 format( ' Matrix types (see DCHKSB for details): ' )
696 *
697  9996 format( / ' Special Matrices:',
698  $ / ' 1=Zero matrix. ',
699  $ ' 5=Diagonal: clustered entries.',
700  $ / ' 2=Identity matrix. ',
701  $ ' 6=Diagonal: large, evenly spaced.',
702  $ / ' 3=Diagonal: evenly spaced entries. ',
703  $ ' 7=Diagonal: small, evenly spaced.',
704  $ / ' 4=Diagonal: geometr. spaced entries.' )
705  9995 format( ' Dense ', a, ' Banded Matrices:',
706  $ / ' 8=Evenly spaced eigenvals. ',
707  $ ' 12=Small, evenly spaced eigenvals.',
708  $ / ' 9=Geometrically spaced eigenvals. ',
709  $ ' 13=Matrix with random O(1) entries.',
710  $ / ' 10=Clustered eigenvalues. ',
711  $ ' 14=Matrix with large random entries.',
712  $ / ' 11=Large, evenly spaced eigenvals. ',
713  $ ' 15=Matrix with small random entries.' )
714 *
715  9994 format( / ' Tests performed: (S is Tridiag, U is ', a, ',',
716  $ / 20x, a, ' means ', a, '.', / ' UPLO=''U'':',
717  $ / ' 1= | A - U S U', a1, ' | / ( |A| n ulp ) ',
718  $ ' 2= | I - U U', a1, ' | / ( n ulp )', / ' UPLO=''L'':',
719  $ / ' 3= | A - U S U', a1, ' | / ( |A| n ulp ) ',
720  $ ' 4= | I - U U', a1, ' | / ( n ulp )' )
721  9993 format( ' N=', i5, ', K=', i4, ', seed=', 4( i4, ',' ), ' type ',
722  $ i2, ', test(', i2, ')=', g10.3 )
723 *
724 * End of DCHKSB
725 *
726  END