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