1 SUBROUTINE CDRVBD( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
2 $ A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S,
3 $ SSAV, E, WORK, LWORK, RWORK, IWORK, NOUNIT,
4 $ 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, LDU, LDVT, LWORK, NOUNIT, NSIZES,
12 $ NTYPES
13 REAL THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL DOTYPE( * )
17 INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
18 REAL E( * ), RWORK( * ), S( * ), SSAV( * )
19 COMPLEX A( LDA, * ), ASAV( LDA, * ), U( LDU, * ),
20 $ USAV( LDU, * ), VT( LDVT, * ),
21 $ VTSAV( LDVT, * ), WORK( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * CDRVBD checks the singular value decomposition (SVD) driver CGESVD
28 * and CGESDD.
29 * CGESVD and CGESDD factors A = U diag(S) VT, where U and VT are
30 * unitary and diag(S) is diagonal with the entries of the array S on
31 * its diagonal. The entries of S are the singular values, nonnegative
32 * and stored in decreasing order. U and VT can be optionally not
33 * computed, overwritten on A, or computed partially.
34 *
35 * A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN.
36 * U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N.
37 *
38 * When CDRVBD is called, a number of matrix "sizes" (M's and N's)
39 * and a number of matrix "types" are specified. For each size (M,N)
40 * and each type of matrix, and for the minimal workspace as well as
41 * workspace adequate to permit blocking, an M x N matrix "A" will be
42 * generated and used to test the SVD routines. For each matrix, A will
43 * be factored as A = U diag(S) VT and the following 12 tests computed:
44 *
45 * Test for CGESVD:
46 *
47 * (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
48 *
49 * (2) | I - U'U | / ( M ulp )
50 *
51 * (3) | I - VT VT' | / ( N ulp )
52 *
53 * (4) S contains MNMIN nonnegative values in decreasing order.
54 * (Return 0 if true, 1/ULP if false.)
55 *
56 * (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
57 * computed U.
58 *
59 * (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
60 * computed VT.
61 *
62 * (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
63 * vector of singular values from the partial SVD
64 *
65 * Test for CGESDD:
66 *
67 * (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
68 *
69 * (2) | I - U'U | / ( M ulp )
70 *
71 * (3) | I - VT VT' | / ( N ulp )
72 *
73 * (4) S contains MNMIN nonnegative values in decreasing order.
74 * (Return 0 if true, 1/ULP if false.)
75 *
76 * (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
77 * computed U.
78 *
79 * (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
80 * computed VT.
81 *
82 * (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
83 * vector of singular values from the partial SVD
84 *
85 * The "sizes" are specified by the arrays MM(1:NSIZES) and
86 * NN(1:NSIZES); the value of each element pair (MM(j),NN(j))
87 * specifies one size. The "types" are specified by a logical array
88 * DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j"
89 * will be generated.
90 * Currently, the list of possible types is:
91 *
92 * (1) The zero matrix.
93 * (2) The identity matrix.
94 * (3) A matrix of the form U D V, where U and V are unitary and
95 * D has evenly spaced entries 1, ..., ULP with random signs
96 * on the diagonal.
97 * (4) Same as (3), but multiplied by the underflow-threshold / ULP.
98 * (5) Same as (3), but multiplied by the overflow-threshold * ULP.
99 *
100 * Arguments
101 * ==========
102 *
103 * NSIZES (input) INTEGER
104 * The number of sizes of matrices to use. If it is zero,
105 * CDRVBD does nothing. It must be at least zero.
106 *
107 * MM (input) INTEGER array, dimension (NSIZES)
108 * An array containing the matrix "heights" to be used. For
109 * each j=1,...,NSIZES, if MM(j) is zero, then MM(j) and NN(j)
110 * will be ignored. The MM(j) values must be at least zero.
111 *
112 * NN (input) INTEGER array, dimension (NSIZES)
113 * An array containing the matrix "widths" to be used. For
114 * each j=1,...,NSIZES, if NN(j) is zero, then MM(j) and NN(j)
115 * will be ignored. The NN(j) values must be at least zero.
116 *
117 * NTYPES (input) INTEGER
118 * The number of elements in DOTYPE. If it is zero, CDRVBD
119 * does nothing. It must be at least zero. If it is MAXTYP+1
120 * and NSIZES is 1, then an additional type, MAXTYP+1 is
121 * defined, which is to use whatever matrices are in A and B.
122 * This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
123 * DOTYPE(MAXTYP+1) is .TRUE. .
124 *
125 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
126 * If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
127 * of type j will be generated. If NTYPES is smaller than the
128 * maximum number of types defined (PARAMETER MAXTYP), then
129 * types NTYPES+1 through MAXTYP will not be generated. If
130 * NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
131 * DOTYPE(NTYPES) will be ignored.
132 *
133 * ISEED (input/output) INTEGER array, dimension (4)
134 * On entry ISEED specifies the seed of the random number
135 * generator. The array elements should be between 0 and 4095;
136 * if not they will be reduced mod 4096. Also, ISEED(4) must
137 * be odd. The random number generator uses a linear
138 * congruential sequence limited to small integers, and so
139 * should produce machine independent random numbers. The
140 * values of ISEED are changed on exit, and can be used in the
141 * next call to CDRVBD to continue the same random number
142 * sequence.
143 *
144 * THRESH (input) REAL
145 * A test will count as "failed" if the "error", computed as
146 * described above, exceeds THRESH. Note that the error
147 * is scaled to be O(1), so THRESH should be a reasonably
148 * small multiple of 1, e.g., 10 or 100. In particular,
149 * it should not depend on the precision (single vs. double)
150 * or the size of the matrix. It must be at least zero.
151 *
152 * NOUNIT (input) INTEGER
153 * The FORTRAN unit number for printing out error messages
154 * (e.g., if a routine returns IINFO not equal to 0.)
155 *
156 * A (output) COMPLEX array, dimension (LDA,max(NN))
157 * Used to hold the matrix whose singular values are to be
158 * computed. On exit, A contains the last matrix actually
159 * used.
160 *
161 * LDA (input) INTEGER
162 * The leading dimension of A. It must be at
163 * least 1 and at least max( MM ).
164 *
165 * U (output) COMPLEX array, dimension (LDU,max(MM))
166 * Used to hold the computed matrix of right singular vectors.
167 * On exit, U contains the last such vectors actually computed.
168 *
169 * LDU (input) INTEGER
170 * The leading dimension of U. It must be at
171 * least 1 and at least max( MM ).
172 *
173 * VT (output) COMPLEX array, dimension (LDVT,max(NN))
174 * Used to hold the computed matrix of left singular vectors.
175 * On exit, VT contains the last such vectors actually computed.
176 *
177 * LDVT (input) INTEGER
178 * The leading dimension of VT. It must be at
179 * least 1 and at least max( NN ).
180 *
181 * ASAV (output) COMPLEX array, dimension (LDA,max(NN))
182 * Used to hold a different copy of the matrix whose singular
183 * values are to be computed. On exit, A contains the last
184 * matrix actually used.
185 *
186 * USAV (output) COMPLEX array, dimension (LDU,max(MM))
187 * Used to hold a different copy of the computed matrix of
188 * right singular vectors. On exit, USAV contains the last such
189 * vectors actually computed.
190 *
191 * VTSAV (output) COMPLEX array, dimension (LDVT,max(NN))
192 * Used to hold a different copy of the computed matrix of
193 * left singular vectors. On exit, VTSAV contains the last such
194 * vectors actually computed.
195 *
196 * S (output) REAL array, dimension (max(min(MM,NN)))
197 * Contains the computed singular values.
198 *
199 * SSAV (output) REAL array, dimension (max(min(MM,NN)))
200 * Contains another copy of the computed singular values.
201 *
202 * E (output) REAL array, dimension (max(min(MM,NN)))
203 * Workspace for CGESVD.
204 *
205 * WORK (workspace) COMPLEX array, dimension (LWORK)
206 *
207 * LWORK (input) INTEGER
208 * The number of entries in WORK. This must be at least
209 * MAX(3*MIN(M,N)+MAX(M,N)**2,5*MIN(M,N),3*MAX(M,N)) for all
210 * pairs (M,N)=(MM(j),NN(j))
211 *
212 * RWORK (workspace) REAL array,
213 * dimension ( 5*max(max(MM,NN)) )
214 *
215 * IWORK (workspace) INTEGER array, dimension at least 8*min(M,N)
216 *
217 * RESULT (output) REAL array, dimension (7)
218 * The values computed by the 7 tests described above.
219 * The values are currently limited to 1/ULP, to avoid
220 * overflow.
221 *
222 * INFO (output) INTEGER
223 * If 0, then everything ran OK.
224 * -1: NSIZES < 0
225 * -2: Some MM(j) < 0
226 * -3: Some NN(j) < 0
227 * -4: NTYPES < 0
228 * -7: THRESH < 0
229 * -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
230 * -12: LDU < 1 or LDU < MMAX.
231 * -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ).
232 * -21: LWORK too small.
233 * If CLATMS, or CGESVD returns an error code, the
234 * absolute value of it is returned.
235 *
236 * =====================================================================
237 *
238 * .. Parameters ..
239 REAL ZERO, ONE
240 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
241 COMPLEX CZERO, CONE
242 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
243 $ CONE = ( 1.0E+0, 0.0E+0 ) )
244 INTEGER MAXTYP
245 PARAMETER ( MAXTYP = 5 )
246 * ..
247 * .. Local Scalars ..
248 LOGICAL BADMM, BADNN
249 CHARACTER JOBQ, JOBU, JOBVT
250 INTEGER I, IINFO, IJQ, IJU, IJVT, IWSPC, IWTMP, J,
251 $ JSIZE, JTYPE, LSWORK, M, MINWRK, MMAX, MNMAX,
252 $ MNMIN, MTYPES, N, NERRS, NFAIL, NMAX, NTEST,
253 $ NTESTF, NTESTT
254 REAL ANORM, DIF, DIV, OVFL, ULP, ULPINV, UNFL
255 * ..
256 * .. Local Arrays ..
257 CHARACTER CJOB( 4 )
258 INTEGER IOLDSD( 4 )
259 REAL RESULT( 14 )
260 * ..
261 * .. External Functions ..
262 REAL SLAMCH
263 EXTERNAL SLAMCH
264 * ..
265 * .. External Subroutines ..
266 EXTERNAL ALASVM, CBDT01, CGESDD, CGESVD, CLACPY, CLASET,
267 $ CLATMS, CUNT01, CUNT03, XERBLA
268 * ..
269 * .. Intrinsic Functions ..
270 INTRINSIC ABS, MAX, MIN, REAL
271 * ..
272 * .. Data statements ..
273 DATA CJOB / 'N', 'O', 'S', 'A' /
274 * ..
275 * .. Executable Statements ..
276 *
277 * Check for errors
278 *
279 INFO = 0
280 *
281 * Important constants
282 *
283 NERRS = 0
284 NTESTT = 0
285 NTESTF = 0
286 BADMM = .FALSE.
287 BADNN = .FALSE.
288 MMAX = 1
289 NMAX = 1
290 MNMAX = 1
291 MINWRK = 1
292 DO 10 J = 1, NSIZES
293 MMAX = MAX( MMAX, MM( J ) )
294 IF( MM( J ).LT.0 )
295 $ BADMM = .TRUE.
296 NMAX = MAX( NMAX, NN( J ) )
297 IF( NN( J ).LT.0 )
298 $ BADNN = .TRUE.
299 MNMAX = MAX( MNMAX, MIN( MM( J ), NN( J ) ) )
300 MINWRK = MAX( MINWRK, MAX( 3*MIN( MM( J ),
301 $ NN( J ) )+MAX( MM( J ), NN( J ) )**2, 5*MIN( MM( J ),
302 $ NN( J ) ), 3*MAX( MM( J ), NN( J ) ) ) )
303 10 CONTINUE
304 *
305 * Check for errors
306 *
307 IF( NSIZES.LT.0 ) THEN
308 INFO = -1
309 ELSE IF( BADMM ) THEN
310 INFO = -2
311 ELSE IF( BADNN ) THEN
312 INFO = -3
313 ELSE IF( NTYPES.LT.0 ) THEN
314 INFO = -4
315 ELSE IF( LDA.LT.MAX( 1, MMAX ) ) THEN
316 INFO = -10
317 ELSE IF( LDU.LT.MAX( 1, MMAX ) ) THEN
318 INFO = -12
319 ELSE IF( LDVT.LT.MAX( 1, NMAX ) ) THEN
320 INFO = -14
321 ELSE IF( MINWRK.GT.LWORK ) THEN
322 INFO = -21
323 END IF
324 *
325 IF( INFO.NE.0 ) THEN
326 CALL XERBLA( 'CDRVBD', -INFO )
327 RETURN
328 END IF
329 *
330 * Quick return if nothing to do
331 *
332 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
333 $ RETURN
334 *
335 * More Important constants
336 *
337 UNFL = SLAMCH( 'S' )
338 OVFL = ONE / UNFL
339 ULP = SLAMCH( 'E' )
340 ULPINV = ONE / ULP
341 *
342 * Loop over sizes, types
343 *
344 NERRS = 0
345 *
346 DO 180 JSIZE = 1, NSIZES
347 M = MM( JSIZE )
348 N = NN( JSIZE )
349 MNMIN = MIN( M, N )
350 *
351 IF( NSIZES.NE.1 ) THEN
352 MTYPES = MIN( MAXTYP, NTYPES )
353 ELSE
354 MTYPES = MIN( MAXTYP+1, NTYPES )
355 END IF
356 *
357 DO 170 JTYPE = 1, MTYPES
358 IF( .NOT.DOTYPE( JTYPE ) )
359 $ GO TO 170
360 NTEST = 0
361 *
362 DO 20 J = 1, 4
363 IOLDSD( J ) = ISEED( J )
364 20 CONTINUE
365 *
366 * Compute "A"
367 *
368 IF( MTYPES.GT.MAXTYP )
369 $ GO TO 50
370 *
371 IF( JTYPE.EQ.1 ) THEN
372 *
373 * Zero matrix
374 *
375 CALL CLASET( 'Full', M, N, CZERO, CZERO, A, LDA )
376 DO 30 I = 1, MIN( M, N )
377 S( I ) = ZERO
378 30 CONTINUE
379 *
380 ELSE IF( JTYPE.EQ.2 ) THEN
381 *
382 * Identity matrix
383 *
384 CALL CLASET( 'Full', M, N, CZERO, CONE, A, LDA )
385 DO 40 I = 1, MIN( M, N )
386 S( I ) = ONE
387 40 CONTINUE
388 *
389 ELSE
390 *
391 * (Scaled) random matrix
392 *
393 IF( JTYPE.EQ.3 )
394 $ ANORM = ONE
395 IF( JTYPE.EQ.4 )
396 $ ANORM = UNFL / ULP
397 IF( JTYPE.EQ.5 )
398 $ ANORM = OVFL*ULP
399 CALL CLATMS( M, N, 'U', ISEED, 'N', S, 4, REAL( MNMIN ),
400 $ ANORM, M-1, N-1, 'N', A, LDA, WORK, IINFO )
401 IF( IINFO.NE.0 ) THEN
402 WRITE( NOUNIT, FMT = 9996 )'Generator', IINFO, M, N,
403 $ JTYPE, IOLDSD
404 INFO = ABS( IINFO )
405 RETURN
406 END IF
407 END IF
408 *
409 50 CONTINUE
410 CALL CLACPY( 'F', M, N, A, LDA, ASAV, LDA )
411 *
412 * Do for minimal and adequate (for blocking) workspace
413 *
414 DO 160 IWSPC = 1, 4
415 *
416 * Test for CGESVD
417 *
418 IWTMP = 2*MIN( M, N )+MAX( M, N )
419 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
420 LSWORK = MIN( LSWORK, LWORK )
421 LSWORK = MAX( LSWORK, 1 )
422 IF( IWSPC.EQ.4 )
423 $ LSWORK = LWORK
424 *
425 DO 60 J = 1, 14
426 RESULT( J ) = -ONE
427 60 CONTINUE
428 *
429 * Factorize A
430 *
431 IF( IWSPC.GT.1 )
432 $ CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
433 CALL CGESVD( 'A', 'A', M, N, A, LDA, SSAV, USAV, LDU,
434 $ VTSAV, LDVT, WORK, LSWORK, RWORK, IINFO )
435 IF( IINFO.NE.0 ) THEN
436 WRITE( NOUNIT, FMT = 9995 )'GESVD', IINFO, M, N,
437 $ JTYPE, LSWORK, IOLDSD
438 INFO = ABS( IINFO )
439 RETURN
440 END IF
441 *
442 * Do tests 1--4
443 *
444 CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
445 $ VTSAV, LDVT, WORK, RWORK, RESULT( 1 ) )
446 IF( M.NE.0 .AND. N.NE.0 ) THEN
447 CALL CUNT01( 'Columns', MNMIN, M, USAV, LDU, WORK,
448 $ LWORK, RWORK, RESULT( 2 ) )
449 CALL CUNT01( 'Rows', MNMIN, N, VTSAV, LDVT, WORK,
450 $ LWORK, RWORK, RESULT( 3 ) )
451 END IF
452 RESULT( 4 ) = 0
453 DO 70 I = 1, MNMIN - 1
454 IF( SSAV( I ).LT.SSAV( I+1 ) )
455 $ RESULT( 4 ) = ULPINV
456 IF( SSAV( I ).LT.ZERO )
457 $ RESULT( 4 ) = ULPINV
458 70 CONTINUE
459 IF( MNMIN.GE.1 ) THEN
460 IF( SSAV( MNMIN ).LT.ZERO )
461 $ RESULT( 4 ) = ULPINV
462 END IF
463 *
464 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
465 *
466 RESULT( 5 ) = ZERO
467 RESULT( 6 ) = ZERO
468 RESULT( 7 ) = ZERO
469 DO 100 IJU = 0, 3
470 DO 90 IJVT = 0, 3
471 IF( ( IJU.EQ.3 .AND. IJVT.EQ.3 ) .OR.
472 $ ( IJU.EQ.1 .AND. IJVT.EQ.1 ) )GO TO 90
473 JOBU = CJOB( IJU+1 )
474 JOBVT = CJOB( IJVT+1 )
475 CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
476 CALL CGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
477 $ VT, LDVT, WORK, LSWORK, RWORK, IINFO )
478 *
479 * Compare U
480 *
481 DIF = ZERO
482 IF( M.GT.0 .AND. N.GT.0 ) THEN
483 IF( IJU.EQ.1 ) THEN
484 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
485 $ LDU, A, LDA, WORK, LWORK, RWORK,
486 $ DIF, IINFO )
487 ELSE IF( IJU.EQ.2 ) THEN
488 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
489 $ LDU, U, LDU, WORK, LWORK, RWORK,
490 $ DIF, IINFO )
491 ELSE IF( IJU.EQ.3 ) THEN
492 CALL CUNT03( 'C', M, M, M, MNMIN, USAV, LDU,
493 $ U, LDU, WORK, LWORK, RWORK, DIF,
494 $ IINFO )
495 END IF
496 END IF
497 RESULT( 5 ) = MAX( RESULT( 5 ), DIF )
498 *
499 * Compare VT
500 *
501 DIF = ZERO
502 IF( M.GT.0 .AND. N.GT.0 ) THEN
503 IF( IJVT.EQ.1 ) THEN
504 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
505 $ LDVT, A, LDA, WORK, LWORK,
506 $ RWORK, DIF, IINFO )
507 ELSE IF( IJVT.EQ.2 ) THEN
508 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
509 $ LDVT, VT, LDVT, WORK, LWORK,
510 $ RWORK, DIF, IINFO )
511 ELSE IF( IJVT.EQ.3 ) THEN
512 CALL CUNT03( 'R', N, N, N, MNMIN, VTSAV,
513 $ LDVT, VT, LDVT, WORK, LWORK,
514 $ RWORK, DIF, IINFO )
515 END IF
516 END IF
517 RESULT( 6 ) = MAX( RESULT( 6 ), DIF )
518 *
519 * Compare S
520 *
521 DIF = ZERO
522 DIV = MAX( REAL( MNMIN )*ULP*S( 1 ),
523 $ SLAMCH( 'Safe minimum' ) )
524 DO 80 I = 1, MNMIN - 1
525 IF( SSAV( I ).LT.SSAV( I+1 ) )
526 $ DIF = ULPINV
527 IF( SSAV( I ).LT.ZERO )
528 $ DIF = ULPINV
529 DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
530 80 CONTINUE
531 RESULT( 7 ) = MAX( RESULT( 7 ), DIF )
532 90 CONTINUE
533 100 CONTINUE
534 *
535 * Test for CGESDD
536 *
537 IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
538 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
539 LSWORK = MIN( LSWORK, LWORK )
540 LSWORK = MAX( LSWORK, 1 )
541 IF( IWSPC.EQ.4 )
542 $ LSWORK = LWORK
543 *
544 * Factorize A
545 *
546 CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
547 CALL CGESDD( 'A', M, N, A, LDA, SSAV, USAV, LDU, VTSAV,
548 $ LDVT, WORK, LSWORK, RWORK, IWORK, IINFO )
549 IF( IINFO.NE.0 ) THEN
550 WRITE( NOUNIT, FMT = 9995 )'GESDD', IINFO, M, N,
551 $ JTYPE, LSWORK, IOLDSD
552 INFO = ABS( IINFO )
553 RETURN
554 END IF
555 *
556 * Do tests 1--4
557 *
558 CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
559 $ VTSAV, LDVT, WORK, RWORK, RESULT( 8 ) )
560 IF( M.NE.0 .AND. N.NE.0 ) THEN
561 CALL CUNT01( 'Columns', MNMIN, M, USAV, LDU, WORK,
562 $ LWORK, RWORK, RESULT( 9 ) )
563 CALL CUNT01( 'Rows', MNMIN, N, VTSAV, LDVT, WORK,
564 $ LWORK, RWORK, RESULT( 10 ) )
565 END IF
566 RESULT( 11 ) = 0
567 DO 110 I = 1, MNMIN - 1
568 IF( SSAV( I ).LT.SSAV( I+1 ) )
569 $ RESULT( 11 ) = ULPINV
570 IF( SSAV( I ).LT.ZERO )
571 $ RESULT( 11 ) = ULPINV
572 110 CONTINUE
573 IF( MNMIN.GE.1 ) THEN
574 IF( SSAV( MNMIN ).LT.ZERO )
575 $ RESULT( 11 ) = ULPINV
576 END IF
577 *
578 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
579 *
580 RESULT( 12 ) = ZERO
581 RESULT( 13 ) = ZERO
582 RESULT( 14 ) = ZERO
583 DO 130 IJQ = 0, 2
584 JOBQ = CJOB( IJQ+1 )
585 CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
586 CALL CGESDD( JOBQ, M, N, A, LDA, S, U, LDU, VT, LDVT,
587 $ WORK, LSWORK, RWORK, IWORK, IINFO )
588 *
589 * Compare U
590 *
591 DIF = ZERO
592 IF( M.GT.0 .AND. N.GT.0 ) THEN
593 IF( IJQ.EQ.1 ) THEN
594 IF( M.GE.N ) THEN
595 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
596 $ LDU, A, LDA, WORK, LWORK, RWORK,
597 $ DIF, IINFO )
598 ELSE
599 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
600 $ LDU, U, LDU, WORK, LWORK, RWORK,
601 $ DIF, IINFO )
602 END IF
603 ELSE IF( IJQ.EQ.2 ) THEN
604 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV, LDU,
605 $ U, LDU, WORK, LWORK, RWORK, DIF,
606 $ IINFO )
607 END IF
608 END IF
609 RESULT( 12 ) = MAX( RESULT( 12 ), DIF )
610 *
611 * Compare VT
612 *
613 DIF = ZERO
614 IF( M.GT.0 .AND. N.GT.0 ) THEN
615 IF( IJQ.EQ.1 ) THEN
616 IF( M.GE.N ) THEN
617 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
618 $ LDVT, VT, LDVT, WORK, LWORK,
619 $ RWORK, DIF, IINFO )
620 ELSE
621 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
622 $ LDVT, A, LDA, WORK, LWORK,
623 $ RWORK, DIF, IINFO )
624 END IF
625 ELSE IF( IJQ.EQ.2 ) THEN
626 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
627 $ LDVT, VT, LDVT, WORK, LWORK, RWORK,
628 $ DIF, IINFO )
629 END IF
630 END IF
631 RESULT( 13 ) = MAX( RESULT( 13 ), DIF )
632 *
633 * Compare S
634 *
635 DIF = ZERO
636 DIV = MAX( REAL( MNMIN )*ULP*S( 1 ),
637 $ SLAMCH( 'Safe minimum' ) )
638 DO 120 I = 1, MNMIN - 1
639 IF( SSAV( I ).LT.SSAV( I+1 ) )
640 $ DIF = ULPINV
641 IF( SSAV( I ).LT.ZERO )
642 $ DIF = ULPINV
643 DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
644 120 CONTINUE
645 RESULT( 14 ) = MAX( RESULT( 14 ), DIF )
646 130 CONTINUE
647 *
648 * End of Loop -- Check for RESULT(j) > THRESH
649 *
650 NTEST = 0
651 NFAIL = 0
652 DO 140 J = 1, 14
653 IF( RESULT( J ).GE.ZERO )
654 $ NTEST = NTEST + 1
655 IF( RESULT( J ).GE.THRESH )
656 $ NFAIL = NFAIL + 1
657 140 CONTINUE
658 *
659 IF( NFAIL.GT.0 )
660 $ NTESTF = NTESTF + 1
661 IF( NTESTF.EQ.1 ) THEN
662 WRITE( NOUNIT, FMT = 9999 )
663 WRITE( NOUNIT, FMT = 9998 )THRESH
664 NTESTF = 2
665 END IF
666 *
667 DO 150 J = 1, 14
668 IF( RESULT( J ).GE.THRESH ) THEN
669 WRITE( NOUNIT, FMT = 9997 )M, N, JTYPE, IWSPC,
670 $ IOLDSD, J, RESULT( J )
671 END IF
672 150 CONTINUE
673 *
674 NERRS = NERRS + NFAIL
675 NTESTT = NTESTT + NTEST
676 *
677 160 CONTINUE
678 *
679 170 CONTINUE
680 180 CONTINUE
681 *
682 * Summary
683 *
684 CALL ALASVM( 'CBD', NOUNIT, NERRS, NTESTT, 0 )
685 *
686 9999 FORMAT( ' SVD -- Complex Singular Value Decomposition Driver ',
687 $ / ' Matrix types (see CDRVBD for details):',
688 $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
689 $ / ' 3 = Evenly spaced singular values near 1',
690 $ / ' 4 = Evenly spaced singular values near underflow',
691 $ / ' 5 = Evenly spaced singular values near overflow',
692 $ / / ' Tests performed: ( A is dense, U and V are unitary,',
693 $ / 19X, ' S is an array, and Upartial, VTpartial, and',
694 $ / 19X, ' Spartial are partially computed U, VT and S),', / )
695 9998 FORMAT( ' Tests performed with Test Threshold = ', F8.2,
696 $ / ' CGESVD: ', /
697 $ ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
698 $ / ' 2 = | I - U**T U | / ( M ulp ) ',
699 $ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
700 $ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
701 $ ' decreasing order, else 1/ulp',
702 $ / ' 5 = | U - Upartial | / ( M ulp )',
703 $ / ' 6 = | VT - VTpartial | / ( N ulp )',
704 $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
705 $ / ' CGESDD: ', /
706 $ ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
707 $ / ' 9 = | I - U**T U | / ( M ulp ) ',
708 $ / '10 = | I - VT VT**T | / ( N ulp ) ',
709 $ / '11 = 0 if S contains min(M,N) nonnegative values in',
710 $ ' decreasing order, else 1/ulp',
711 $ / '12 = | U - Upartial | / ( M ulp )',
712 $ / '13 = | VT - VTpartial | / ( N ulp )',
713 $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )', / / )
714 9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1,
715 $ ', seed=', 4( I4, ',' ), ' test(', I1, ')=', G11.4 )
716 9996 FORMAT( ' CDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
717 $ I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
718 $ I5, ')' )
719 9995 FORMAT( ' CDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
720 $ I6, ', N=', I6, ', JTYPE=', I6, ', LSWORK=', I6, / 9X,
721 $ 'ISEED=(', 3( I5, ',' ), I5, ')' )
722 *
723 RETURN
724 *
725 * End of CDRVBD
726 *
727 END
2 $ A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S,
3 $ SSAV, E, WORK, LWORK, RWORK, IWORK, NOUNIT,
4 $ 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, LDU, LDVT, LWORK, NOUNIT, NSIZES,
12 $ NTYPES
13 REAL THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL DOTYPE( * )
17 INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
18 REAL E( * ), RWORK( * ), S( * ), SSAV( * )
19 COMPLEX A( LDA, * ), ASAV( LDA, * ), U( LDU, * ),
20 $ USAV( LDU, * ), VT( LDVT, * ),
21 $ VTSAV( LDVT, * ), WORK( * )
22 * ..
23 *
24 * Purpose
25 * =======
26 *
27 * CDRVBD checks the singular value decomposition (SVD) driver CGESVD
28 * and CGESDD.
29 * CGESVD and CGESDD factors A = U diag(S) VT, where U and VT are
30 * unitary and diag(S) is diagonal with the entries of the array S on
31 * its diagonal. The entries of S are the singular values, nonnegative
32 * and stored in decreasing order. U and VT can be optionally not
33 * computed, overwritten on A, or computed partially.
34 *
35 * A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN.
36 * U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N.
37 *
38 * When CDRVBD is called, a number of matrix "sizes" (M's and N's)
39 * and a number of matrix "types" are specified. For each size (M,N)
40 * and each type of matrix, and for the minimal workspace as well as
41 * workspace adequate to permit blocking, an M x N matrix "A" will be
42 * generated and used to test the SVD routines. For each matrix, A will
43 * be factored as A = U diag(S) VT and the following 12 tests computed:
44 *
45 * Test for CGESVD:
46 *
47 * (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
48 *
49 * (2) | I - U'U | / ( M ulp )
50 *
51 * (3) | I - VT VT' | / ( N ulp )
52 *
53 * (4) S contains MNMIN nonnegative values in decreasing order.
54 * (Return 0 if true, 1/ULP if false.)
55 *
56 * (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
57 * computed U.
58 *
59 * (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
60 * computed VT.
61 *
62 * (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
63 * vector of singular values from the partial SVD
64 *
65 * Test for CGESDD:
66 *
67 * (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
68 *
69 * (2) | I - U'U | / ( M ulp )
70 *
71 * (3) | I - VT VT' | / ( N ulp )
72 *
73 * (4) S contains MNMIN nonnegative values in decreasing order.
74 * (Return 0 if true, 1/ULP if false.)
75 *
76 * (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
77 * computed U.
78 *
79 * (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
80 * computed VT.
81 *
82 * (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
83 * vector of singular values from the partial SVD
84 *
85 * The "sizes" are specified by the arrays MM(1:NSIZES) and
86 * NN(1:NSIZES); the value of each element pair (MM(j),NN(j))
87 * specifies one size. The "types" are specified by a logical array
88 * DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j"
89 * will be generated.
90 * Currently, the list of possible types is:
91 *
92 * (1) The zero matrix.
93 * (2) The identity matrix.
94 * (3) A matrix of the form U D V, where U and V are unitary and
95 * D has evenly spaced entries 1, ..., ULP with random signs
96 * on the diagonal.
97 * (4) Same as (3), but multiplied by the underflow-threshold / ULP.
98 * (5) Same as (3), but multiplied by the overflow-threshold * ULP.
99 *
100 * Arguments
101 * ==========
102 *
103 * NSIZES (input) INTEGER
104 * The number of sizes of matrices to use. If it is zero,
105 * CDRVBD does nothing. It must be at least zero.
106 *
107 * MM (input) INTEGER array, dimension (NSIZES)
108 * An array containing the matrix "heights" to be used. For
109 * each j=1,...,NSIZES, if MM(j) is zero, then MM(j) and NN(j)
110 * will be ignored. The MM(j) values must be at least zero.
111 *
112 * NN (input) INTEGER array, dimension (NSIZES)
113 * An array containing the matrix "widths" to be used. For
114 * each j=1,...,NSIZES, if NN(j) is zero, then MM(j) and NN(j)
115 * will be ignored. The NN(j) values must be at least zero.
116 *
117 * NTYPES (input) INTEGER
118 * The number of elements in DOTYPE. If it is zero, CDRVBD
119 * does nothing. It must be at least zero. If it is MAXTYP+1
120 * and NSIZES is 1, then an additional type, MAXTYP+1 is
121 * defined, which is to use whatever matrices are in A and B.
122 * This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
123 * DOTYPE(MAXTYP+1) is .TRUE. .
124 *
125 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
126 * If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
127 * of type j will be generated. If NTYPES is smaller than the
128 * maximum number of types defined (PARAMETER MAXTYP), then
129 * types NTYPES+1 through MAXTYP will not be generated. If
130 * NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
131 * DOTYPE(NTYPES) will be ignored.
132 *
133 * ISEED (input/output) INTEGER array, dimension (4)
134 * On entry ISEED specifies the seed of the random number
135 * generator. The array elements should be between 0 and 4095;
136 * if not they will be reduced mod 4096. Also, ISEED(4) must
137 * be odd. The random number generator uses a linear
138 * congruential sequence limited to small integers, and so
139 * should produce machine independent random numbers. The
140 * values of ISEED are changed on exit, and can be used in the
141 * next call to CDRVBD to continue the same random number
142 * sequence.
143 *
144 * THRESH (input) REAL
145 * A test will count as "failed" if the "error", computed as
146 * described above, exceeds THRESH. Note that the error
147 * is scaled to be O(1), so THRESH should be a reasonably
148 * small multiple of 1, e.g., 10 or 100. In particular,
149 * it should not depend on the precision (single vs. double)
150 * or the size of the matrix. It must be at least zero.
151 *
152 * NOUNIT (input) INTEGER
153 * The FORTRAN unit number for printing out error messages
154 * (e.g., if a routine returns IINFO not equal to 0.)
155 *
156 * A (output) COMPLEX array, dimension (LDA,max(NN))
157 * Used to hold the matrix whose singular values are to be
158 * computed. On exit, A contains the last matrix actually
159 * used.
160 *
161 * LDA (input) INTEGER
162 * The leading dimension of A. It must be at
163 * least 1 and at least max( MM ).
164 *
165 * U (output) COMPLEX array, dimension (LDU,max(MM))
166 * Used to hold the computed matrix of right singular vectors.
167 * On exit, U contains the last such vectors actually computed.
168 *
169 * LDU (input) INTEGER
170 * The leading dimension of U. It must be at
171 * least 1 and at least max( MM ).
172 *
173 * VT (output) COMPLEX array, dimension (LDVT,max(NN))
174 * Used to hold the computed matrix of left singular vectors.
175 * On exit, VT contains the last such vectors actually computed.
176 *
177 * LDVT (input) INTEGER
178 * The leading dimension of VT. It must be at
179 * least 1 and at least max( NN ).
180 *
181 * ASAV (output) COMPLEX array, dimension (LDA,max(NN))
182 * Used to hold a different copy of the matrix whose singular
183 * values are to be computed. On exit, A contains the last
184 * matrix actually used.
185 *
186 * USAV (output) COMPLEX array, dimension (LDU,max(MM))
187 * Used to hold a different copy of the computed matrix of
188 * right singular vectors. On exit, USAV contains the last such
189 * vectors actually computed.
190 *
191 * VTSAV (output) COMPLEX array, dimension (LDVT,max(NN))
192 * Used to hold a different copy of the computed matrix of
193 * left singular vectors. On exit, VTSAV contains the last such
194 * vectors actually computed.
195 *
196 * S (output) REAL array, dimension (max(min(MM,NN)))
197 * Contains the computed singular values.
198 *
199 * SSAV (output) REAL array, dimension (max(min(MM,NN)))
200 * Contains another copy of the computed singular values.
201 *
202 * E (output) REAL array, dimension (max(min(MM,NN)))
203 * Workspace for CGESVD.
204 *
205 * WORK (workspace) COMPLEX array, dimension (LWORK)
206 *
207 * LWORK (input) INTEGER
208 * The number of entries in WORK. This must be at least
209 * MAX(3*MIN(M,N)+MAX(M,N)**2,5*MIN(M,N),3*MAX(M,N)) for all
210 * pairs (M,N)=(MM(j),NN(j))
211 *
212 * RWORK (workspace) REAL array,
213 * dimension ( 5*max(max(MM,NN)) )
214 *
215 * IWORK (workspace) INTEGER array, dimension at least 8*min(M,N)
216 *
217 * RESULT (output) REAL array, dimension (7)
218 * The values computed by the 7 tests described above.
219 * The values are currently limited to 1/ULP, to avoid
220 * overflow.
221 *
222 * INFO (output) INTEGER
223 * If 0, then everything ran OK.
224 * -1: NSIZES < 0
225 * -2: Some MM(j) < 0
226 * -3: Some NN(j) < 0
227 * -4: NTYPES < 0
228 * -7: THRESH < 0
229 * -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
230 * -12: LDU < 1 or LDU < MMAX.
231 * -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ).
232 * -21: LWORK too small.
233 * If CLATMS, or CGESVD returns an error code, the
234 * absolute value of it is returned.
235 *
236 * =====================================================================
237 *
238 * .. Parameters ..
239 REAL ZERO, ONE
240 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
241 COMPLEX CZERO, CONE
242 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
243 $ CONE = ( 1.0E+0, 0.0E+0 ) )
244 INTEGER MAXTYP
245 PARAMETER ( MAXTYP = 5 )
246 * ..
247 * .. Local Scalars ..
248 LOGICAL BADMM, BADNN
249 CHARACTER JOBQ, JOBU, JOBVT
250 INTEGER I, IINFO, IJQ, IJU, IJVT, IWSPC, IWTMP, J,
251 $ JSIZE, JTYPE, LSWORK, M, MINWRK, MMAX, MNMAX,
252 $ MNMIN, MTYPES, N, NERRS, NFAIL, NMAX, NTEST,
253 $ NTESTF, NTESTT
254 REAL ANORM, DIF, DIV, OVFL, ULP, ULPINV, UNFL
255 * ..
256 * .. Local Arrays ..
257 CHARACTER CJOB( 4 )
258 INTEGER IOLDSD( 4 )
259 REAL RESULT( 14 )
260 * ..
261 * .. External Functions ..
262 REAL SLAMCH
263 EXTERNAL SLAMCH
264 * ..
265 * .. External Subroutines ..
266 EXTERNAL ALASVM, CBDT01, CGESDD, CGESVD, CLACPY, CLASET,
267 $ CLATMS, CUNT01, CUNT03, XERBLA
268 * ..
269 * .. Intrinsic Functions ..
270 INTRINSIC ABS, MAX, MIN, REAL
271 * ..
272 * .. Data statements ..
273 DATA CJOB / 'N', 'O', 'S', 'A' /
274 * ..
275 * .. Executable Statements ..
276 *
277 * Check for errors
278 *
279 INFO = 0
280 *
281 * Important constants
282 *
283 NERRS = 0
284 NTESTT = 0
285 NTESTF = 0
286 BADMM = .FALSE.
287 BADNN = .FALSE.
288 MMAX = 1
289 NMAX = 1
290 MNMAX = 1
291 MINWRK = 1
292 DO 10 J = 1, NSIZES
293 MMAX = MAX( MMAX, MM( J ) )
294 IF( MM( J ).LT.0 )
295 $ BADMM = .TRUE.
296 NMAX = MAX( NMAX, NN( J ) )
297 IF( NN( J ).LT.0 )
298 $ BADNN = .TRUE.
299 MNMAX = MAX( MNMAX, MIN( MM( J ), NN( J ) ) )
300 MINWRK = MAX( MINWRK, MAX( 3*MIN( MM( J ),
301 $ NN( J ) )+MAX( MM( J ), NN( J ) )**2, 5*MIN( MM( J ),
302 $ NN( J ) ), 3*MAX( MM( J ), NN( J ) ) ) )
303 10 CONTINUE
304 *
305 * Check for errors
306 *
307 IF( NSIZES.LT.0 ) THEN
308 INFO = -1
309 ELSE IF( BADMM ) THEN
310 INFO = -2
311 ELSE IF( BADNN ) THEN
312 INFO = -3
313 ELSE IF( NTYPES.LT.0 ) THEN
314 INFO = -4
315 ELSE IF( LDA.LT.MAX( 1, MMAX ) ) THEN
316 INFO = -10
317 ELSE IF( LDU.LT.MAX( 1, MMAX ) ) THEN
318 INFO = -12
319 ELSE IF( LDVT.LT.MAX( 1, NMAX ) ) THEN
320 INFO = -14
321 ELSE IF( MINWRK.GT.LWORK ) THEN
322 INFO = -21
323 END IF
324 *
325 IF( INFO.NE.0 ) THEN
326 CALL XERBLA( 'CDRVBD', -INFO )
327 RETURN
328 END IF
329 *
330 * Quick return if nothing to do
331 *
332 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
333 $ RETURN
334 *
335 * More Important constants
336 *
337 UNFL = SLAMCH( 'S' )
338 OVFL = ONE / UNFL
339 ULP = SLAMCH( 'E' )
340 ULPINV = ONE / ULP
341 *
342 * Loop over sizes, types
343 *
344 NERRS = 0
345 *
346 DO 180 JSIZE = 1, NSIZES
347 M = MM( JSIZE )
348 N = NN( JSIZE )
349 MNMIN = MIN( M, N )
350 *
351 IF( NSIZES.NE.1 ) THEN
352 MTYPES = MIN( MAXTYP, NTYPES )
353 ELSE
354 MTYPES = MIN( MAXTYP+1, NTYPES )
355 END IF
356 *
357 DO 170 JTYPE = 1, MTYPES
358 IF( .NOT.DOTYPE( JTYPE ) )
359 $ GO TO 170
360 NTEST = 0
361 *
362 DO 20 J = 1, 4
363 IOLDSD( J ) = ISEED( J )
364 20 CONTINUE
365 *
366 * Compute "A"
367 *
368 IF( MTYPES.GT.MAXTYP )
369 $ GO TO 50
370 *
371 IF( JTYPE.EQ.1 ) THEN
372 *
373 * Zero matrix
374 *
375 CALL CLASET( 'Full', M, N, CZERO, CZERO, A, LDA )
376 DO 30 I = 1, MIN( M, N )
377 S( I ) = ZERO
378 30 CONTINUE
379 *
380 ELSE IF( JTYPE.EQ.2 ) THEN
381 *
382 * Identity matrix
383 *
384 CALL CLASET( 'Full', M, N, CZERO, CONE, A, LDA )
385 DO 40 I = 1, MIN( M, N )
386 S( I ) = ONE
387 40 CONTINUE
388 *
389 ELSE
390 *
391 * (Scaled) random matrix
392 *
393 IF( JTYPE.EQ.3 )
394 $ ANORM = ONE
395 IF( JTYPE.EQ.4 )
396 $ ANORM = UNFL / ULP
397 IF( JTYPE.EQ.5 )
398 $ ANORM = OVFL*ULP
399 CALL CLATMS( M, N, 'U', ISEED, 'N', S, 4, REAL( MNMIN ),
400 $ ANORM, M-1, N-1, 'N', A, LDA, WORK, IINFO )
401 IF( IINFO.NE.0 ) THEN
402 WRITE( NOUNIT, FMT = 9996 )'Generator', IINFO, M, N,
403 $ JTYPE, IOLDSD
404 INFO = ABS( IINFO )
405 RETURN
406 END IF
407 END IF
408 *
409 50 CONTINUE
410 CALL CLACPY( 'F', M, N, A, LDA, ASAV, LDA )
411 *
412 * Do for minimal and adequate (for blocking) workspace
413 *
414 DO 160 IWSPC = 1, 4
415 *
416 * Test for CGESVD
417 *
418 IWTMP = 2*MIN( M, N )+MAX( M, N )
419 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
420 LSWORK = MIN( LSWORK, LWORK )
421 LSWORK = MAX( LSWORK, 1 )
422 IF( IWSPC.EQ.4 )
423 $ LSWORK = LWORK
424 *
425 DO 60 J = 1, 14
426 RESULT( J ) = -ONE
427 60 CONTINUE
428 *
429 * Factorize A
430 *
431 IF( IWSPC.GT.1 )
432 $ CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
433 CALL CGESVD( 'A', 'A', M, N, A, LDA, SSAV, USAV, LDU,
434 $ VTSAV, LDVT, WORK, LSWORK, RWORK, IINFO )
435 IF( IINFO.NE.0 ) THEN
436 WRITE( NOUNIT, FMT = 9995 )'GESVD', IINFO, M, N,
437 $ JTYPE, LSWORK, IOLDSD
438 INFO = ABS( IINFO )
439 RETURN
440 END IF
441 *
442 * Do tests 1--4
443 *
444 CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
445 $ VTSAV, LDVT, WORK, RWORK, RESULT( 1 ) )
446 IF( M.NE.0 .AND. N.NE.0 ) THEN
447 CALL CUNT01( 'Columns', MNMIN, M, USAV, LDU, WORK,
448 $ LWORK, RWORK, RESULT( 2 ) )
449 CALL CUNT01( 'Rows', MNMIN, N, VTSAV, LDVT, WORK,
450 $ LWORK, RWORK, RESULT( 3 ) )
451 END IF
452 RESULT( 4 ) = 0
453 DO 70 I = 1, MNMIN - 1
454 IF( SSAV( I ).LT.SSAV( I+1 ) )
455 $ RESULT( 4 ) = ULPINV
456 IF( SSAV( I ).LT.ZERO )
457 $ RESULT( 4 ) = ULPINV
458 70 CONTINUE
459 IF( MNMIN.GE.1 ) THEN
460 IF( SSAV( MNMIN ).LT.ZERO )
461 $ RESULT( 4 ) = ULPINV
462 END IF
463 *
464 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
465 *
466 RESULT( 5 ) = ZERO
467 RESULT( 6 ) = ZERO
468 RESULT( 7 ) = ZERO
469 DO 100 IJU = 0, 3
470 DO 90 IJVT = 0, 3
471 IF( ( IJU.EQ.3 .AND. IJVT.EQ.3 ) .OR.
472 $ ( IJU.EQ.1 .AND. IJVT.EQ.1 ) )GO TO 90
473 JOBU = CJOB( IJU+1 )
474 JOBVT = CJOB( IJVT+1 )
475 CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
476 CALL CGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
477 $ VT, LDVT, WORK, LSWORK, RWORK, IINFO )
478 *
479 * Compare U
480 *
481 DIF = ZERO
482 IF( M.GT.0 .AND. N.GT.0 ) THEN
483 IF( IJU.EQ.1 ) THEN
484 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
485 $ LDU, A, LDA, WORK, LWORK, RWORK,
486 $ DIF, IINFO )
487 ELSE IF( IJU.EQ.2 ) THEN
488 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
489 $ LDU, U, LDU, WORK, LWORK, RWORK,
490 $ DIF, IINFO )
491 ELSE IF( IJU.EQ.3 ) THEN
492 CALL CUNT03( 'C', M, M, M, MNMIN, USAV, LDU,
493 $ U, LDU, WORK, LWORK, RWORK, DIF,
494 $ IINFO )
495 END IF
496 END IF
497 RESULT( 5 ) = MAX( RESULT( 5 ), DIF )
498 *
499 * Compare VT
500 *
501 DIF = ZERO
502 IF( M.GT.0 .AND. N.GT.0 ) THEN
503 IF( IJVT.EQ.1 ) THEN
504 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
505 $ LDVT, A, LDA, WORK, LWORK,
506 $ RWORK, DIF, IINFO )
507 ELSE IF( IJVT.EQ.2 ) THEN
508 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
509 $ LDVT, VT, LDVT, WORK, LWORK,
510 $ RWORK, DIF, IINFO )
511 ELSE IF( IJVT.EQ.3 ) THEN
512 CALL CUNT03( 'R', N, N, N, MNMIN, VTSAV,
513 $ LDVT, VT, LDVT, WORK, LWORK,
514 $ RWORK, DIF, IINFO )
515 END IF
516 END IF
517 RESULT( 6 ) = MAX( RESULT( 6 ), DIF )
518 *
519 * Compare S
520 *
521 DIF = ZERO
522 DIV = MAX( REAL( MNMIN )*ULP*S( 1 ),
523 $ SLAMCH( 'Safe minimum' ) )
524 DO 80 I = 1, MNMIN - 1
525 IF( SSAV( I ).LT.SSAV( I+1 ) )
526 $ DIF = ULPINV
527 IF( SSAV( I ).LT.ZERO )
528 $ DIF = ULPINV
529 DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
530 80 CONTINUE
531 RESULT( 7 ) = MAX( RESULT( 7 ), DIF )
532 90 CONTINUE
533 100 CONTINUE
534 *
535 * Test for CGESDD
536 *
537 IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
538 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
539 LSWORK = MIN( LSWORK, LWORK )
540 LSWORK = MAX( LSWORK, 1 )
541 IF( IWSPC.EQ.4 )
542 $ LSWORK = LWORK
543 *
544 * Factorize A
545 *
546 CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
547 CALL CGESDD( 'A', M, N, A, LDA, SSAV, USAV, LDU, VTSAV,
548 $ LDVT, WORK, LSWORK, RWORK, IWORK, IINFO )
549 IF( IINFO.NE.0 ) THEN
550 WRITE( NOUNIT, FMT = 9995 )'GESDD', IINFO, M, N,
551 $ JTYPE, LSWORK, IOLDSD
552 INFO = ABS( IINFO )
553 RETURN
554 END IF
555 *
556 * Do tests 1--4
557 *
558 CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
559 $ VTSAV, LDVT, WORK, RWORK, RESULT( 8 ) )
560 IF( M.NE.0 .AND. N.NE.0 ) THEN
561 CALL CUNT01( 'Columns', MNMIN, M, USAV, LDU, WORK,
562 $ LWORK, RWORK, RESULT( 9 ) )
563 CALL CUNT01( 'Rows', MNMIN, N, VTSAV, LDVT, WORK,
564 $ LWORK, RWORK, RESULT( 10 ) )
565 END IF
566 RESULT( 11 ) = 0
567 DO 110 I = 1, MNMIN - 1
568 IF( SSAV( I ).LT.SSAV( I+1 ) )
569 $ RESULT( 11 ) = ULPINV
570 IF( SSAV( I ).LT.ZERO )
571 $ RESULT( 11 ) = ULPINV
572 110 CONTINUE
573 IF( MNMIN.GE.1 ) THEN
574 IF( SSAV( MNMIN ).LT.ZERO )
575 $ RESULT( 11 ) = ULPINV
576 END IF
577 *
578 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
579 *
580 RESULT( 12 ) = ZERO
581 RESULT( 13 ) = ZERO
582 RESULT( 14 ) = ZERO
583 DO 130 IJQ = 0, 2
584 JOBQ = CJOB( IJQ+1 )
585 CALL CLACPY( 'F', M, N, ASAV, LDA, A, LDA )
586 CALL CGESDD( JOBQ, M, N, A, LDA, S, U, LDU, VT, LDVT,
587 $ WORK, LSWORK, RWORK, IWORK, IINFO )
588 *
589 * Compare U
590 *
591 DIF = ZERO
592 IF( M.GT.0 .AND. N.GT.0 ) THEN
593 IF( IJQ.EQ.1 ) THEN
594 IF( M.GE.N ) THEN
595 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
596 $ LDU, A, LDA, WORK, LWORK, RWORK,
597 $ DIF, IINFO )
598 ELSE
599 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV,
600 $ LDU, U, LDU, WORK, LWORK, RWORK,
601 $ DIF, IINFO )
602 END IF
603 ELSE IF( IJQ.EQ.2 ) THEN
604 CALL CUNT03( 'C', M, MNMIN, M, MNMIN, USAV, LDU,
605 $ U, LDU, WORK, LWORK, RWORK, DIF,
606 $ IINFO )
607 END IF
608 END IF
609 RESULT( 12 ) = MAX( RESULT( 12 ), DIF )
610 *
611 * Compare VT
612 *
613 DIF = ZERO
614 IF( M.GT.0 .AND. N.GT.0 ) THEN
615 IF( IJQ.EQ.1 ) THEN
616 IF( M.GE.N ) THEN
617 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
618 $ LDVT, VT, LDVT, WORK, LWORK,
619 $ RWORK, DIF, IINFO )
620 ELSE
621 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
622 $ LDVT, A, LDA, WORK, LWORK,
623 $ RWORK, DIF, IINFO )
624 END IF
625 ELSE IF( IJQ.EQ.2 ) THEN
626 CALL CUNT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
627 $ LDVT, VT, LDVT, WORK, LWORK, RWORK,
628 $ DIF, IINFO )
629 END IF
630 END IF
631 RESULT( 13 ) = MAX( RESULT( 13 ), DIF )
632 *
633 * Compare S
634 *
635 DIF = ZERO
636 DIV = MAX( REAL( MNMIN )*ULP*S( 1 ),
637 $ SLAMCH( 'Safe minimum' ) )
638 DO 120 I = 1, MNMIN - 1
639 IF( SSAV( I ).LT.SSAV( I+1 ) )
640 $ DIF = ULPINV
641 IF( SSAV( I ).LT.ZERO )
642 $ DIF = ULPINV
643 DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
644 120 CONTINUE
645 RESULT( 14 ) = MAX( RESULT( 14 ), DIF )
646 130 CONTINUE
647 *
648 * End of Loop -- Check for RESULT(j) > THRESH
649 *
650 NTEST = 0
651 NFAIL = 0
652 DO 140 J = 1, 14
653 IF( RESULT( J ).GE.ZERO )
654 $ NTEST = NTEST + 1
655 IF( RESULT( J ).GE.THRESH )
656 $ NFAIL = NFAIL + 1
657 140 CONTINUE
658 *
659 IF( NFAIL.GT.0 )
660 $ NTESTF = NTESTF + 1
661 IF( NTESTF.EQ.1 ) THEN
662 WRITE( NOUNIT, FMT = 9999 )
663 WRITE( NOUNIT, FMT = 9998 )THRESH
664 NTESTF = 2
665 END IF
666 *
667 DO 150 J = 1, 14
668 IF( RESULT( J ).GE.THRESH ) THEN
669 WRITE( NOUNIT, FMT = 9997 )M, N, JTYPE, IWSPC,
670 $ IOLDSD, J, RESULT( J )
671 END IF
672 150 CONTINUE
673 *
674 NERRS = NERRS + NFAIL
675 NTESTT = NTESTT + NTEST
676 *
677 160 CONTINUE
678 *
679 170 CONTINUE
680 180 CONTINUE
681 *
682 * Summary
683 *
684 CALL ALASVM( 'CBD', NOUNIT, NERRS, NTESTT, 0 )
685 *
686 9999 FORMAT( ' SVD -- Complex Singular Value Decomposition Driver ',
687 $ / ' Matrix types (see CDRVBD for details):',
688 $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
689 $ / ' 3 = Evenly spaced singular values near 1',
690 $ / ' 4 = Evenly spaced singular values near underflow',
691 $ / ' 5 = Evenly spaced singular values near overflow',
692 $ / / ' Tests performed: ( A is dense, U and V are unitary,',
693 $ / 19X, ' S is an array, and Upartial, VTpartial, and',
694 $ / 19X, ' Spartial are partially computed U, VT and S),', / )
695 9998 FORMAT( ' Tests performed with Test Threshold = ', F8.2,
696 $ / ' CGESVD: ', /
697 $ ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
698 $ / ' 2 = | I - U**T U | / ( M ulp ) ',
699 $ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
700 $ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
701 $ ' decreasing order, else 1/ulp',
702 $ / ' 5 = | U - Upartial | / ( M ulp )',
703 $ / ' 6 = | VT - VTpartial | / ( N ulp )',
704 $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
705 $ / ' CGESDD: ', /
706 $ ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
707 $ / ' 9 = | I - U**T U | / ( M ulp ) ',
708 $ / '10 = | I - VT VT**T | / ( N ulp ) ',
709 $ / '11 = 0 if S contains min(M,N) nonnegative values in',
710 $ ' decreasing order, else 1/ulp',
711 $ / '12 = | U - Upartial | / ( M ulp )',
712 $ / '13 = | VT - VTpartial | / ( N ulp )',
713 $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )', / / )
714 9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1,
715 $ ', seed=', 4( I4, ',' ), ' test(', I1, ')=', G11.4 )
716 9996 FORMAT( ' CDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
717 $ I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
718 $ I5, ')' )
719 9995 FORMAT( ' CDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
720 $ I6, ', N=', I6, ', JTYPE=', I6, ', LSWORK=', I6, / 9X,
721 $ 'ISEED=(', 3( I5, ',' ), I5, ')' )
722 *
723 RETURN
724 *
725 * End of CDRVBD
726 *
727 END