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