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