LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
cdrvbd.f
Go to the documentation of this file.
1 *> \brief \b CDRVBD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE CDRVBD( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
12 * A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S,
13 * SSAV, E, WORK, LWORK, RWORK, IWORK, NOUNIT,
14 * INFO )
15 *
16 * .. Scalar Arguments ..
17 * INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUNIT, NSIZES,
18 * $ NTYPES
19 * REAL THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL DOTYPE( * )
23 * INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
24 * REAL E( * ), RWORK( * ), S( * ), SSAV( * )
25 * COMPLEX A( LDA, * ), ASAV( LDA, * ), U( LDU, * ),
26 * $ USAV( LDU, * ), VT( LDVT, * ),
27 * $ VTSAV( LDVT, * ), WORK( * )
28 * ..
29 *
30 *
31 *> \par Purpose:
32 * =============
33 *>
34 *> \verbatim
35 *>
36 *> CDRVBD checks the singular value decomposition (SVD) driver CGESVD
37 *> and CGESDD.
38 *> CGESVD and CGESDD factors A = U diag(S) VT, where U and VT are
39 *> unitary and diag(S) is diagonal with the entries of the array S on
40 *> its diagonal. The entries of S are the singular values, nonnegative
41 *> and stored in decreasing order. U and VT can be optionally not
42 *> computed, overwritten on A, or computed partially.
43 *>
44 *> A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN.
45 *> U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N.
46 *>
47 *> When CDRVBD is called, a number of matrix "sizes" (M's and N's)
48 *> and a number of matrix "types" are specified. For each size (M,N)
49 *> and each type of matrix, and for the minimal workspace as well as
50 *> workspace adequate to permit blocking, an M x N matrix "A" will be
51 *> generated and used to test the SVD routines. For each matrix, A will
52 *> be factored as A = U diag(S) VT and the following 12 tests computed:
53 *>
54 *> Test for CGESVD:
55 *>
56 *> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
57 *>
58 *> (2) | I - U'U | / ( M ulp )
59 *>
60 *> (3) | I - VT VT' | / ( N ulp )
61 *>
62 *> (4) S contains MNMIN nonnegative values in decreasing order.
63 *> (Return 0 if true, 1/ULP if false.)
64 *>
65 *> (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
66 *> computed U.
67 *>
68 *> (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
69 *> computed VT.
70 *>
71 *> (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
72 *> vector of singular values from the partial SVD
73 *>
74 *> Test for CGESDD:
75 *>
76 *> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
77 *>
78 *> (2) | I - U'U | / ( M ulp )
79 *>
80 *> (3) | I - VT VT' | / ( N ulp )
81 *>
82 *> (4) S contains MNMIN nonnegative values in decreasing order.
83 *> (Return 0 if true, 1/ULP if false.)
84 *>
85 *> (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
86 *> computed U.
87 *>
88 *> (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
89 *> computed VT.
90 *>
91 *> (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
92 *> vector of singular values from the partial SVD
93 *>
94 *> The "sizes" are specified by the arrays MM(1:NSIZES) and
95 *> NN(1:NSIZES); the value of each element pair (MM(j),NN(j))
96 *> specifies one size. The "types" are specified by a logical array
97 *> DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j"
98 *> will be generated.
99 *> Currently, the list of possible types is:
100 *>
101 *> (1) The zero matrix.
102 *> (2) The identity matrix.
103 *> (3) A matrix of the form U D V, where U and V are unitary and
104 *> D has evenly spaced entries 1, ..., ULP with random signs
105 *> on the diagonal.
106 *> (4) Same as (3), but multiplied by the underflow-threshold / ULP.
107 *> (5) Same as (3), but multiplied by the overflow-threshold * ULP.
108 *> \endverbatim
109 *
110 * Arguments:
111 * ==========
112 *
113 *> \param[in] NSIZES
114 *> \verbatim
115 *> NSIZES is INTEGER
116 *> The number of sizes of matrices to use. If it is zero,
117 *> CDRVBD does nothing. It must be at least zero.
118 *> \endverbatim
119 *>
120 *> \param[in] MM
121 *> \verbatim
122 *> MM is INTEGER array, dimension (NSIZES)
123 *> An array containing the matrix "heights" to be used. For
124 *> each j=1,...,NSIZES, if MM(j) is zero, then MM(j) and NN(j)
125 *> will be ignored. The MM(j) values must be at least zero.
126 *> \endverbatim
127 *>
128 *> \param[in] NN
129 *> \verbatim
130 *> NN is INTEGER array, dimension (NSIZES)
131 *> An array containing the matrix "widths" to be used. For
132 *> each j=1,...,NSIZES, if NN(j) is zero, then MM(j) and NN(j)
133 *> will be ignored. The NN(j) values must be at least zero.
134 *> \endverbatim
135 *>
136 *> \param[in] NTYPES
137 *> \verbatim
138 *> NTYPES is INTEGER
139 *> The number of elements in DOTYPE. If it is zero, CDRVBD
140 *> does nothing. It must be at least zero. If it is MAXTYP+1
141 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
142 *> defined, which is to use whatever matrices are in A and B.
143 *> This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
144 *> DOTYPE(MAXTYP+1) is .TRUE. .
145 *> \endverbatim
146 *>
147 *> \param[in] DOTYPE
148 *> \verbatim
149 *> DOTYPE is LOGICAL array, dimension (NTYPES)
150 *> If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
151 *> of type j will be generated. If NTYPES is smaller than the
152 *> maximum number of types defined (PARAMETER MAXTYP), then
153 *> types NTYPES+1 through MAXTYP will not be generated. If
154 *> NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
155 *> DOTYPE(NTYPES) will be ignored.
156 *> \endverbatim
157 *>
158 *> \param[in,out] ISEED
159 *> \verbatim
160 *> ISEED is INTEGER array, dimension (4)
161 *> On entry ISEED specifies the seed of the random number
162 *> generator. The array elements should be between 0 and 4095;
163 *> if not they will be reduced mod 4096. Also, ISEED(4) must
164 *> be odd. The random number generator uses a linear
165 *> congruential sequence limited to small integers, and so
166 *> should produce machine independent random numbers. The
167 *> values of ISEED are changed on exit, and can be used in the
168 *> next call to CDRVBD to continue the same random number
169 *> sequence.
170 *> \endverbatim
171 *>
172 *> \param[in] THRESH
173 *> \verbatim
174 *> THRESH is REAL
175 *> A test will count as "failed" if the "error", computed as
176 *> described above, exceeds THRESH. Note that the error
177 *> is scaled to be O(1), so THRESH should be a reasonably
178 *> small multiple of 1, e.g., 10 or 100. In particular,
179 *> it should not depend on the precision (single vs. double)
180 *> or the size of the matrix. It must be at least zero.
181 *> \endverbatim
182 *>
183 *> \param[out] A
184 *> \verbatim
185 *> A is COMPLEX array, dimension (LDA,max(NN))
186 *> Used to hold the matrix whose singular values are to be
187 *> computed. On exit, A contains the last matrix actually
188 *> used.
189 *> \endverbatim
190 *>
191 *> \param[in] LDA
192 *> \verbatim
193 *> LDA is INTEGER
194 *> The leading dimension of A. It must be at
195 *> least 1 and at least max( MM ).
196 *> \endverbatim
197 *>
198 *> \param[out] U
199 *> \verbatim
200 *> U is COMPLEX array, dimension (LDU,max(MM))
201 *> Used to hold the computed matrix of right singular vectors.
202 *> On exit, U contains the last such vectors actually computed.
203 *> \endverbatim
204 *>
205 *> \param[in] LDU
206 *> \verbatim
207 *> LDU is INTEGER
208 *> The leading dimension of U. It must be at
209 *> least 1 and at least max( MM ).
210 *> \endverbatim
211 *>
212 *> \param[out] VT
213 *> \verbatim
214 *> VT is COMPLEX array, dimension (LDVT,max(NN))
215 *> Used to hold the computed matrix of left singular vectors.
216 *> On exit, VT contains the last such vectors actually computed.
217 *> \endverbatim
218 *>
219 *> \param[in] LDVT
220 *> \verbatim
221 *> LDVT is INTEGER
222 *> The leading dimension of VT. It must be at
223 *> least 1 and at least max( NN ).
224 *> \endverbatim
225 *>
226 *> \param[out] ASAV
227 *> \verbatim
228 *> ASAV is COMPLEX array, dimension (LDA,max(NN))
229 *> Used to hold a different copy of the matrix whose singular
230 *> values are to be computed. On exit, A contains the last
231 *> matrix actually used.
232 *> \endverbatim
233 *>
234 *> \param[out] USAV
235 *> \verbatim
236 *> USAV is COMPLEX array, dimension (LDU,max(MM))
237 *> Used to hold a different copy of the computed matrix of
238 *> right singular vectors. On exit, USAV contains the last such
239 *> vectors actually computed.
240 *> \endverbatim
241 *>
242 *> \param[out] VTSAV
243 *> \verbatim
244 *> VTSAV is COMPLEX array, dimension (LDVT,max(NN))
245 *> Used to hold a different copy of the computed matrix of
246 *> left singular vectors. On exit, VTSAV contains the last such
247 *> vectors actually computed.
248 *> \endverbatim
249 *>
250 *> \param[out] S
251 *> \verbatim
252 *> S is REAL array, dimension (max(min(MM,NN)))
253 *> Contains the computed singular values.
254 *> \endverbatim
255 *>
256 *> \param[out] SSAV
257 *> \verbatim
258 *> SSAV is REAL array, dimension (max(min(MM,NN)))
259 *> Contains another copy of the computed singular values.
260 *> \endverbatim
261 *>
262 *> \param[out] E
263 *> \verbatim
264 *> E is REAL array, dimension (max(min(MM,NN)))
265 *> Workspace for CGESVD.
266 *> \endverbatim
267 *>
268 *> \param[out] WORK
269 *> \verbatim
270 *> WORK is COMPLEX array, dimension (LWORK)
271 *> \endverbatim
272 *>
273 *> \param[in] LWORK
274 *> \verbatim
275 *> LWORK is INTEGER
276 *> The number of entries in WORK. This must be at least
277 *> MAX(3*MIN(M,N)+MAX(M,N)**2,5*MIN(M,N),3*MAX(M,N)) for all
278 *> pairs (M,N)=(MM(j),NN(j))
279 *> \endverbatim
280 *>
281 *> \param[out] RWORK
282 *> \verbatim
283 *> RWORK is REAL array,
284 *> dimension ( 5*max(max(MM,NN)) )
285 *> \endverbatim
286 *>
287 *> \param[out] IWORK
288 *> \verbatim
289 *> IWORK is INTEGER array, dimension at least 8*min(M,N)
290 *> \endverbatim
291 *>
292 *> \param[in] NOUNIT
293 *> \verbatim
294 *> NOUNIT is INTEGER
295 *> The FORTRAN unit number for printing out error messages
296 *> (e.g., if a routine returns IINFO not equal to 0.)
297 *> \endverbatim
298 *>
299 *> \param[out] INFO
300 *> \verbatim
301 *> INFO is INTEGER
302 *> If 0, then everything ran OK.
303 *> -1: NSIZES < 0
304 *> -2: Some MM(j) < 0
305 *> -3: Some NN(j) < 0
306 *> -4: NTYPES < 0
307 *> -7: THRESH < 0
308 *> -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
309 *> -12: LDU < 1 or LDU < MMAX.
310 *> -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ).
311 *> -21: LWORK too small.
312 *> If CLATMS, or CGESVD returns an error code, the
313 *> absolute value of it is returned.
314 *> \endverbatim
315 *
316 * Authors:
317 * ========
318 *
319 *> \author Univ. of Tennessee
320 *> \author Univ. of California Berkeley
321 *> \author Univ. of Colorado Denver
322 *> \author NAG Ltd.
323 *
324 *> \date November 2011
325 *
326 *> \ingroup complex_eig
327 *
328 * =====================================================================
329  SUBROUTINE cdrvbd( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
330  $ a, lda, u, ldu, vt, ldvt, asav, usav, vtsav, s,
331  $ ssav, e, work, lwork, rwork, iwork, nounit,
332  $ info )
333 *
334 * -- LAPACK test routine (version 3.4.0) --
335 * -- LAPACK is a software package provided by Univ. of Tennessee, --
336 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
337 * November 2011
338 *
339 * .. Scalar Arguments ..
340  INTEGER info, lda, ldu, ldvt, lwork, nounit, nsizes,
341  $ ntypes
342  REAL thresh
343 * ..
344 * .. Array Arguments ..
345  LOGICAL dotype( * )
346  INTEGER iseed( 4 ), iwork( * ), mm( * ), nn( * )
347  REAL e( * ), rwork( * ), s( * ), ssav( * )
348  COMPLEX a( lda, * ), asav( lda, * ), u( ldu, * ),
349  $ usav( ldu, * ), vt( ldvt, * ),
350  $ vtsav( ldvt, * ), work( * )
351 * ..
352 *
353 * =====================================================================
354 *
355 * .. Parameters ..
356  REAL zero, one
357  parameter( zero = 0.0e+0, one = 1.0e+0 )
358  COMPLEX czero, cone
359  parameter( czero = ( 0.0e+0, 0.0e+0 ),
360  $ cone = ( 1.0e+0, 0.0e+0 ) )
361  INTEGER maxtyp
362  parameter( maxtyp = 5 )
363 * ..
364 * .. Local Scalars ..
365  LOGICAL badmm, badnn
366  CHARACTER jobq, jobu, jobvt
367  INTEGER i, iinfo, ijq, iju, ijvt, iwspc, iwtmp, j,
368  $ jsize, jtype, lswork, m, minwrk, mmax, mnmax,
369  $ mnmin, mtypes, n, nerrs, nfail, nmax, ntest,
370  $ ntestf, ntestt
371  REAL anorm, dif, div, ovfl, ulp, ulpinv, unfl
372 * ..
373 * .. Local Arrays ..
374  CHARACTER cjob( 4 )
375  INTEGER ioldsd( 4 )
376  REAL result( 14 )
377 * ..
378 * .. External Functions ..
379  REAL slamch
380  EXTERNAL slamch
381 * ..
382 * .. External Subroutines ..
383  EXTERNAL alasvm, cbdt01, cgesdd, cgesvd, clacpy, claset,
384  $ clatms, cunt01, cunt03, xerbla
385 * ..
386 * .. Intrinsic Functions ..
387  INTRINSIC abs, max, min, real
388 * ..
389 * .. Data statements ..
390  DATA cjob / 'N', 'O', 'S', 'A' /
391 * ..
392 * .. Executable Statements ..
393 *
394 * Check for errors
395 *
396  info = 0
397 *
398 * Important constants
399 *
400  nerrs = 0
401  ntestt = 0
402  ntestf = 0
403  badmm = .false.
404  badnn = .false.
405  mmax = 1
406  nmax = 1
407  mnmax = 1
408  minwrk = 1
409  DO 10 j = 1, nsizes
410  mmax = max( mmax, mm( j ) )
411  IF( mm( j ).LT.0 )
412  $ badmm = .true.
413  nmax = max( nmax, nn( j ) )
414  IF( nn( j ).LT.0 )
415  $ badnn = .true.
416  mnmax = max( mnmax, min( mm( j ), nn( j ) ) )
417  minwrk = max( minwrk, max( 3*min( mm( j ),
418  $ nn( j ) )+max( mm( j ), nn( j ) )**2, 5*min( mm( j ),
419  $ nn( j ) ), 3*max( mm( j ), nn( j ) ) ) )
420  10 continue
421 *
422 * Check for errors
423 *
424  IF( nsizes.LT.0 ) THEN
425  info = -1
426  ELSE IF( badmm ) THEN
427  info = -2
428  ELSE IF( badnn ) THEN
429  info = -3
430  ELSE IF( ntypes.LT.0 ) THEN
431  info = -4
432  ELSE IF( lda.LT.max( 1, mmax ) ) THEN
433  info = -10
434  ELSE IF( ldu.LT.max( 1, mmax ) ) THEN
435  info = -12
436  ELSE IF( ldvt.LT.max( 1, nmax ) ) THEN
437  info = -14
438  ELSE IF( minwrk.GT.lwork ) THEN
439  info = -21
440  END IF
441 *
442  IF( info.NE.0 ) THEN
443  CALL xerbla( 'CDRVBD', -info )
444  return
445  END IF
446 *
447 * Quick return if nothing to do
448 *
449  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
450  $ return
451 *
452 * More Important constants
453 *
454  unfl = slamch( 'S' )
455  ovfl = one / unfl
456  ulp = slamch( 'E' )
457  ulpinv = one / ulp
458 *
459 * Loop over sizes, types
460 *
461  nerrs = 0
462 *
463  DO 180 jsize = 1, nsizes
464  m = mm( jsize )
465  n = nn( jsize )
466  mnmin = min( m, n )
467 *
468  IF( nsizes.NE.1 ) THEN
469  mtypes = min( maxtyp, ntypes )
470  ELSE
471  mtypes = min( maxtyp+1, ntypes )
472  END IF
473 *
474  DO 170 jtype = 1, mtypes
475  IF( .NOT.dotype( jtype ) )
476  $ go to 170
477  ntest = 0
478 *
479  DO 20 j = 1, 4
480  ioldsd( j ) = iseed( j )
481  20 continue
482 *
483 * Compute "A"
484 *
485  IF( mtypes.GT.maxtyp )
486  $ go to 50
487 *
488  IF( jtype.EQ.1 ) THEN
489 *
490 * Zero matrix
491 *
492  CALL claset( 'Full', m, n, czero, czero, a, lda )
493  DO 30 i = 1, min( m, n )
494  s( i ) = zero
495  30 continue
496 *
497  ELSE IF( jtype.EQ.2 ) THEN
498 *
499 * Identity matrix
500 *
501  CALL claset( 'Full', m, n, czero, cone, a, lda )
502  DO 40 i = 1, min( m, n )
503  s( i ) = one
504  40 continue
505 *
506  ELSE
507 *
508 * (Scaled) random matrix
509 *
510  IF( jtype.EQ.3 )
511  $ anorm = one
512  IF( jtype.EQ.4 )
513  $ anorm = unfl / ulp
514  IF( jtype.EQ.5 )
515  $ anorm = ovfl*ulp
516  CALL clatms( m, n, 'U', iseed, 'N', s, 4, REAL( MNMIN ),
517  $ anorm, m-1, n-1, 'N', a, lda, work, iinfo )
518  IF( iinfo.NE.0 ) THEN
519  WRITE( nounit, fmt = 9996 )'Generator', iinfo, m, n,
520  $ jtype, ioldsd
521  info = abs( iinfo )
522  return
523  END IF
524  END IF
525 *
526  50 continue
527  CALL clacpy( 'F', m, n, a, lda, asav, lda )
528 *
529 * Do for minimal and adequate (for blocking) workspace
530 *
531  DO 160 iwspc = 1, 4
532 *
533 * Test for CGESVD
534 *
535  iwtmp = 2*min( m, n )+max( m, n )
536  lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
537  lswork = min( lswork, lwork )
538  lswork = max( lswork, 1 )
539  IF( iwspc.EQ.4 )
540  $ lswork = lwork
541 *
542  DO 60 j = 1, 14
543  result( j ) = -one
544  60 continue
545 *
546 * Factorize A
547 *
548  IF( iwspc.GT.1 )
549  $ CALL clacpy( 'F', m, n, asav, lda, a, lda )
550  CALL cgesvd( 'A', 'A', m, n, a, lda, ssav, usav, ldu,
551  $ vtsav, ldvt, work, lswork, rwork, iinfo )
552  IF( iinfo.NE.0 ) THEN
553  WRITE( nounit, fmt = 9995 )'GESVD', iinfo, m, n,
554  $ jtype, lswork, ioldsd
555  info = abs( iinfo )
556  return
557  END IF
558 *
559 * Do tests 1--4
560 *
561  CALL cbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
562  $ vtsav, ldvt, work, rwork, result( 1 ) )
563  IF( m.NE.0 .AND. n.NE.0 ) THEN
564  CALL cunt01( 'Columns', mnmin, m, usav, ldu, work,
565  $ lwork, rwork, result( 2 ) )
566  CALL cunt01( 'Rows', mnmin, n, vtsav, ldvt, work,
567  $ lwork, rwork, result( 3 ) )
568  END IF
569  result( 4 ) = 0
570  DO 70 i = 1, mnmin - 1
571  IF( ssav( i ).LT.ssav( i+1 ) )
572  $ result( 4 ) = ulpinv
573  IF( ssav( i ).LT.zero )
574  $ result( 4 ) = ulpinv
575  70 continue
576  IF( mnmin.GE.1 ) THEN
577  IF( ssav( mnmin ).LT.zero )
578  $ result( 4 ) = ulpinv
579  END IF
580 *
581 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
582 *
583  result( 5 ) = zero
584  result( 6 ) = zero
585  result( 7 ) = zero
586  DO 100 iju = 0, 3
587  DO 90 ijvt = 0, 3
588  IF( ( iju.EQ.3 .AND. ijvt.EQ.3 ) .OR.
589  $ ( iju.EQ.1 .AND. ijvt.EQ.1 ) )go to 90
590  jobu = cjob( iju+1 )
591  jobvt = cjob( ijvt+1 )
592  CALL clacpy( 'F', m, n, asav, lda, a, lda )
593  CALL cgesvd( jobu, jobvt, m, n, a, lda, s, u, ldu,
594  $ vt, ldvt, work, lswork, rwork, iinfo )
595 *
596 * Compare U
597 *
598  dif = zero
599  IF( m.GT.0 .AND. n.GT.0 ) THEN
600  IF( iju.EQ.1 ) THEN
601  CALL cunt03( 'C', m, mnmin, m, mnmin, usav,
602  $ ldu, a, lda, work, lwork, rwork,
603  $ dif, iinfo )
604  ELSE IF( iju.EQ.2 ) THEN
605  CALL cunt03( 'C', m, mnmin, m, mnmin, usav,
606  $ ldu, u, ldu, work, lwork, rwork,
607  $ dif, iinfo )
608  ELSE IF( iju.EQ.3 ) THEN
609  CALL cunt03( 'C', m, m, m, mnmin, usav, ldu,
610  $ u, ldu, work, lwork, rwork, dif,
611  $ iinfo )
612  END IF
613  END IF
614  result( 5 ) = max( result( 5 ), dif )
615 *
616 * Compare VT
617 *
618  dif = zero
619  IF( m.GT.0 .AND. n.GT.0 ) THEN
620  IF( ijvt.EQ.1 ) THEN
621  CALL cunt03( 'R', n, mnmin, n, mnmin, vtsav,
622  $ ldvt, a, lda, work, lwork,
623  $ rwork, dif, iinfo )
624  ELSE IF( ijvt.EQ.2 ) THEN
625  CALL cunt03( 'R', n, mnmin, n, mnmin, vtsav,
626  $ ldvt, vt, ldvt, work, lwork,
627  $ rwork, dif, iinfo )
628  ELSE IF( ijvt.EQ.3 ) THEN
629  CALL cunt03( 'R', n, n, n, mnmin, vtsav,
630  $ ldvt, vt, ldvt, work, lwork,
631  $ rwork, dif, iinfo )
632  END IF
633  END IF
634  result( 6 ) = max( result( 6 ), dif )
635 *
636 * Compare S
637 *
638  dif = zero
639  div = max( REAL( mnmin )*ulp*s( 1 ),
640  $ slamch( 'Safe minimum' ) )
641  DO 80 i = 1, mnmin - 1
642  IF( ssav( i ).LT.ssav( i+1 ) )
643  $ dif = ulpinv
644  IF( ssav( i ).LT.zero )
645  $ dif = ulpinv
646  dif = max( dif, abs( ssav( i )-s( i ) ) / div )
647  80 continue
648  result( 7 ) = max( result( 7 ), dif )
649  90 continue
650  100 continue
651 *
652 * Test for CGESDD
653 *
654  iwtmp = 2*mnmin*mnmin + 2*mnmin + max( m, n )
655  lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
656  lswork = min( lswork, lwork )
657  lswork = max( lswork, 1 )
658  IF( iwspc.EQ.4 )
659  $ lswork = lwork
660 *
661 * Factorize A
662 *
663  CALL clacpy( 'F', m, n, asav, lda, a, lda )
664  CALL cgesdd( 'A', m, n, a, lda, ssav, usav, ldu, vtsav,
665  $ ldvt, work, lswork, rwork, iwork, iinfo )
666  IF( iinfo.NE.0 ) THEN
667  WRITE( nounit, fmt = 9995 )'GESDD', iinfo, m, n,
668  $ jtype, lswork, ioldsd
669  info = abs( iinfo )
670  return
671  END IF
672 *
673 * Do tests 1--4
674 *
675  CALL cbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
676  $ vtsav, ldvt, work, rwork, result( 8 ) )
677  IF( m.NE.0 .AND. n.NE.0 ) THEN
678  CALL cunt01( 'Columns', mnmin, m, usav, ldu, work,
679  $ lwork, rwork, result( 9 ) )
680  CALL cunt01( 'Rows', mnmin, n, vtsav, ldvt, work,
681  $ lwork, rwork, result( 10 ) )
682  END IF
683  result( 11 ) = 0
684  DO 110 i = 1, mnmin - 1
685  IF( ssav( i ).LT.ssav( i+1 ) )
686  $ result( 11 ) = ulpinv
687  IF( ssav( i ).LT.zero )
688  $ result( 11 ) = ulpinv
689  110 continue
690  IF( mnmin.GE.1 ) THEN
691  IF( ssav( mnmin ).LT.zero )
692  $ result( 11 ) = ulpinv
693  END IF
694 *
695 * Do partial SVDs, comparing to SSAV, USAV, and VTSAV
696 *
697  result( 12 ) = zero
698  result( 13 ) = zero
699  result( 14 ) = zero
700  DO 130 ijq = 0, 2
701  jobq = cjob( ijq+1 )
702  CALL clacpy( 'F', m, n, asav, lda, a, lda )
703  CALL cgesdd( jobq, m, n, a, lda, s, u, ldu, vt, ldvt,
704  $ work, lswork, rwork, iwork, iinfo )
705 *
706 * Compare U
707 *
708  dif = zero
709  IF( m.GT.0 .AND. n.GT.0 ) THEN
710  IF( ijq.EQ.1 ) THEN
711  IF( m.GE.n ) THEN
712  CALL cunt03( 'C', m, mnmin, m, mnmin, usav,
713  $ ldu, a, lda, work, lwork, rwork,
714  $ dif, iinfo )
715  ELSE
716  CALL cunt03( 'C', m, mnmin, m, mnmin, usav,
717  $ ldu, u, ldu, work, lwork, rwork,
718  $ dif, iinfo )
719  END IF
720  ELSE IF( ijq.EQ.2 ) THEN
721  CALL cunt03( 'C', m, mnmin, m, mnmin, usav, ldu,
722  $ u, ldu, work, lwork, rwork, dif,
723  $ iinfo )
724  END IF
725  END IF
726  result( 12 ) = max( result( 12 ), dif )
727 *
728 * Compare VT
729 *
730  dif = zero
731  IF( m.GT.0 .AND. n.GT.0 ) THEN
732  IF( ijq.EQ.1 ) THEN
733  IF( m.GE.n ) THEN
734  CALL cunt03( 'R', n, mnmin, n, mnmin, vtsav,
735  $ ldvt, vt, ldvt, work, lwork,
736  $ rwork, dif, iinfo )
737  ELSE
738  CALL cunt03( 'R', n, mnmin, n, mnmin, vtsav,
739  $ ldvt, a, lda, work, lwork,
740  $ rwork, dif, iinfo )
741  END IF
742  ELSE IF( ijq.EQ.2 ) THEN
743  CALL cunt03( 'R', n, mnmin, n, mnmin, vtsav,
744  $ ldvt, vt, ldvt, work, lwork, rwork,
745  $ dif, iinfo )
746  END IF
747  END IF
748  result( 13 ) = max( result( 13 ), dif )
749 *
750 * Compare S
751 *
752  dif = zero
753  div = max( REAL( mnmin )*ulp*s( 1 ),
754  $ slamch( 'Safe minimum' ) )
755  DO 120 i = 1, mnmin - 1
756  IF( ssav( i ).LT.ssav( i+1 ) )
757  $ dif = ulpinv
758  IF( ssav( i ).LT.zero )
759  $ dif = ulpinv
760  dif = max( dif, abs( ssav( i )-s( i ) ) / div )
761  120 continue
762  result( 14 ) = max( result( 14 ), dif )
763  130 continue
764 *
765 * End of Loop -- Check for RESULT(j) > THRESH
766 *
767  ntest = 0
768  nfail = 0
769  DO 140 j = 1, 14
770  IF( result( j ).GE.zero )
771  $ ntest = ntest + 1
772  IF( result( j ).GE.thresh )
773  $ nfail = nfail + 1
774  140 continue
775 *
776  IF( nfail.GT.0 )
777  $ ntestf = ntestf + 1
778  IF( ntestf.EQ.1 ) THEN
779  WRITE( nounit, fmt = 9999 )
780  WRITE( nounit, fmt = 9998 )thresh
781  ntestf = 2
782  END IF
783 *
784  DO 150 j = 1, 14
785  IF( result( j ).GE.thresh ) THEN
786  WRITE( nounit, fmt = 9997 )m, n, jtype, iwspc,
787  $ ioldsd, j, result( j )
788  END IF
789  150 continue
790 *
791  nerrs = nerrs + nfail
792  ntestt = ntestt + ntest
793 *
794  160 continue
795 *
796  170 continue
797  180 continue
798 *
799 * Summary
800 *
801  CALL alasvm( 'CBD', nounit, nerrs, ntestt, 0 )
802 *
803  9999 format( ' SVD -- Complex Singular Value Decomposition Driver ',
804  $ / ' Matrix types (see CDRVBD for details):',
805  $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
806  $ / ' 3 = Evenly spaced singular values near 1',
807  $ / ' 4 = Evenly spaced singular values near underflow',
808  $ / ' 5 = Evenly spaced singular values near overflow',
809  $ / / ' Tests performed: ( A is dense, U and V are unitary,',
810  $ / 19x, ' S is an array, and Upartial, VTpartial, and',
811  $ / 19x, ' Spartial are partially computed U, VT and S),', / )
812  9998 format( ' Tests performed with Test Threshold = ', f8.2,
813  $ / ' CGESVD: ', /
814  $ ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
815  $ / ' 2 = | I - U**T U | / ( M ulp ) ',
816  $ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
817  $ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
818  $ ' decreasing order, else 1/ulp',
819  $ / ' 5 = | U - Upartial | / ( M ulp )',
820  $ / ' 6 = | VT - VTpartial | / ( N ulp )',
821  $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
822  $ / ' CGESDD: ', /
823  $ ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
824  $ / ' 9 = | I - U**T U | / ( M ulp ) ',
825  $ / '10 = | I - VT VT**T | / ( N ulp ) ',
826  $ / '11 = 0 if S contains min(M,N) nonnegative values in',
827  $ ' decreasing order, else 1/ulp',
828  $ / '12 = | U - Upartial | / ( M ulp )',
829  $ / '13 = | VT - VTpartial | / ( N ulp )',
830  $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )', / / )
831  9997 format( ' M=', i5, ', N=', i5, ', type ', i1, ', IWS=', i1,
832  $ ', seed=', 4( i4, ',' ), ' test(', i1, ')=', g11.4 )
833  9996 format( ' CDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
834  $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
835  $ i5, ')' )
836  9995 format( ' CDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
837  $ i6, ', N=', i6, ', JTYPE=', i6, ', LSWORK=', i6, / 9x,
838  $ 'ISEED=(', 3( i5, ',' ), i5, ')' )
839 *
840  return
841 *
842 * End of CDRVBD
843 *
844  END