1 SUBROUTINE CCHKBD( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
2 $ ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX,
3 $ Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK,
4 $ RWORK, NOUT, INFO )
5 *
6 * -- LAPACK test routine (version 3.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
12 $ NSIZES, NTYPES
13 REAL THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL DOTYPE( * )
17 INTEGER ISEED( 4 ), MVAL( * ), NVAL( * )
18 REAL BD( * ), BE( * ), RWORK( * ), S1( * ), S2( * )
19 COMPLEX A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
20 $ U( LDPT, * ), VT( LDPT, * ), WORK( * ),
21 $ X( LDX, * ), Y( LDX, * ), Z( LDX, * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * CCHKBD checks the singular value decomposition (SVD) routines.
28 *
29 * CGEBRD reduces a complex general m by n matrix A to real upper or
30 * lower bidiagonal form by an orthogonal transformation: Q' * A * P = B
31 * (or A = Q * B * P'). The matrix B is upper bidiagonal if m >= n
32 * and lower bidiagonal if m < n.
33 *
34 * CUNGBR generates the orthogonal matrices Q and P' from CGEBRD.
35 * Note that Q and P are not necessarily square.
36 *
37 * CBDSQR computes the singular value decomposition of the bidiagonal
38 * matrix B as B = U S V'. It is called three times to compute
39 * 1) B = U S1 V', where S1 is the diagonal matrix of singular
40 * values and the columns of the matrices U and V are the left
41 * and right singular vectors, respectively, of B.
42 * 2) Same as 1), but the singular values are stored in S2 and the
43 * singular vectors are not computed.
44 * 3) A = (UQ) S (P'V'), the SVD of the original matrix A.
45 * In addition, CBDSQR has an option to apply the left orthogonal matrix
46 * U to a matrix X, useful in least squares applications.
47 *
48 * For each pair of matrix dimensions (M,N) and each selected matrix
49 * type, an M by N matrix A and an M by NRHS matrix X are generated.
50 * The problem dimensions are as follows
51 * A: M x N
52 * Q: M x min(M,N) (but M x M if NRHS > 0)
53 * P: min(M,N) x N
54 * B: min(M,N) x min(M,N)
55 * U, V: min(M,N) x min(M,N)
56 * S1, S2 diagonal, order min(M,N)
57 * X: M x NRHS
58 *
59 * For each generated matrix, 14 tests are performed:
60 *
61 * Test CGEBRD and CUNGBR
62 *
63 * (1) | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P'
64 *
65 * (2) | I - Q' Q | / ( M ulp )
66 *
67 * (3) | I - PT PT' | / ( N ulp )
68 *
69 * Test CBDSQR on bidiagonal matrix B
70 *
71 * (4) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'
72 *
73 * (5) | Y - U Z | / ( |Y| max(min(M,N),k) ulp ), where Y = Q' X
74 * and Z = U' Y.
75 * (6) | I - U' U | / ( min(M,N) ulp )
76 *
77 * (7) | I - VT VT' | / ( min(M,N) ulp )
78 *
79 * (8) S1 contains min(M,N) nonnegative values in decreasing order.
80 * (Return 0 if true, 1/ULP if false.)
81 *
82 * (9) 0 if the true singular values of B are within THRESH of
83 * those in S1. 2*THRESH if they are not. (Tested using
84 * SSVDCH)
85 *
86 * (10) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
87 * computing U and V.
88 *
89 * Test CBDSQR on matrix A
90 *
91 * (11) | A - (QU) S (VT PT) | / ( |A| max(M,N) ulp )
92 *
93 * (12) | X - (QU) Z | / ( |X| max(M,k) ulp )
94 *
95 * (13) | I - (QU)'(QU) | / ( M ulp )
96 *
97 * (14) | I - (VT PT) (PT'VT') | / ( N ulp )
98 *
99 * The possible matrix types are
100 *
101 * (1) The zero matrix.
102 * (2) The identity matrix.
103 *
104 * (3) A diagonal matrix with evenly spaced entries
105 * 1, ..., ULP and random signs.
106 * (ULP = (first number larger than 1) - 1 )
107 * (4) A diagonal matrix with geometrically spaced entries
108 * 1, ..., ULP and random signs.
109 * (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
110 * and random signs.
111 *
112 * (6) Same as (3), but multiplied by SQRT( overflow threshold )
113 * (7) Same as (3), but multiplied by SQRT( underflow threshold )
114 *
115 * (8) A matrix of the form U D V, where U and V are orthogonal and
116 * D has evenly spaced entries 1, ..., ULP with random signs
117 * on the diagonal.
118 *
119 * (9) A matrix of the form U D V, where U and V are orthogonal and
120 * D has geometrically spaced entries 1, ..., ULP with random
121 * signs on the diagonal.
122 *
123 * (10) A matrix of the form U D V, where U and V are orthogonal and
124 * D has "clustered" entries 1, ULP,..., ULP with random
125 * signs on the diagonal.
126 *
127 * (11) Same as (8), but multiplied by SQRT( overflow threshold )
128 * (12) Same as (8), but multiplied by SQRT( underflow threshold )
129 *
130 * (13) Rectangular matrix with random entries chosen from (-1,1).
131 * (14) Same as (13), but multiplied by SQRT( overflow threshold )
132 * (15) Same as (13), but multiplied by SQRT( underflow threshold )
133 *
134 * Special case:
135 * (16) A bidiagonal matrix with random entries chosen from a
136 * logarithmic distribution on [ulp^2,ulp^(-2)] (I.e., each
137 * entry is e^x, where x is chosen uniformly on
138 * [ 2 log(ulp), -2 log(ulp) ] .) For *this* type:
139 * (a) CGEBRD is not called to reduce it to bidiagonal form.
140 * (b) the bidiagonal is min(M,N) x min(M,N); if M<N, the
141 * matrix will be lower bidiagonal, otherwise upper.
142 * (c) only tests 5--8 and 14 are performed.
143 *
144 * A subset of the full set of matrix types may be selected through
145 * the logical array DOTYPE.
146 *
147 * Arguments
148 * ==========
149 *
150 * NSIZES (input) INTEGER
151 * The number of values of M and N contained in the vectors
152 * MVAL and NVAL. The matrix sizes are used in pairs (M,N).
153 *
154 * MVAL (input) INTEGER array, dimension (NM)
155 * The values of the matrix row dimension M.
156 *
157 * NVAL (input) INTEGER array, dimension (NM)
158 * The values of the matrix column dimension N.
159 *
160 * NTYPES (input) INTEGER
161 * The number of elements in DOTYPE. If it is zero, CCHKBD
162 * does nothing. It must be at least zero. If it is MAXTYP+1
163 * and NSIZES is 1, then an additional type, MAXTYP+1 is
164 * defined, which is to use whatever matrices are in A and B.
165 * This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
166 * DOTYPE(MAXTYP+1) is .TRUE. .
167 *
168 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
169 * If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
170 * of type j will be generated. If NTYPES is smaller than the
171 * maximum number of types defined (PARAMETER MAXTYP), then
172 * types NTYPES+1 through MAXTYP will not be generated. If
173 * NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
174 * DOTYPE(NTYPES) will be ignored.
175 *
176 * NRHS (input) INTEGER
177 * The number of columns in the "right-hand side" matrices X, Y,
178 * and Z, used in testing CBDSQR. If NRHS = 0, then the
179 * operations on the right-hand side will not be tested.
180 * NRHS must be at least 0.
181 *
182 * ISEED (input/output) INTEGER array, dimension (4)
183 * On entry ISEED specifies the seed of the random number
184 * generator. The array elements should be between 0 and 4095;
185 * if not they will be reduced mod 4096. Also, ISEED(4) must
186 * be odd. The values of ISEED are changed on exit, and can be
187 * used in the next call to CCHKBD to continue the same random
188 * number sequence.
189 *
190 * THRESH (input) REAL
191 * The threshold value for the test ratios. A result is
192 * included in the output file if RESULT >= THRESH. To have
193 * every test ratio printed, use THRESH = 0. Note that the
194 * expected value of the test ratios is O(1), so THRESH should
195 * be a reasonably small multiple of 1, e.g., 10 or 100.
196 *
197 * A (workspace) COMPLEX array, dimension (LDA,NMAX)
198 * where NMAX is the maximum value of N in NVAL.
199 *
200 * LDA (input) INTEGER
201 * The leading dimension of the array A. LDA >= max(1,MMAX),
202 * where MMAX is the maximum value of M in MVAL.
203 *
204 * BD (workspace) REAL array, dimension
205 * (max(min(MVAL(j),NVAL(j))))
206 *
207 * BE (workspace) REAL array, dimension
208 * (max(min(MVAL(j),NVAL(j))))
209 *
210 * S1 (workspace) REAL array, dimension
211 * (max(min(MVAL(j),NVAL(j))))
212 *
213 * S2 (workspace) REAL array, dimension
214 * (max(min(MVAL(j),NVAL(j))))
215 *
216 * X (workspace) COMPLEX array, dimension (LDX,NRHS)
217 *
218 * LDX (input) INTEGER
219 * The leading dimension of the arrays X, Y, and Z.
220 * LDX >= max(1,MMAX).
221 *
222 * Y (workspace) COMPLEX array, dimension (LDX,NRHS)
223 *
224 * Z (workspace) COMPLEX array, dimension (LDX,NRHS)
225 *
226 * Q (workspace) COMPLEX array, dimension (LDQ,MMAX)
227 *
228 * LDQ (input) INTEGER
229 * The leading dimension of the array Q. LDQ >= max(1,MMAX).
230 *
231 * PT (workspace) COMPLEX array, dimension (LDPT,NMAX)
232 *
233 * LDPT (input) INTEGER
234 * The leading dimension of the arrays PT, U, and V.
235 * LDPT >= max(1, max(min(MVAL(j),NVAL(j)))).
236 *
237 * U (workspace) COMPLEX array, dimension
238 * (LDPT,max(min(MVAL(j),NVAL(j))))
239 *
240 * V (workspace) COMPLEX array, dimension
241 * (LDPT,max(min(MVAL(j),NVAL(j))))
242 *
243 * WORK (workspace) COMPLEX array, dimension (LWORK)
244 *
245 * LWORK (input) INTEGER
246 * The number of entries in WORK. This must be at least
247 * 3(M+N) and M(M + max(M,N,k) + 1) + N*min(M,N) for all
248 * pairs (M,N)=(MM(j),NN(j))
249 *
250 * RWORK (workspace) REAL array, dimension
251 * (5*max(min(M,N)))
252 *
253 * NOUT (input) INTEGER
254 * The FORTRAN unit number for printing out error messages
255 * (e.g., if a routine returns IINFO not equal to 0.)
256 *
257 * INFO (output) INTEGER
258 * If 0, then everything ran OK.
259 * -1: NSIZES < 0
260 * -2: Some MM(j) < 0
261 * -3: Some NN(j) < 0
262 * -4: NTYPES < 0
263 * -6: NRHS < 0
264 * -8: THRESH < 0
265 * -11: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
266 * -17: LDB < 1 or LDB < MMAX.
267 * -21: LDQ < 1 or LDQ < MMAX.
268 * -23: LDP < 1 or LDP < MNMAX.
269 * -27: LWORK too small.
270 * If CLATMR, CLATMS, CGEBRD, CUNGBR, or CBDSQR,
271 * returns an error code, the
272 * absolute value of it is returned.
273 *
274 *-----------------------------------------------------------------------
275 *
276 * Some Local Variables and Parameters:
277 * ---- ----- --------- --- ----------
278 *
279 * ZERO, ONE Real 0 and 1.
280 * MAXTYP The number of types defined.
281 * NTEST The number of tests performed, or which can
282 * be performed so far, for the current matrix.
283 * MMAX Largest value in NN.
284 * NMAX Largest value in NN.
285 * MNMIN min(MM(j), NN(j)) (the dimension of the bidiagonal
286 * matrix.)
287 * MNMAX The maximum value of MNMIN for j=1,...,NSIZES.
288 * NFAIL The number of tests which have exceeded THRESH
289 * COND, IMODE Values to be passed to the matrix generators.
290 * ANORM Norm of A; passed to matrix generators.
291 *
292 * OVFL, UNFL Overflow and underflow thresholds.
293 * RTOVFL, RTUNFL Square roots of the previous 2 values.
294 * ULP, ULPINV Finest relative precision and its inverse.
295 *
296 * The following four arrays decode JTYPE:
297 * KTYPE(j) The general type (1-10) for type "j".
298 * KMODE(j) The MODE value to be passed to the matrix
299 * generator for type "j".
300 * KMAGN(j) The order of magnitude ( O(1),
301 * O(overflow^(1/2) ), O(underflow^(1/2) )
302 *
303 * ======================================================================
304 *
305 * .. Parameters ..
306 REAL ZERO, ONE, TWO, HALF
307 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
308 $ HALF = 0.5E0 )
309 COMPLEX CZERO, CONE
310 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
311 $ CONE = ( 1.0E+0, 0.0E+0 ) )
312 INTEGER MAXTYP
313 PARAMETER ( MAXTYP = 16 )
314 * ..
315 * .. Local Scalars ..
316 LOGICAL BADMM, BADNN, BIDIAG
317 CHARACTER UPLO
318 CHARACTER*3 PATH
319 INTEGER I, IINFO, IMODE, ITYPE, J, JCOL, JSIZE, JTYPE,
320 $ LOG2UI, M, MINWRK, MMAX, MNMAX, MNMIN, MQ,
321 $ MTYPES, N, NFAIL, NMAX, NTEST
322 REAL AMNINV, ANORM, COND, OVFL, RTOVFL, RTUNFL,
323 $ TEMP1, TEMP2, ULP, ULPINV, UNFL
324 * ..
325 * .. Local Arrays ..
326 INTEGER IOLDSD( 4 ), IWORK( 1 ), KMAGN( MAXTYP ),
327 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
328 REAL DUMMA( 1 ), RESULT( 14 )
329 * ..
330 * .. External Functions ..
331 REAL SLAMCH, SLARND
332 EXTERNAL SLAMCH, SLARND
333 * ..
334 * .. External Subroutines ..
335 EXTERNAL ALASUM, CBDSQR, CBDT01, CBDT02, CBDT03, CGEBRD,
336 $ CGEMM, CLACPY, CLASET, CLATMR, CLATMS, CUNGBR,
337 $ CUNT01, SCOPY, SLABAD, SLAHD2, SSVDCH, XERBLA
338 * ..
339 * .. Intrinsic Functions ..
340 INTRINSIC ABS, EXP, INT, LOG, MAX, MIN, SQRT
341 * ..
342 * .. Scalars in Common ..
343 LOGICAL LERR, OK
344 CHARACTER*32 SRNAMT
345 INTEGER INFOT, NUNIT
346 * ..
347 * .. Common blocks ..
348 COMMON / INFOC / INFOT, NUNIT, OK, LERR
349 COMMON / SRNAMC / SRNAMT
350 * ..
351 * .. Data statements ..
352 DATA KTYPE / 1, 2, 5*4, 5*6, 3*9, 10 /
353 DATA KMAGN / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
354 DATA KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
355 $ 0, 0, 0 /
356 * ..
357 * .. Executable Statements ..
358 *
359 * Check for errors
360 *
361 INFO = 0
362 *
363 BADMM = .FALSE.
364 BADNN = .FALSE.
365 MMAX = 1
366 NMAX = 1
367 MNMAX = 1
368 MINWRK = 1
369 DO 10 J = 1, NSIZES
370 MMAX = MAX( MMAX, MVAL( J ) )
371 IF( MVAL( J ).LT.0 )
372 $ BADMM = .TRUE.
373 NMAX = MAX( NMAX, NVAL( J ) )
374 IF( NVAL( J ).LT.0 )
375 $ BADNN = .TRUE.
376 MNMAX = MAX( MNMAX, MIN( MVAL( J ), NVAL( J ) ) )
377 MINWRK = MAX( MINWRK, 3*( MVAL( J )+NVAL( J ) ),
378 $ MVAL( J )*( MVAL( J )+MAX( MVAL( J ), NVAL( J ),
379 $ NRHS )+1 )+NVAL( J )*MIN( NVAL( J ), MVAL( J ) ) )
380 10 CONTINUE
381 *
382 * Check for errors
383 *
384 IF( NSIZES.LT.0 ) THEN
385 INFO = -1
386 ELSE IF( BADMM ) THEN
387 INFO = -2
388 ELSE IF( BADNN ) THEN
389 INFO = -3
390 ELSE IF( NTYPES.LT.0 ) THEN
391 INFO = -4
392 ELSE IF( NRHS.LT.0 ) THEN
393 INFO = -6
394 ELSE IF( LDA.LT.MMAX ) THEN
395 INFO = -11
396 ELSE IF( LDX.LT.MMAX ) THEN
397 INFO = -17
398 ELSE IF( LDQ.LT.MMAX ) THEN
399 INFO = -21
400 ELSE IF( LDPT.LT.MNMAX ) THEN
401 INFO = -23
402 ELSE IF( MINWRK.GT.LWORK ) THEN
403 INFO = -27
404 END IF
405 *
406 IF( INFO.NE.0 ) THEN
407 CALL XERBLA( 'CCHKBD', -INFO )
408 RETURN
409 END IF
410 *
411 * Initialize constants
412 *
413 PATH( 1: 1 ) = 'Complex precision'
414 PATH( 2: 3 ) = 'BD'
415 NFAIL = 0
416 NTEST = 0
417 UNFL = SLAMCH( 'Safe minimum' )
418 OVFL = SLAMCH( 'Overflow' )
419 CALL SLABAD( UNFL, OVFL )
420 ULP = SLAMCH( 'Precision' )
421 ULPINV = ONE / ULP
422 LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) )
423 RTUNFL = SQRT( UNFL )
424 RTOVFL = SQRT( OVFL )
425 INFOT = 0
426 *
427 * Loop over sizes, types
428 *
429 DO 180 JSIZE = 1, NSIZES
430 M = MVAL( JSIZE )
431 N = NVAL( JSIZE )
432 MNMIN = MIN( M, N )
433 AMNINV = ONE / MAX( M, N, 1 )
434 *
435 IF( NSIZES.NE.1 ) THEN
436 MTYPES = MIN( MAXTYP, NTYPES )
437 ELSE
438 MTYPES = MIN( MAXTYP+1, NTYPES )
439 END IF
440 *
441 DO 170 JTYPE = 1, MTYPES
442 IF( .NOT.DOTYPE( JTYPE ) )
443 $ GO TO 170
444 *
445 DO 20 J = 1, 4
446 IOLDSD( J ) = ISEED( J )
447 20 CONTINUE
448 *
449 DO 30 J = 1, 14
450 RESULT( J ) = -ONE
451 30 CONTINUE
452 *
453 UPLO = ' '
454 *
455 * Compute "A"
456 *
457 * Control parameters:
458 *
459 * KMAGN KMODE KTYPE
460 * =1 O(1) clustered 1 zero
461 * =2 large clustered 2 identity
462 * =3 small exponential (none)
463 * =4 arithmetic diagonal, (w/ eigenvalues)
464 * =5 random symmetric, w/ eigenvalues
465 * =6 nonsymmetric, w/ singular values
466 * =7 random diagonal
467 * =8 random symmetric
468 * =9 random nonsymmetric
469 * =10 random bidiagonal (log. distrib.)
470 *
471 IF( MTYPES.GT.MAXTYP )
472 $ GO TO 100
473 *
474 ITYPE = KTYPE( JTYPE )
475 IMODE = KMODE( JTYPE )
476 *
477 * Compute norm
478 *
479 GO TO ( 40, 50, 60 )KMAGN( JTYPE )
480 *
481 40 CONTINUE
482 ANORM = ONE
483 GO TO 70
484 *
485 50 CONTINUE
486 ANORM = ( RTOVFL*ULP )*AMNINV
487 GO TO 70
488 *
489 60 CONTINUE
490 ANORM = RTUNFL*MAX( M, N )*ULPINV
491 GO TO 70
492 *
493 70 CONTINUE
494 *
495 CALL CLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
496 IINFO = 0
497 COND = ULPINV
498 *
499 BIDIAG = .FALSE.
500 IF( ITYPE.EQ.1 ) THEN
501 *
502 * Zero matrix
503 *
504 IINFO = 0
505 *
506 ELSE IF( ITYPE.EQ.2 ) THEN
507 *
508 * Identity
509 *
510 DO 80 JCOL = 1, MNMIN
511 A( JCOL, JCOL ) = ANORM
512 80 CONTINUE
513 *
514 ELSE IF( ITYPE.EQ.4 ) THEN
515 *
516 * Diagonal Matrix, [Eigen]values Specified
517 *
518 CALL CLATMS( MNMIN, MNMIN, 'S', ISEED, 'N', RWORK, IMODE,
519 $ COND, ANORM, 0, 0, 'N', A, LDA, WORK,
520 $ IINFO )
521 *
522 ELSE IF( ITYPE.EQ.5 ) THEN
523 *
524 * Symmetric, eigenvalues specified
525 *
526 CALL CLATMS( MNMIN, MNMIN, 'S', ISEED, 'S', RWORK, IMODE,
527 $ COND, ANORM, M, N, 'N', A, LDA, WORK,
528 $ IINFO )
529 *
530 ELSE IF( ITYPE.EQ.6 ) THEN
531 *
532 * Nonsymmetric, singular values specified
533 *
534 CALL CLATMS( M, N, 'S', ISEED, 'N', RWORK, IMODE, COND,
535 $ ANORM, M, N, 'N', A, LDA, WORK, IINFO )
536 *
537 ELSE IF( ITYPE.EQ.7 ) THEN
538 *
539 * Diagonal, random entries
540 *
541 CALL CLATMR( MNMIN, MNMIN, 'S', ISEED, 'N', WORK, 6, ONE,
542 $ CONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE,
543 $ WORK( 2*MNMIN+1 ), 1, ONE, 'N', IWORK, 0, 0,
544 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
545 *
546 ELSE IF( ITYPE.EQ.8 ) THEN
547 *
548 * Symmetric, random entries
549 *
550 CALL CLATMR( MNMIN, MNMIN, 'S', ISEED, 'S', WORK, 6, ONE,
551 $ CONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE,
552 $ WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N,
553 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
554 *
555 ELSE IF( ITYPE.EQ.9 ) THEN
556 *
557 * Nonsymmetric, random entries
558 *
559 CALL CLATMR( M, N, 'S', ISEED, 'N', WORK, 6, ONE, CONE,
560 $ 'T', 'N', WORK( MNMIN+1 ), 1, ONE,
561 $ WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N,
562 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
563 *
564 ELSE IF( ITYPE.EQ.10 ) THEN
565 *
566 * Bidiagonal, random entries
567 *
568 TEMP1 = -TWO*LOG( ULP )
569 DO 90 J = 1, MNMIN
570 BD( J ) = EXP( TEMP1*SLARND( 2, ISEED ) )
571 IF( J.LT.MNMIN )
572 $ BE( J ) = EXP( TEMP1*SLARND( 2, ISEED ) )
573 90 CONTINUE
574 *
575 IINFO = 0
576 BIDIAG = .TRUE.
577 IF( M.GE.N ) THEN
578 UPLO = 'U'
579 ELSE
580 UPLO = 'L'
581 END IF
582 ELSE
583 IINFO = 1
584 END IF
585 *
586 IF( IINFO.EQ.0 ) THEN
587 *
588 * Generate Right-Hand Side
589 *
590 IF( BIDIAG ) THEN
591 CALL CLATMR( MNMIN, NRHS, 'S', ISEED, 'N', WORK, 6,
592 $ ONE, CONE, 'T', 'N', WORK( MNMIN+1 ), 1,
593 $ ONE, WORK( 2*MNMIN+1 ), 1, ONE, 'N',
594 $ IWORK, MNMIN, NRHS, ZERO, ONE, 'NO', Y,
595 $ LDX, IWORK, IINFO )
596 ELSE
597 CALL CLATMR( M, NRHS, 'S', ISEED, 'N', WORK, 6, ONE,
598 $ CONE, 'T', 'N', WORK( M+1 ), 1, ONE,
599 $ WORK( 2*M+1 ), 1, ONE, 'N', IWORK, M,
600 $ NRHS, ZERO, ONE, 'NO', X, LDX, IWORK,
601 $ IINFO )
602 END IF
603 END IF
604 *
605 * Error Exit
606 *
607 IF( IINFO.NE.0 ) THEN
608 WRITE( NOUT, FMT = 9998 )'Generator', IINFO, M, N,
609 $ JTYPE, IOLDSD
610 INFO = ABS( IINFO )
611 RETURN
612 END IF
613 *
614 100 CONTINUE
615 *
616 * Call CGEBRD and CUNGBR to compute B, Q, and P, do tests.
617 *
618 IF( .NOT.BIDIAG ) THEN
619 *
620 * Compute transformations to reduce A to bidiagonal form:
621 * B := Q' * A * P.
622 *
623 CALL CLACPY( ' ', M, N, A, LDA, Q, LDQ )
624 CALL CGEBRD( M, N, Q, LDQ, BD, BE, WORK, WORK( MNMIN+1 ),
625 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
626 *
627 * Check error code from CGEBRD.
628 *
629 IF( IINFO.NE.0 ) THEN
630 WRITE( NOUT, FMT = 9998 )'CGEBRD', IINFO, M, N,
631 $ JTYPE, IOLDSD
632 INFO = ABS( IINFO )
633 RETURN
634 END IF
635 *
636 CALL CLACPY( ' ', M, N, Q, LDQ, PT, LDPT )
637 IF( M.GE.N ) THEN
638 UPLO = 'U'
639 ELSE
640 UPLO = 'L'
641 END IF
642 *
643 * Generate Q
644 *
645 MQ = M
646 IF( NRHS.LE.0 )
647 $ MQ = MNMIN
648 CALL CUNGBR( 'Q', M, MQ, N, Q, LDQ, WORK,
649 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
650 *
651 * Check error code from CUNGBR.
652 *
653 IF( IINFO.NE.0 ) THEN
654 WRITE( NOUT, FMT = 9998 )'CUNGBR(Q)', IINFO, M, N,
655 $ JTYPE, IOLDSD
656 INFO = ABS( IINFO )
657 RETURN
658 END IF
659 *
660 * Generate P'
661 *
662 CALL CUNGBR( 'P', MNMIN, N, M, PT, LDPT, WORK( MNMIN+1 ),
663 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
664 *
665 * Check error code from CUNGBR.
666 *
667 IF( IINFO.NE.0 ) THEN
668 WRITE( NOUT, FMT = 9998 )'CUNGBR(P)', IINFO, M, N,
669 $ JTYPE, IOLDSD
670 INFO = ABS( IINFO )
671 RETURN
672 END IF
673 *
674 * Apply Q' to an M by NRHS matrix X: Y := Q' * X.
675 *
676 CALL CGEMM( 'Conjugate transpose', 'No transpose', M,
677 $ NRHS, M, CONE, Q, LDQ, X, LDX, CZERO, Y,
678 $ LDX )
679 *
680 * Test 1: Check the decomposition A := Q * B * PT
681 * 2: Check the orthogonality of Q
682 * 3: Check the orthogonality of PT
683 *
684 CALL CBDT01( M, N, 1, A, LDA, Q, LDQ, BD, BE, PT, LDPT,
685 $ WORK, RWORK, RESULT( 1 ) )
686 CALL CUNT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
687 $ RWORK, RESULT( 2 ) )
688 CALL CUNT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
689 $ RWORK, RESULT( 3 ) )
690 END IF
691 *
692 * Use CBDSQR to form the SVD of the bidiagonal matrix B:
693 * B := U * S1 * VT, and compute Z = U' * Y.
694 *
695 CALL SCOPY( MNMIN, BD, 1, S1, 1 )
696 IF( MNMIN.GT.0 )
697 $ CALL SCOPY( MNMIN-1, BE, 1, RWORK, 1 )
698 CALL CLACPY( ' ', M, NRHS, Y, LDX, Z, LDX )
699 CALL CLASET( 'Full', MNMIN, MNMIN, CZERO, CONE, U, LDPT )
700 CALL CLASET( 'Full', MNMIN, MNMIN, CZERO, CONE, VT, LDPT )
701 *
702 CALL CBDSQR( UPLO, MNMIN, MNMIN, MNMIN, NRHS, S1, RWORK, VT,
703 $ LDPT, U, LDPT, Z, LDX, RWORK( MNMIN+1 ),
704 $ IINFO )
705 *
706 * Check error code from CBDSQR.
707 *
708 IF( IINFO.NE.0 ) THEN
709 WRITE( NOUT, FMT = 9998 )'CBDSQR(vects)', IINFO, M, N,
710 $ JTYPE, IOLDSD
711 INFO = ABS( IINFO )
712 IF( IINFO.LT.0 ) THEN
713 RETURN
714 ELSE
715 RESULT( 4 ) = ULPINV
716 GO TO 150
717 END IF
718 END IF
719 *
720 * Use CBDSQR to compute only the singular values of the
721 * bidiagonal matrix B; U, VT, and Z should not be modified.
722 *
723 CALL SCOPY( MNMIN, BD, 1, S2, 1 )
724 IF( MNMIN.GT.0 )
725 $ CALL SCOPY( MNMIN-1, BE, 1, RWORK, 1 )
726 *
727 CALL CBDSQR( UPLO, MNMIN, 0, 0, 0, S2, RWORK, VT, LDPT, U,
728 $ LDPT, Z, LDX, RWORK( MNMIN+1 ), IINFO )
729 *
730 * Check error code from CBDSQR.
731 *
732 IF( IINFO.NE.0 ) THEN
733 WRITE( NOUT, FMT = 9998 )'CBDSQR(values)', IINFO, M, N,
734 $ JTYPE, IOLDSD
735 INFO = ABS( IINFO )
736 IF( IINFO.LT.0 ) THEN
737 RETURN
738 ELSE
739 RESULT( 9 ) = ULPINV
740 GO TO 150
741 END IF
742 END IF
743 *
744 * Test 4: Check the decomposition B := U * S1 * VT
745 * 5: Check the computation Z := U' * Y
746 * 6: Check the orthogonality of U
747 * 7: Check the orthogonality of VT
748 *
749 CALL CBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT,
750 $ WORK, RESULT( 4 ) )
751 CALL CBDT02( MNMIN, NRHS, Y, LDX, Z, LDX, U, LDPT, WORK,
752 $ RWORK, RESULT( 5 ) )
753 CALL CUNT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK,
754 $ RWORK, RESULT( 6 ) )
755 CALL CUNT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK,
756 $ RWORK, RESULT( 7 ) )
757 *
758 * Test 8: Check that the singular values are sorted in
759 * non-increasing order and are non-negative
760 *
761 RESULT( 8 ) = ZERO
762 DO 110 I = 1, MNMIN - 1
763 IF( S1( I ).LT.S1( I+1 ) )
764 $ RESULT( 8 ) = ULPINV
765 IF( S1( I ).LT.ZERO )
766 $ RESULT( 8 ) = ULPINV
767 110 CONTINUE
768 IF( MNMIN.GE.1 ) THEN
769 IF( S1( MNMIN ).LT.ZERO )
770 $ RESULT( 8 ) = ULPINV
771 END IF
772 *
773 * Test 9: Compare CBDSQR with and without singular vectors
774 *
775 TEMP2 = ZERO
776 *
777 DO 120 J = 1, MNMIN
778 TEMP1 = ABS( S1( J )-S2( J ) ) /
779 $ MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ),
780 $ ULP*MAX( ABS( S1( J ) ), ABS( S2( J ) ) ) )
781 TEMP2 = MAX( TEMP1, TEMP2 )
782 120 CONTINUE
783 *
784 RESULT( 9 ) = TEMP2
785 *
786 * Test 10: Sturm sequence test of singular values
787 * Go up by factors of two until it succeeds
788 *
789 TEMP1 = THRESH*( HALF-ULP )
790 *
791 DO 130 J = 0, LOG2UI
792 CALL SSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO )
793 IF( IINFO.EQ.0 )
794 $ GO TO 140
795 TEMP1 = TEMP1*TWO
796 130 CONTINUE
797 *
798 140 CONTINUE
799 RESULT( 10 ) = TEMP1
800 *
801 * Use CBDSQR to form the decomposition A := (QU) S (VT PT)
802 * from the bidiagonal form A := Q B PT.
803 *
804 IF( .NOT.BIDIAG ) THEN
805 CALL SCOPY( MNMIN, BD, 1, S2, 1 )
806 IF( MNMIN.GT.0 )
807 $ CALL SCOPY( MNMIN-1, BE, 1, RWORK, 1 )
808 *
809 CALL CBDSQR( UPLO, MNMIN, N, M, NRHS, S2, RWORK, PT,
810 $ LDPT, Q, LDQ, Y, LDX, RWORK( MNMIN+1 ),
811 $ IINFO )
812 *
813 * Test 11: Check the decomposition A := Q*U * S2 * VT*PT
814 * 12: Check the computation Z := U' * Q' * X
815 * 13: Check the orthogonality of Q*U
816 * 14: Check the orthogonality of VT*PT
817 *
818 CALL CBDT01( M, N, 0, A, LDA, Q, LDQ, S2, DUMMA, PT,
819 $ LDPT, WORK, RWORK, RESULT( 11 ) )
820 CALL CBDT02( M, NRHS, X, LDX, Y, LDX, Q, LDQ, WORK,
821 $ RWORK, RESULT( 12 ) )
822 CALL CUNT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
823 $ RWORK, RESULT( 13 ) )
824 CALL CUNT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
825 $ RWORK, RESULT( 14 ) )
826 END IF
827 *
828 * End of Loop -- Check for RESULT(j) > THRESH
829 *
830 150 CONTINUE
831 DO 160 J = 1, 14
832 IF( RESULT( J ).GE.THRESH ) THEN
833 IF( NFAIL.EQ.0 )
834 $ CALL SLAHD2( NOUT, PATH )
835 WRITE( NOUT, FMT = 9999 )M, N, JTYPE, IOLDSD, J,
836 $ RESULT( J )
837 NFAIL = NFAIL + 1
838 END IF
839 160 CONTINUE
840 IF( .NOT.BIDIAG ) THEN
841 NTEST = NTEST + 14
842 ELSE
843 NTEST = NTEST + 5
844 END IF
845 *
846 170 CONTINUE
847 180 CONTINUE
848 *
849 * Summary
850 *
851 CALL ALASUM( PATH, NOUT, NFAIL, NTEST, 0 )
852 *
853 RETURN
854 *
855 * End of CCHKBD
856 *
857 9999 FORMAT( ' M=', I5, ', N=', I5, ', type ', I2, ', seed=',
858 $ 4( I4, ',' ), ' test(', I2, ')=', G11.4 )
859 9998 FORMAT( ' CCHKBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
860 $ I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
861 $ I5, ')' )
862 *
863 END
2 $ ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX,
3 $ Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK,
4 $ RWORK, NOUT, INFO )
5 *
6 * -- LAPACK test routine (version 3.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
12 $ NSIZES, NTYPES
13 REAL THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL DOTYPE( * )
17 INTEGER ISEED( 4 ), MVAL( * ), NVAL( * )
18 REAL BD( * ), BE( * ), RWORK( * ), S1( * ), S2( * )
19 COMPLEX A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
20 $ U( LDPT, * ), VT( LDPT, * ), WORK( * ),
21 $ X( LDX, * ), Y( LDX, * ), Z( LDX, * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * CCHKBD checks the singular value decomposition (SVD) routines.
28 *
29 * CGEBRD reduces a complex general m by n matrix A to real upper or
30 * lower bidiagonal form by an orthogonal transformation: Q' * A * P = B
31 * (or A = Q * B * P'). The matrix B is upper bidiagonal if m >= n
32 * and lower bidiagonal if m < n.
33 *
34 * CUNGBR generates the orthogonal matrices Q and P' from CGEBRD.
35 * Note that Q and P are not necessarily square.
36 *
37 * CBDSQR computes the singular value decomposition of the bidiagonal
38 * matrix B as B = U S V'. It is called three times to compute
39 * 1) B = U S1 V', where S1 is the diagonal matrix of singular
40 * values and the columns of the matrices U and V are the left
41 * and right singular vectors, respectively, of B.
42 * 2) Same as 1), but the singular values are stored in S2 and the
43 * singular vectors are not computed.
44 * 3) A = (UQ) S (P'V'), the SVD of the original matrix A.
45 * In addition, CBDSQR has an option to apply the left orthogonal matrix
46 * U to a matrix X, useful in least squares applications.
47 *
48 * For each pair of matrix dimensions (M,N) and each selected matrix
49 * type, an M by N matrix A and an M by NRHS matrix X are generated.
50 * The problem dimensions are as follows
51 * A: M x N
52 * Q: M x min(M,N) (but M x M if NRHS > 0)
53 * P: min(M,N) x N
54 * B: min(M,N) x min(M,N)
55 * U, V: min(M,N) x min(M,N)
56 * S1, S2 diagonal, order min(M,N)
57 * X: M x NRHS
58 *
59 * For each generated matrix, 14 tests are performed:
60 *
61 * Test CGEBRD and CUNGBR
62 *
63 * (1) | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P'
64 *
65 * (2) | I - Q' Q | / ( M ulp )
66 *
67 * (3) | I - PT PT' | / ( N ulp )
68 *
69 * Test CBDSQR on bidiagonal matrix B
70 *
71 * (4) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'
72 *
73 * (5) | Y - U Z | / ( |Y| max(min(M,N),k) ulp ), where Y = Q' X
74 * and Z = U' Y.
75 * (6) | I - U' U | / ( min(M,N) ulp )
76 *
77 * (7) | I - VT VT' | / ( min(M,N) ulp )
78 *
79 * (8) S1 contains min(M,N) nonnegative values in decreasing order.
80 * (Return 0 if true, 1/ULP if false.)
81 *
82 * (9) 0 if the true singular values of B are within THRESH of
83 * those in S1. 2*THRESH if they are not. (Tested using
84 * SSVDCH)
85 *
86 * (10) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
87 * computing U and V.
88 *
89 * Test CBDSQR on matrix A
90 *
91 * (11) | A - (QU) S (VT PT) | / ( |A| max(M,N) ulp )
92 *
93 * (12) | X - (QU) Z | / ( |X| max(M,k) ulp )
94 *
95 * (13) | I - (QU)'(QU) | / ( M ulp )
96 *
97 * (14) | I - (VT PT) (PT'VT') | / ( N ulp )
98 *
99 * The possible matrix types are
100 *
101 * (1) The zero matrix.
102 * (2) The identity matrix.
103 *
104 * (3) A diagonal matrix with evenly spaced entries
105 * 1, ..., ULP and random signs.
106 * (ULP = (first number larger than 1) - 1 )
107 * (4) A diagonal matrix with geometrically spaced entries
108 * 1, ..., ULP and random signs.
109 * (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
110 * and random signs.
111 *
112 * (6) Same as (3), but multiplied by SQRT( overflow threshold )
113 * (7) Same as (3), but multiplied by SQRT( underflow threshold )
114 *
115 * (8) A matrix of the form U D V, where U and V are orthogonal and
116 * D has evenly spaced entries 1, ..., ULP with random signs
117 * on the diagonal.
118 *
119 * (9) A matrix of the form U D V, where U and V are orthogonal and
120 * D has geometrically spaced entries 1, ..., ULP with random
121 * signs on the diagonal.
122 *
123 * (10) A matrix of the form U D V, where U and V are orthogonal and
124 * D has "clustered" entries 1, ULP,..., ULP with random
125 * signs on the diagonal.
126 *
127 * (11) Same as (8), but multiplied by SQRT( overflow threshold )
128 * (12) Same as (8), but multiplied by SQRT( underflow threshold )
129 *
130 * (13) Rectangular matrix with random entries chosen from (-1,1).
131 * (14) Same as (13), but multiplied by SQRT( overflow threshold )
132 * (15) Same as (13), but multiplied by SQRT( underflow threshold )
133 *
134 * Special case:
135 * (16) A bidiagonal matrix with random entries chosen from a
136 * logarithmic distribution on [ulp^2,ulp^(-2)] (I.e., each
137 * entry is e^x, where x is chosen uniformly on
138 * [ 2 log(ulp), -2 log(ulp) ] .) For *this* type:
139 * (a) CGEBRD is not called to reduce it to bidiagonal form.
140 * (b) the bidiagonal is min(M,N) x min(M,N); if M<N, the
141 * matrix will be lower bidiagonal, otherwise upper.
142 * (c) only tests 5--8 and 14 are performed.
143 *
144 * A subset of the full set of matrix types may be selected through
145 * the logical array DOTYPE.
146 *
147 * Arguments
148 * ==========
149 *
150 * NSIZES (input) INTEGER
151 * The number of values of M and N contained in the vectors
152 * MVAL and NVAL. The matrix sizes are used in pairs (M,N).
153 *
154 * MVAL (input) INTEGER array, dimension (NM)
155 * The values of the matrix row dimension M.
156 *
157 * NVAL (input) INTEGER array, dimension (NM)
158 * The values of the matrix column dimension N.
159 *
160 * NTYPES (input) INTEGER
161 * The number of elements in DOTYPE. If it is zero, CCHKBD
162 * does nothing. It must be at least zero. If it is MAXTYP+1
163 * and NSIZES is 1, then an additional type, MAXTYP+1 is
164 * defined, which is to use whatever matrices are in A and B.
165 * This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
166 * DOTYPE(MAXTYP+1) is .TRUE. .
167 *
168 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
169 * If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
170 * of type j will be generated. If NTYPES is smaller than the
171 * maximum number of types defined (PARAMETER MAXTYP), then
172 * types NTYPES+1 through MAXTYP will not be generated. If
173 * NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
174 * DOTYPE(NTYPES) will be ignored.
175 *
176 * NRHS (input) INTEGER
177 * The number of columns in the "right-hand side" matrices X, Y,
178 * and Z, used in testing CBDSQR. If NRHS = 0, then the
179 * operations on the right-hand side will not be tested.
180 * NRHS must be at least 0.
181 *
182 * ISEED (input/output) INTEGER array, dimension (4)
183 * On entry ISEED specifies the seed of the random number
184 * generator. The array elements should be between 0 and 4095;
185 * if not they will be reduced mod 4096. Also, ISEED(4) must
186 * be odd. The values of ISEED are changed on exit, and can be
187 * used in the next call to CCHKBD to continue the same random
188 * number sequence.
189 *
190 * THRESH (input) REAL
191 * The threshold value for the test ratios. A result is
192 * included in the output file if RESULT >= THRESH. To have
193 * every test ratio printed, use THRESH = 0. Note that the
194 * expected value of the test ratios is O(1), so THRESH should
195 * be a reasonably small multiple of 1, e.g., 10 or 100.
196 *
197 * A (workspace) COMPLEX array, dimension (LDA,NMAX)
198 * where NMAX is the maximum value of N in NVAL.
199 *
200 * LDA (input) INTEGER
201 * The leading dimension of the array A. LDA >= max(1,MMAX),
202 * where MMAX is the maximum value of M in MVAL.
203 *
204 * BD (workspace) REAL array, dimension
205 * (max(min(MVAL(j),NVAL(j))))
206 *
207 * BE (workspace) REAL array, dimension
208 * (max(min(MVAL(j),NVAL(j))))
209 *
210 * S1 (workspace) REAL array, dimension
211 * (max(min(MVAL(j),NVAL(j))))
212 *
213 * S2 (workspace) REAL array, dimension
214 * (max(min(MVAL(j),NVAL(j))))
215 *
216 * X (workspace) COMPLEX array, dimension (LDX,NRHS)
217 *
218 * LDX (input) INTEGER
219 * The leading dimension of the arrays X, Y, and Z.
220 * LDX >= max(1,MMAX).
221 *
222 * Y (workspace) COMPLEX array, dimension (LDX,NRHS)
223 *
224 * Z (workspace) COMPLEX array, dimension (LDX,NRHS)
225 *
226 * Q (workspace) COMPLEX array, dimension (LDQ,MMAX)
227 *
228 * LDQ (input) INTEGER
229 * The leading dimension of the array Q. LDQ >= max(1,MMAX).
230 *
231 * PT (workspace) COMPLEX array, dimension (LDPT,NMAX)
232 *
233 * LDPT (input) INTEGER
234 * The leading dimension of the arrays PT, U, and V.
235 * LDPT >= max(1, max(min(MVAL(j),NVAL(j)))).
236 *
237 * U (workspace) COMPLEX array, dimension
238 * (LDPT,max(min(MVAL(j),NVAL(j))))
239 *
240 * V (workspace) COMPLEX array, dimension
241 * (LDPT,max(min(MVAL(j),NVAL(j))))
242 *
243 * WORK (workspace) COMPLEX array, dimension (LWORK)
244 *
245 * LWORK (input) INTEGER
246 * The number of entries in WORK. This must be at least
247 * 3(M+N) and M(M + max(M,N,k) + 1) + N*min(M,N) for all
248 * pairs (M,N)=(MM(j),NN(j))
249 *
250 * RWORK (workspace) REAL array, dimension
251 * (5*max(min(M,N)))
252 *
253 * NOUT (input) INTEGER
254 * The FORTRAN unit number for printing out error messages
255 * (e.g., if a routine returns IINFO not equal to 0.)
256 *
257 * INFO (output) INTEGER
258 * If 0, then everything ran OK.
259 * -1: NSIZES < 0
260 * -2: Some MM(j) < 0
261 * -3: Some NN(j) < 0
262 * -4: NTYPES < 0
263 * -6: NRHS < 0
264 * -8: THRESH < 0
265 * -11: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
266 * -17: LDB < 1 or LDB < MMAX.
267 * -21: LDQ < 1 or LDQ < MMAX.
268 * -23: LDP < 1 or LDP < MNMAX.
269 * -27: LWORK too small.
270 * If CLATMR, CLATMS, CGEBRD, CUNGBR, or CBDSQR,
271 * returns an error code, the
272 * absolute value of it is returned.
273 *
274 *-----------------------------------------------------------------------
275 *
276 * Some Local Variables and Parameters:
277 * ---- ----- --------- --- ----------
278 *
279 * ZERO, ONE Real 0 and 1.
280 * MAXTYP The number of types defined.
281 * NTEST The number of tests performed, or which can
282 * be performed so far, for the current matrix.
283 * MMAX Largest value in NN.
284 * NMAX Largest value in NN.
285 * MNMIN min(MM(j), NN(j)) (the dimension of the bidiagonal
286 * matrix.)
287 * MNMAX The maximum value of MNMIN for j=1,...,NSIZES.
288 * NFAIL The number of tests which have exceeded THRESH
289 * COND, IMODE Values to be passed to the matrix generators.
290 * ANORM Norm of A; passed to matrix generators.
291 *
292 * OVFL, UNFL Overflow and underflow thresholds.
293 * RTOVFL, RTUNFL Square roots of the previous 2 values.
294 * ULP, ULPINV Finest relative precision and its inverse.
295 *
296 * The following four arrays decode JTYPE:
297 * KTYPE(j) The general type (1-10) for type "j".
298 * KMODE(j) The MODE value to be passed to the matrix
299 * generator for type "j".
300 * KMAGN(j) The order of magnitude ( O(1),
301 * O(overflow^(1/2) ), O(underflow^(1/2) )
302 *
303 * ======================================================================
304 *
305 * .. Parameters ..
306 REAL ZERO, ONE, TWO, HALF
307 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
308 $ HALF = 0.5E0 )
309 COMPLEX CZERO, CONE
310 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
311 $ CONE = ( 1.0E+0, 0.0E+0 ) )
312 INTEGER MAXTYP
313 PARAMETER ( MAXTYP = 16 )
314 * ..
315 * .. Local Scalars ..
316 LOGICAL BADMM, BADNN, BIDIAG
317 CHARACTER UPLO
318 CHARACTER*3 PATH
319 INTEGER I, IINFO, IMODE, ITYPE, J, JCOL, JSIZE, JTYPE,
320 $ LOG2UI, M, MINWRK, MMAX, MNMAX, MNMIN, MQ,
321 $ MTYPES, N, NFAIL, NMAX, NTEST
322 REAL AMNINV, ANORM, COND, OVFL, RTOVFL, RTUNFL,
323 $ TEMP1, TEMP2, ULP, ULPINV, UNFL
324 * ..
325 * .. Local Arrays ..
326 INTEGER IOLDSD( 4 ), IWORK( 1 ), KMAGN( MAXTYP ),
327 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
328 REAL DUMMA( 1 ), RESULT( 14 )
329 * ..
330 * .. External Functions ..
331 REAL SLAMCH, SLARND
332 EXTERNAL SLAMCH, SLARND
333 * ..
334 * .. External Subroutines ..
335 EXTERNAL ALASUM, CBDSQR, CBDT01, CBDT02, CBDT03, CGEBRD,
336 $ CGEMM, CLACPY, CLASET, CLATMR, CLATMS, CUNGBR,
337 $ CUNT01, SCOPY, SLABAD, SLAHD2, SSVDCH, XERBLA
338 * ..
339 * .. Intrinsic Functions ..
340 INTRINSIC ABS, EXP, INT, LOG, MAX, MIN, SQRT
341 * ..
342 * .. Scalars in Common ..
343 LOGICAL LERR, OK
344 CHARACTER*32 SRNAMT
345 INTEGER INFOT, NUNIT
346 * ..
347 * .. Common blocks ..
348 COMMON / INFOC / INFOT, NUNIT, OK, LERR
349 COMMON / SRNAMC / SRNAMT
350 * ..
351 * .. Data statements ..
352 DATA KTYPE / 1, 2, 5*4, 5*6, 3*9, 10 /
353 DATA KMAGN / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
354 DATA KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
355 $ 0, 0, 0 /
356 * ..
357 * .. Executable Statements ..
358 *
359 * Check for errors
360 *
361 INFO = 0
362 *
363 BADMM = .FALSE.
364 BADNN = .FALSE.
365 MMAX = 1
366 NMAX = 1
367 MNMAX = 1
368 MINWRK = 1
369 DO 10 J = 1, NSIZES
370 MMAX = MAX( MMAX, MVAL( J ) )
371 IF( MVAL( J ).LT.0 )
372 $ BADMM = .TRUE.
373 NMAX = MAX( NMAX, NVAL( J ) )
374 IF( NVAL( J ).LT.0 )
375 $ BADNN = .TRUE.
376 MNMAX = MAX( MNMAX, MIN( MVAL( J ), NVAL( J ) ) )
377 MINWRK = MAX( MINWRK, 3*( MVAL( J )+NVAL( J ) ),
378 $ MVAL( J )*( MVAL( J )+MAX( MVAL( J ), NVAL( J ),
379 $ NRHS )+1 )+NVAL( J )*MIN( NVAL( J ), MVAL( J ) ) )
380 10 CONTINUE
381 *
382 * Check for errors
383 *
384 IF( NSIZES.LT.0 ) THEN
385 INFO = -1
386 ELSE IF( BADMM ) THEN
387 INFO = -2
388 ELSE IF( BADNN ) THEN
389 INFO = -3
390 ELSE IF( NTYPES.LT.0 ) THEN
391 INFO = -4
392 ELSE IF( NRHS.LT.0 ) THEN
393 INFO = -6
394 ELSE IF( LDA.LT.MMAX ) THEN
395 INFO = -11
396 ELSE IF( LDX.LT.MMAX ) THEN
397 INFO = -17
398 ELSE IF( LDQ.LT.MMAX ) THEN
399 INFO = -21
400 ELSE IF( LDPT.LT.MNMAX ) THEN
401 INFO = -23
402 ELSE IF( MINWRK.GT.LWORK ) THEN
403 INFO = -27
404 END IF
405 *
406 IF( INFO.NE.0 ) THEN
407 CALL XERBLA( 'CCHKBD', -INFO )
408 RETURN
409 END IF
410 *
411 * Initialize constants
412 *
413 PATH( 1: 1 ) = 'Complex precision'
414 PATH( 2: 3 ) = 'BD'
415 NFAIL = 0
416 NTEST = 0
417 UNFL = SLAMCH( 'Safe minimum' )
418 OVFL = SLAMCH( 'Overflow' )
419 CALL SLABAD( UNFL, OVFL )
420 ULP = SLAMCH( 'Precision' )
421 ULPINV = ONE / ULP
422 LOG2UI = INT( LOG( ULPINV ) / LOG( TWO ) )
423 RTUNFL = SQRT( UNFL )
424 RTOVFL = SQRT( OVFL )
425 INFOT = 0
426 *
427 * Loop over sizes, types
428 *
429 DO 180 JSIZE = 1, NSIZES
430 M = MVAL( JSIZE )
431 N = NVAL( JSIZE )
432 MNMIN = MIN( M, N )
433 AMNINV = ONE / MAX( M, N, 1 )
434 *
435 IF( NSIZES.NE.1 ) THEN
436 MTYPES = MIN( MAXTYP, NTYPES )
437 ELSE
438 MTYPES = MIN( MAXTYP+1, NTYPES )
439 END IF
440 *
441 DO 170 JTYPE = 1, MTYPES
442 IF( .NOT.DOTYPE( JTYPE ) )
443 $ GO TO 170
444 *
445 DO 20 J = 1, 4
446 IOLDSD( J ) = ISEED( J )
447 20 CONTINUE
448 *
449 DO 30 J = 1, 14
450 RESULT( J ) = -ONE
451 30 CONTINUE
452 *
453 UPLO = ' '
454 *
455 * Compute "A"
456 *
457 * Control parameters:
458 *
459 * KMAGN KMODE KTYPE
460 * =1 O(1) clustered 1 zero
461 * =2 large clustered 2 identity
462 * =3 small exponential (none)
463 * =4 arithmetic diagonal, (w/ eigenvalues)
464 * =5 random symmetric, w/ eigenvalues
465 * =6 nonsymmetric, w/ singular values
466 * =7 random diagonal
467 * =8 random symmetric
468 * =9 random nonsymmetric
469 * =10 random bidiagonal (log. distrib.)
470 *
471 IF( MTYPES.GT.MAXTYP )
472 $ GO TO 100
473 *
474 ITYPE = KTYPE( JTYPE )
475 IMODE = KMODE( JTYPE )
476 *
477 * Compute norm
478 *
479 GO TO ( 40, 50, 60 )KMAGN( JTYPE )
480 *
481 40 CONTINUE
482 ANORM = ONE
483 GO TO 70
484 *
485 50 CONTINUE
486 ANORM = ( RTOVFL*ULP )*AMNINV
487 GO TO 70
488 *
489 60 CONTINUE
490 ANORM = RTUNFL*MAX( M, N )*ULPINV
491 GO TO 70
492 *
493 70 CONTINUE
494 *
495 CALL CLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
496 IINFO = 0
497 COND = ULPINV
498 *
499 BIDIAG = .FALSE.
500 IF( ITYPE.EQ.1 ) THEN
501 *
502 * Zero matrix
503 *
504 IINFO = 0
505 *
506 ELSE IF( ITYPE.EQ.2 ) THEN
507 *
508 * Identity
509 *
510 DO 80 JCOL = 1, MNMIN
511 A( JCOL, JCOL ) = ANORM
512 80 CONTINUE
513 *
514 ELSE IF( ITYPE.EQ.4 ) THEN
515 *
516 * Diagonal Matrix, [Eigen]values Specified
517 *
518 CALL CLATMS( MNMIN, MNMIN, 'S', ISEED, 'N', RWORK, IMODE,
519 $ COND, ANORM, 0, 0, 'N', A, LDA, WORK,
520 $ IINFO )
521 *
522 ELSE IF( ITYPE.EQ.5 ) THEN
523 *
524 * Symmetric, eigenvalues specified
525 *
526 CALL CLATMS( MNMIN, MNMIN, 'S', ISEED, 'S', RWORK, IMODE,
527 $ COND, ANORM, M, N, 'N', A, LDA, WORK,
528 $ IINFO )
529 *
530 ELSE IF( ITYPE.EQ.6 ) THEN
531 *
532 * Nonsymmetric, singular values specified
533 *
534 CALL CLATMS( M, N, 'S', ISEED, 'N', RWORK, IMODE, COND,
535 $ ANORM, M, N, 'N', A, LDA, WORK, IINFO )
536 *
537 ELSE IF( ITYPE.EQ.7 ) THEN
538 *
539 * Diagonal, random entries
540 *
541 CALL CLATMR( MNMIN, MNMIN, 'S', ISEED, 'N', WORK, 6, ONE,
542 $ CONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE,
543 $ WORK( 2*MNMIN+1 ), 1, ONE, 'N', IWORK, 0, 0,
544 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
545 *
546 ELSE IF( ITYPE.EQ.8 ) THEN
547 *
548 * Symmetric, random entries
549 *
550 CALL CLATMR( MNMIN, MNMIN, 'S', ISEED, 'S', WORK, 6, ONE,
551 $ CONE, 'T', 'N', WORK( MNMIN+1 ), 1, ONE,
552 $ WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N,
553 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
554 *
555 ELSE IF( ITYPE.EQ.9 ) THEN
556 *
557 * Nonsymmetric, random entries
558 *
559 CALL CLATMR( M, N, 'S', ISEED, 'N', WORK, 6, ONE, CONE,
560 $ 'T', 'N', WORK( MNMIN+1 ), 1, ONE,
561 $ WORK( M+MNMIN+1 ), 1, ONE, 'N', IWORK, M, N,
562 $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
563 *
564 ELSE IF( ITYPE.EQ.10 ) THEN
565 *
566 * Bidiagonal, random entries
567 *
568 TEMP1 = -TWO*LOG( ULP )
569 DO 90 J = 1, MNMIN
570 BD( J ) = EXP( TEMP1*SLARND( 2, ISEED ) )
571 IF( J.LT.MNMIN )
572 $ BE( J ) = EXP( TEMP1*SLARND( 2, ISEED ) )
573 90 CONTINUE
574 *
575 IINFO = 0
576 BIDIAG = .TRUE.
577 IF( M.GE.N ) THEN
578 UPLO = 'U'
579 ELSE
580 UPLO = 'L'
581 END IF
582 ELSE
583 IINFO = 1
584 END IF
585 *
586 IF( IINFO.EQ.0 ) THEN
587 *
588 * Generate Right-Hand Side
589 *
590 IF( BIDIAG ) THEN
591 CALL CLATMR( MNMIN, NRHS, 'S', ISEED, 'N', WORK, 6,
592 $ ONE, CONE, 'T', 'N', WORK( MNMIN+1 ), 1,
593 $ ONE, WORK( 2*MNMIN+1 ), 1, ONE, 'N',
594 $ IWORK, MNMIN, NRHS, ZERO, ONE, 'NO', Y,
595 $ LDX, IWORK, IINFO )
596 ELSE
597 CALL CLATMR( M, NRHS, 'S', ISEED, 'N', WORK, 6, ONE,
598 $ CONE, 'T', 'N', WORK( M+1 ), 1, ONE,
599 $ WORK( 2*M+1 ), 1, ONE, 'N', IWORK, M,
600 $ NRHS, ZERO, ONE, 'NO', X, LDX, IWORK,
601 $ IINFO )
602 END IF
603 END IF
604 *
605 * Error Exit
606 *
607 IF( IINFO.NE.0 ) THEN
608 WRITE( NOUT, FMT = 9998 )'Generator', IINFO, M, N,
609 $ JTYPE, IOLDSD
610 INFO = ABS( IINFO )
611 RETURN
612 END IF
613 *
614 100 CONTINUE
615 *
616 * Call CGEBRD and CUNGBR to compute B, Q, and P, do tests.
617 *
618 IF( .NOT.BIDIAG ) THEN
619 *
620 * Compute transformations to reduce A to bidiagonal form:
621 * B := Q' * A * P.
622 *
623 CALL CLACPY( ' ', M, N, A, LDA, Q, LDQ )
624 CALL CGEBRD( M, N, Q, LDQ, BD, BE, WORK, WORK( MNMIN+1 ),
625 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
626 *
627 * Check error code from CGEBRD.
628 *
629 IF( IINFO.NE.0 ) THEN
630 WRITE( NOUT, FMT = 9998 )'CGEBRD', IINFO, M, N,
631 $ JTYPE, IOLDSD
632 INFO = ABS( IINFO )
633 RETURN
634 END IF
635 *
636 CALL CLACPY( ' ', M, N, Q, LDQ, PT, LDPT )
637 IF( M.GE.N ) THEN
638 UPLO = 'U'
639 ELSE
640 UPLO = 'L'
641 END IF
642 *
643 * Generate Q
644 *
645 MQ = M
646 IF( NRHS.LE.0 )
647 $ MQ = MNMIN
648 CALL CUNGBR( 'Q', M, MQ, N, Q, LDQ, WORK,
649 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
650 *
651 * Check error code from CUNGBR.
652 *
653 IF( IINFO.NE.0 ) THEN
654 WRITE( NOUT, FMT = 9998 )'CUNGBR(Q)', IINFO, M, N,
655 $ JTYPE, IOLDSD
656 INFO = ABS( IINFO )
657 RETURN
658 END IF
659 *
660 * Generate P'
661 *
662 CALL CUNGBR( 'P', MNMIN, N, M, PT, LDPT, WORK( MNMIN+1 ),
663 $ WORK( 2*MNMIN+1 ), LWORK-2*MNMIN, IINFO )
664 *
665 * Check error code from CUNGBR.
666 *
667 IF( IINFO.NE.0 ) THEN
668 WRITE( NOUT, FMT = 9998 )'CUNGBR(P)', IINFO, M, N,
669 $ JTYPE, IOLDSD
670 INFO = ABS( IINFO )
671 RETURN
672 END IF
673 *
674 * Apply Q' to an M by NRHS matrix X: Y := Q' * X.
675 *
676 CALL CGEMM( 'Conjugate transpose', 'No transpose', M,
677 $ NRHS, M, CONE, Q, LDQ, X, LDX, CZERO, Y,
678 $ LDX )
679 *
680 * Test 1: Check the decomposition A := Q * B * PT
681 * 2: Check the orthogonality of Q
682 * 3: Check the orthogonality of PT
683 *
684 CALL CBDT01( M, N, 1, A, LDA, Q, LDQ, BD, BE, PT, LDPT,
685 $ WORK, RWORK, RESULT( 1 ) )
686 CALL CUNT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
687 $ RWORK, RESULT( 2 ) )
688 CALL CUNT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
689 $ RWORK, RESULT( 3 ) )
690 END IF
691 *
692 * Use CBDSQR to form the SVD of the bidiagonal matrix B:
693 * B := U * S1 * VT, and compute Z = U' * Y.
694 *
695 CALL SCOPY( MNMIN, BD, 1, S1, 1 )
696 IF( MNMIN.GT.0 )
697 $ CALL SCOPY( MNMIN-1, BE, 1, RWORK, 1 )
698 CALL CLACPY( ' ', M, NRHS, Y, LDX, Z, LDX )
699 CALL CLASET( 'Full', MNMIN, MNMIN, CZERO, CONE, U, LDPT )
700 CALL CLASET( 'Full', MNMIN, MNMIN, CZERO, CONE, VT, LDPT )
701 *
702 CALL CBDSQR( UPLO, MNMIN, MNMIN, MNMIN, NRHS, S1, RWORK, VT,
703 $ LDPT, U, LDPT, Z, LDX, RWORK( MNMIN+1 ),
704 $ IINFO )
705 *
706 * Check error code from CBDSQR.
707 *
708 IF( IINFO.NE.0 ) THEN
709 WRITE( NOUT, FMT = 9998 )'CBDSQR(vects)', IINFO, M, N,
710 $ JTYPE, IOLDSD
711 INFO = ABS( IINFO )
712 IF( IINFO.LT.0 ) THEN
713 RETURN
714 ELSE
715 RESULT( 4 ) = ULPINV
716 GO TO 150
717 END IF
718 END IF
719 *
720 * Use CBDSQR to compute only the singular values of the
721 * bidiagonal matrix B; U, VT, and Z should not be modified.
722 *
723 CALL SCOPY( MNMIN, BD, 1, S2, 1 )
724 IF( MNMIN.GT.0 )
725 $ CALL SCOPY( MNMIN-1, BE, 1, RWORK, 1 )
726 *
727 CALL CBDSQR( UPLO, MNMIN, 0, 0, 0, S2, RWORK, VT, LDPT, U,
728 $ LDPT, Z, LDX, RWORK( MNMIN+1 ), IINFO )
729 *
730 * Check error code from CBDSQR.
731 *
732 IF( IINFO.NE.0 ) THEN
733 WRITE( NOUT, FMT = 9998 )'CBDSQR(values)', IINFO, M, N,
734 $ JTYPE, IOLDSD
735 INFO = ABS( IINFO )
736 IF( IINFO.LT.0 ) THEN
737 RETURN
738 ELSE
739 RESULT( 9 ) = ULPINV
740 GO TO 150
741 END IF
742 END IF
743 *
744 * Test 4: Check the decomposition B := U * S1 * VT
745 * 5: Check the computation Z := U' * Y
746 * 6: Check the orthogonality of U
747 * 7: Check the orthogonality of VT
748 *
749 CALL CBDT03( UPLO, MNMIN, 1, BD, BE, U, LDPT, S1, VT, LDPT,
750 $ WORK, RESULT( 4 ) )
751 CALL CBDT02( MNMIN, NRHS, Y, LDX, Z, LDX, U, LDPT, WORK,
752 $ RWORK, RESULT( 5 ) )
753 CALL CUNT01( 'Columns', MNMIN, MNMIN, U, LDPT, WORK, LWORK,
754 $ RWORK, RESULT( 6 ) )
755 CALL CUNT01( 'Rows', MNMIN, MNMIN, VT, LDPT, WORK, LWORK,
756 $ RWORK, RESULT( 7 ) )
757 *
758 * Test 8: Check that the singular values are sorted in
759 * non-increasing order and are non-negative
760 *
761 RESULT( 8 ) = ZERO
762 DO 110 I = 1, MNMIN - 1
763 IF( S1( I ).LT.S1( I+1 ) )
764 $ RESULT( 8 ) = ULPINV
765 IF( S1( I ).LT.ZERO )
766 $ RESULT( 8 ) = ULPINV
767 110 CONTINUE
768 IF( MNMIN.GE.1 ) THEN
769 IF( S1( MNMIN ).LT.ZERO )
770 $ RESULT( 8 ) = ULPINV
771 END IF
772 *
773 * Test 9: Compare CBDSQR with and without singular vectors
774 *
775 TEMP2 = ZERO
776 *
777 DO 120 J = 1, MNMIN
778 TEMP1 = ABS( S1( J )-S2( J ) ) /
779 $ MAX( SQRT( UNFL )*MAX( S1( 1 ), ONE ),
780 $ ULP*MAX( ABS( S1( J ) ), ABS( S2( J ) ) ) )
781 TEMP2 = MAX( TEMP1, TEMP2 )
782 120 CONTINUE
783 *
784 RESULT( 9 ) = TEMP2
785 *
786 * Test 10: Sturm sequence test of singular values
787 * Go up by factors of two until it succeeds
788 *
789 TEMP1 = THRESH*( HALF-ULP )
790 *
791 DO 130 J = 0, LOG2UI
792 CALL SSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO )
793 IF( IINFO.EQ.0 )
794 $ GO TO 140
795 TEMP1 = TEMP1*TWO
796 130 CONTINUE
797 *
798 140 CONTINUE
799 RESULT( 10 ) = TEMP1
800 *
801 * Use CBDSQR to form the decomposition A := (QU) S (VT PT)
802 * from the bidiagonal form A := Q B PT.
803 *
804 IF( .NOT.BIDIAG ) THEN
805 CALL SCOPY( MNMIN, BD, 1, S2, 1 )
806 IF( MNMIN.GT.0 )
807 $ CALL SCOPY( MNMIN-1, BE, 1, RWORK, 1 )
808 *
809 CALL CBDSQR( UPLO, MNMIN, N, M, NRHS, S2, RWORK, PT,
810 $ LDPT, Q, LDQ, Y, LDX, RWORK( MNMIN+1 ),
811 $ IINFO )
812 *
813 * Test 11: Check the decomposition A := Q*U * S2 * VT*PT
814 * 12: Check the computation Z := U' * Q' * X
815 * 13: Check the orthogonality of Q*U
816 * 14: Check the orthogonality of VT*PT
817 *
818 CALL CBDT01( M, N, 0, A, LDA, Q, LDQ, S2, DUMMA, PT,
819 $ LDPT, WORK, RWORK, RESULT( 11 ) )
820 CALL CBDT02( M, NRHS, X, LDX, Y, LDX, Q, LDQ, WORK,
821 $ RWORK, RESULT( 12 ) )
822 CALL CUNT01( 'Columns', M, MQ, Q, LDQ, WORK, LWORK,
823 $ RWORK, RESULT( 13 ) )
824 CALL CUNT01( 'Rows', MNMIN, N, PT, LDPT, WORK, LWORK,
825 $ RWORK, RESULT( 14 ) )
826 END IF
827 *
828 * End of Loop -- Check for RESULT(j) > THRESH
829 *
830 150 CONTINUE
831 DO 160 J = 1, 14
832 IF( RESULT( J ).GE.THRESH ) THEN
833 IF( NFAIL.EQ.0 )
834 $ CALL SLAHD2( NOUT, PATH )
835 WRITE( NOUT, FMT = 9999 )M, N, JTYPE, IOLDSD, J,
836 $ RESULT( J )
837 NFAIL = NFAIL + 1
838 END IF
839 160 CONTINUE
840 IF( .NOT.BIDIAG ) THEN
841 NTEST = NTEST + 14
842 ELSE
843 NTEST = NTEST + 5
844 END IF
845 *
846 170 CONTINUE
847 180 CONTINUE
848 *
849 * Summary
850 *
851 CALL ALASUM( PATH, NOUT, NFAIL, NTEST, 0 )
852 *
853 RETURN
854 *
855 * End of CCHKBD
856 *
857 9999 FORMAT( ' M=', I5, ', N=', I5, ', type ', I2, ', seed=',
858 $ 4( I4, ',' ), ' test(', I2, ')=', G11.4 )
859 9998 FORMAT( ' CCHKBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
860 $ I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
861 $ I5, ')' )
862 *
863 END