OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
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*> CGESDD, CGESVJ, CGEJSV, CGESVDX, and CGESVDQ.
38*>
39*> CGESVD and CGESDD factors A = U diag(S) VT, where U and VT are
40*> unitary and diag(S) is diagonal with the entries of the array S on
41*> its diagonal. The entries of S are the singular values, nonnegative
42*> and stored in decreasing order. U and VT can be optionally not
43*> computed, overwritten on A, or computed partially.
44*>
45*> A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN.
46*> U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N.
47*>
48*> When CDRVBD is called, a number of matrix "sizes" (M's and N's)
49*> and a number of matrix "types" are specified. For each size (M,N)
50*> and each type of matrix, and for the minimal workspace as well as
51*> workspace adequate to permit blocking, an M x N matrix "A" will be
52*> generated and used to test the SVD routines. For each matrix, A will
53*> be factored as A = U diag(S) VT and the following 12 tests computed:
54*>
55*> Test for CGESVD:
56*>
57*> (1) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
58*>
59*> (2) | I - U'U | / ( M ulp )
60*>
61*> (3) | I - VT VT' | / ( N ulp )
62*>
63*> (4) S contains MNMIN nonnegative values in decreasing order.
64*> (Return 0 if true, 1/ULP if false.)
65*>
66*> (5) | U - Upartial | / ( M ulp ) where Upartial is a partially
67*> computed U.
68*>
69*> (6) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
70*> computed VT.
71*>
72*> (7) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
73*> vector of singular values from the partial SVD
74*>
75*> Test for CGESDD:
76*>
77*> (8) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
78*>
79*> (9) | I - U'U | / ( M ulp )
80*>
81*> (10) | I - VT VT' | / ( N ulp )
82*>
83*> (11) S contains MNMIN nonnegative values in decreasing order.
84*> (Return 0 if true, 1/ULP if false.)
85*>
86*> (12) | U - Upartial | / ( M ulp ) where Upartial is a partially
87*> computed U.
88*>
89*> (13) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
90*> computed VT.
91*>
92*> (14) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
93*> vector of singular values from the partial SVD
94*>
95*> Test for CGESVDQ:
96*>
97*> (36) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
98*>
99*> (37) | I - U'U | / ( M ulp )
100*>
101*> (38) | I - VT VT' | / ( N ulp )
102*>
103*> (39) S contains MNMIN nonnegative values in decreasing order.
104*> (Return 0 if true, 1/ULP if false.)
105*>
106*> Test for CGESVJ:
107*>
108*> (15) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
109*>
110*> (16) | I - U'U | / ( M ulp )
111*>
112*> (17) | I - VT VT' | / ( N ulp )
113*>
114*> (18) S contains MNMIN nonnegative values in decreasing order.
115*> (Return 0 if true, 1/ULP if false.)
116*>
117*> Test for CGEJSV:
118*>
119*> (19) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
120*>
121*> (20) | I - U'U | / ( M ulp )
122*>
123*> (21) | I - VT VT' | / ( N ulp )
124*>
125*> (22) S contains MNMIN nonnegative values in decreasing order.
126*> (Return 0 if true, 1/ULP if false.)
127*>
128*> Test for CGESVDX( 'V', 'V', 'A' )/CGESVDX( 'N', 'N', 'A' )
129*>
130*> (23) | A - U diag(S) VT | / ( |A| max(M,N) ulp )
131*>
132*> (24) | I - U'U | / ( M ulp )
133*>
134*> (25) | I - VT VT' | / ( N ulp )
135*>
136*> (26) S contains MNMIN nonnegative values in decreasing order.
137*> (Return 0 if true, 1/ULP if false.)
138*>
139*> (27) | U - Upartial | / ( M ulp ) where Upartial is a partially
140*> computed U.
141*>
142*> (28) | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
143*> computed VT.
144*>
145*> (29) | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
146*> vector of singular values from the partial SVD
147*>
148*> Test for CGESVDX( 'V', 'V', 'I' )
149*>
150*> (30) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
151*>
152*> (31) | I - U'U | / ( M ulp )
153*>
154*> (32) | I - VT VT' | / ( N ulp )
155*>
156*> Test for CGESVDX( 'V', 'V', 'V' )
157*>
158*> (33) | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
159*>
160*> (34) | I - U'U | / ( M ulp )
161*>
162*> (35) | I - VT VT' | / ( N ulp )
163*>
164*> The "sizes" are specified by the arrays MM(1:NSIZES) and
165*> NN(1:NSIZES); the value of each element pair (MM(j),NN(j))
166*> specifies one size. The "types" are specified by a logical array
167*> DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j"
168*> will be generated.
169*> Currently, the list of possible types is:
170*>
171*> (1) The zero matrix.
172*> (2) The identity matrix.
173*> (3) A matrix of the form U D V, where U and V are unitary and
174*> D has evenly spaced entries 1, ..., ULP with random signs
175*> on the diagonal.
176*> (4) Same as (3), but multiplied by the underflow-threshold / ULP.
177*> (5) Same as (3), but multiplied by the overflow-threshold * ULP.
178*> \endverbatim
179*
180* Arguments:
181* ==========
182*
183*> \param[in] NSIZES
184*> \verbatim
185*> NSIZES is INTEGER
186*> The number of sizes of matrices to use. If it is zero,
187*> CDRVBD does nothing. It must be at least zero.
188*> \endverbatim
189*>
190*> \param[in] MM
191*> \verbatim
192*> MM is INTEGER array, dimension (NSIZES)
193*> An array containing the matrix "heights" to be used. For
194*> each j=1,...,NSIZES, if MM(j) is zero, then MM(j) and NN(j)
195*> will be ignored. The MM(j) values must be at least zero.
196*> \endverbatim
197*>
198*> \param[in] NN
199*> \verbatim
200*> NN is INTEGER array, dimension (NSIZES)
201*> An array containing the matrix "widths" to be used. For
202*> each j=1,...,NSIZES, if NN(j) is zero, then MM(j) and NN(j)
203*> will be ignored. The NN(j) values must be at least zero.
204*> \endverbatim
205*>
206*> \param[in] NTYPES
207*> \verbatim
208*> NTYPES is INTEGER
209*> The number of elements in DOTYPE. If it is zero, CDRVBD
210*> does nothing. It must be at least zero. If it is MAXTYP+1
211*> and NSIZES is 1, then an additional type, MAXTYP+1 is
212*> defined, which is to use whatever matrices are in A and B.
213*> This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
214*> DOTYPE(MAXTYP+1) is .TRUE. .
215*> \endverbatim
216*>
217*> \param[in] DOTYPE
218*> \verbatim
219*> DOTYPE is LOGICAL array, dimension (NTYPES)
220*> If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
221*> of type j will be generated. If NTYPES is smaller than the
222*> maximum number of types defined (PARAMETER MAXTYP), then
223*> types NTYPES+1 through MAXTYP will not be generated. If
224*> NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
225*> DOTYPE(NTYPES) will be ignored.
226*> \endverbatim
227*>
228*> \param[in,out] ISEED
229*> \verbatim
230*> ISEED is INTEGER array, dimension (4)
231*> On entry ISEED specifies the seed of the random number
232*> generator. The array elements should be between 0 and 4095;
233*> if not they will be reduced mod 4096. Also, ISEED(4) must
234*> be odd. The random number generator uses a linear
235*> congruential sequence limited to small integers, and so
236*> should produce machine independent random numbers. The
237*> values of ISEED are changed on exit, and can be used in the
238*> next call to CDRVBD to continue the same random number
239*> sequence.
240*> \endverbatim
241*>
242*> \param[in] THRESH
243*> \verbatim
244*> THRESH is REAL
245*> A test will count as "failed" if the "error", computed as
246*> described above, exceeds THRESH. Note that the error
247*> is scaled to be O(1), so THRESH should be a reasonably
248*> small multiple of 1, e.g., 10 or 100. In particular,
249*> it should not depend on the precision (single vs. double)
250*> or the size of the matrix. It must be at least zero.
251*> \endverbatim
252*>
253*> \param[out] A
254*> \verbatim
255*> A is COMPLEX array, dimension (LDA,max(NN))
256*> Used to hold the matrix whose singular values are to be
257*> computed. On exit, A contains the last matrix actually
258*> used.
259*> \endverbatim
260*>
261*> \param[in] LDA
262*> \verbatim
263*> LDA is INTEGER
264*> The leading dimension of A. It must be at
265*> least 1 and at least max( MM ).
266*> \endverbatim
267*>
268*> \param[out] U
269*> \verbatim
270*> U is COMPLEX array, dimension (LDU,max(MM))
271*> Used to hold the computed matrix of right singular vectors.
272*> On exit, U contains the last such vectors actually computed.
273*> \endverbatim
274*>
275*> \param[in] LDU
276*> \verbatim
277*> LDU is INTEGER
278*> The leading dimension of U. It must be at
279*> least 1 and at least max( MM ).
280*> \endverbatim
281*>
282*> \param[out] VT
283*> \verbatim
284*> VT is COMPLEX array, dimension (LDVT,max(NN))
285*> Used to hold the computed matrix of left singular vectors.
286*> On exit, VT contains the last such vectors actually computed.
287*> \endverbatim
288*>
289*> \param[in] LDVT
290*> \verbatim
291*> LDVT is INTEGER
292*> The leading dimension of VT. It must be at
293*> least 1 and at least max( NN ).
294*> \endverbatim
295*>
296*> \param[out] ASAV
297*> \verbatim
298*> ASAV is COMPLEX array, dimension (LDA,max(NN))
299*> Used to hold a different copy of the matrix whose singular
300*> values are to be computed. On exit, A contains the last
301*> matrix actually used.
302*> \endverbatim
303*>
304*> \param[out] USAV
305*> \verbatim
306*> USAV is COMPLEX array, dimension (LDU,max(MM))
307*> Used to hold a different copy of the computed matrix of
308*> right singular vectors. On exit, USAV contains the last such
309*> vectors actually computed.
310*> \endverbatim
311*>
312*> \param[out] VTSAV
313*> \verbatim
314*> VTSAV is COMPLEX array, dimension (LDVT,max(NN))
315*> Used to hold a different copy of the computed matrix of
316*> left singular vectors. On exit, VTSAV contains the last such
317*> vectors actually computed.
318*> \endverbatim
319*>
320*> \param[out] S
321*> \verbatim
322*> S is REAL array, dimension (max(min(MM,NN)))
323*> Contains the computed singular values.
324*> \endverbatim
325*>
326*> \param[out] SSAV
327*> \verbatim
328*> SSAV is REAL array, dimension (max(min(MM,NN)))
329*> Contains another copy of the computed singular values.
330*> \endverbatim
331*>
332*> \param[out] E
333*> \verbatim
334*> E is REAL array, dimension (max(min(MM,NN)))
335*> Workspace for CGESVD.
336*> \endverbatim
337*>
338*> \param[out] WORK
339*> \verbatim
340*> WORK is COMPLEX array, dimension (LWORK)
341*> \endverbatim
342*>
343*> \param[in] LWORK
344*> \verbatim
345*> LWORK is INTEGER
346*> The number of entries in WORK. This must be at least
347*> MAX(3*MIN(M,N)+MAX(M,N)**2,5*MIN(M,N),3*MAX(M,N)) for all
348*> pairs (M,N)=(MM(j),NN(j))
349*> \endverbatim
350*>
351*> \param[out] RWORK
352*> \verbatim
353*> RWORK is REAL array,
354*> dimension ( 5*max(max(MM,NN)) )
355*> \endverbatim
356*>
357*> \param[out] IWORK
358*> \verbatim
359*> IWORK is INTEGER array, dimension at least 8*min(M,N)
360*> \endverbatim
361*>
362*> \param[in] NOUNIT
363*> \verbatim
364*> NOUNIT is INTEGER
365*> The FORTRAN unit number for printing out error messages
366*> (e.g., if a routine returns IINFO not equal to 0.)
367*> \endverbatim
368*>
369*> \param[out] INFO
370*> \verbatim
371*> INFO is INTEGER
372*> If 0, then everything ran OK.
373*> -1: NSIZES < 0
374*> -2: Some MM(j) < 0
375*> -3: Some NN(j) < 0
376*> -4: NTYPES < 0
377*> -7: THRESH < 0
378*> -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
379*> -12: LDU < 1 or LDU < MMAX.
380*> -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ).
381*> -29: LWORK too small.
382*> If CLATMS, or CGESVD returns an error code, the
383*> absolute value of it is returned.
384*> \endverbatim
385*
386* Authors:
387* ========
388*
389*> \author Univ. of Tennessee
390*> \author Univ. of California Berkeley
391*> \author Univ. of Colorado Denver
392*> \author NAG Ltd.
393*
394*> \ingroup complex_eig
395*
396* =====================================================================
397 SUBROUTINE cdrvbd( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
398 $ A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S,
399 $ SSAV, E, WORK, LWORK, RWORK, IWORK, NOUNIT,
400 $ INFO )
401*
402* -- LAPACK test routine --
403* -- LAPACK is a software package provided by Univ. of Tennessee, --
404* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
405*
406 IMPLICIT NONE
407*
408* .. Scalar Arguments ..
409 INTEGER INFO, LDA, LDU, LDVT, LWORK, NOUNIT, NSIZES,
410 $ NTYPES
411 REAL THRESH
412* ..
413* .. Array Arguments ..
414 LOGICAL DOTYPE( * )
415 INTEGER ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
416 REAL E( * ), RWORK( * ), S( * ), SSAV( * )
417 COMPLEX A( LDA, * ), ASAV( LDA, * ), U( LDU, * ),
418 $ usav( ldu, * ), vt( ldvt, * ),
419 $ vtsav( ldvt, * ), work( * )
420* ..
421*
422* =====================================================================
423*
424* .. Parameters ..
425 REAL ZERO, ONE, TWO, HALF
426 PARAMETER ( ZERO = 0.0e0, one = 1.0e0, two = 2.0e0,
427 $ half = 0.5e0 )
428 COMPLEX CZERO, CONE
429 parameter( czero = ( 0.0e+0, 0.0e+0 ),
430 $ cone = ( 1.0e+0, 0.0e+0 ) )
431 INTEGER MAXTYP
432 parameter( maxtyp = 5 )
433* ..
434* .. Local Scalars ..
435 LOGICAL BADMM, BADNN
436 CHARACTER JOBQ, JOBU, JOBVT, RANGE
437 INTEGER I, IINFO, IJQ, IJU, IJVT, IL, IU, ITEMP,
438 $ iwspc, iwtmp, j, jsize, jtype, lswork, m,
439 $ minwrk, mmax, mnmax, mnmin, mtypes, n,
440 $ nerrs, nfail, nmax, ns, nsi, nsv, ntest,
441 $ ntestf, ntestt, lrwork
442 REAL ANORM, DIF, DIV, OVFL, RTUNFL, ULP, ULPINV,
443 $ UNFL, VL, VU
444* ..
445* .. Local Scalars for CGESVDQ ..
446 INTEGER LIWORK, NUMRANK
447* ..
448* .. Local Arrays ..
449 CHARACTER CJOB( 4 ), CJOBR( 3 ), CJOBV( 2 )
450 INTEGER IOLDSD( 4 ), ISEED2( 4 )
451 REAL RESULT( 39 )
452* ..
453* .. External Functions ..
454 REAL SLAMCH, SLARND
455 EXTERNAL SLAMCH, SLARND
456* ..
457* .. External Subroutines ..
458 EXTERNAL alasvm, xerbla, cbdt01, cbdt05, cgesdd,
459 $ cgesvd, cgesvdq, cgesvj, cgejsv, cgesvdx,
460 $ clacpy, claset, clatms, cunt01, cunt03
461* ..
462* .. Intrinsic Functions ..
463 INTRINSIC abs, real, max, min
464* ..
465* .. Scalars in Common ..
466 CHARACTER*32 SRNAMT
467* ..
468* .. Common blocks ..
469 COMMON / srnamc / srnamt
470* ..
471* .. Data statements ..
472 DATA cjob / 'N', 'O', 'S', 'A' /
473 DATA cjobr / 'A', 'V', 'I' /
474 DATA cjobv / 'N', 'V' /
475* ..
476* .. Executable Statements ..
477*
478* Check for errors
479*
480 info = 0
481*
482* Important constants
483*
484 nerrs = 0
485 ntestt = 0
486 ntestf = 0
487 badmm = .false.
488 badnn = .false.
489 mmax = 1
490 nmax = 1
491 mnmax = 1
492 minwrk = 1
493 DO 10 j = 1, nsizes
494 mmax = max( mmax, mm( j ) )
495 IF( mm( j ).LT.0 )
496 $ badmm = .true.
497 nmax = max( nmax, nn( j ) )
498 IF( nn( j ).LT.0 )
499 $ badnn = .true.
500 mnmax = max( mnmax, min( mm( j ), nn( j ) ) )
501 minwrk = max( minwrk, max( 3*min( mm( j ),
502 $ nn( j ) )+max( mm( j ), nn( j ) )**2, 5*min( mm( j ),
503 $ nn( j ) ), 3*max( mm( j ), nn( j ) ) ) )
504 10 CONTINUE
505*
506* Check for errors
507*
508 IF( nsizes.LT.0 ) THEN
509 info = -1
510 ELSE IF( badmm ) THEN
511 info = -2
512 ELSE IF( badnn ) THEN
513 info = -3
514 ELSE IF( ntypes.LT.0 ) THEN
515 info = -4
516 ELSE IF( lda.LT.max( 1, mmax ) ) THEN
517 info = -10
518 ELSE IF( ldu.LT.max( 1, mmax ) ) THEN
519 info = -12
520 ELSE IF( ldvt.LT.max( 1, nmax ) ) THEN
521 info = -14
522 ELSE IF( minwrk.GT.lwork ) THEN
523 info = -21
524 END IF
525*
526 IF( info.NE.0 ) THEN
527 CALL xerbla( 'CDRVBD', -info )
528 RETURN
529 END IF
530*
531* Quick return if nothing to do
532*
533 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
534 $ RETURN
535*
536* More Important constants
537*
538 unfl = slamch( 'S' )
539 ovfl = one / unfl
540 ulp = slamch( 'E' )
541 ulpinv = one / ulp
542 rtunfl = sqrt( unfl )
543*
544* Loop over sizes, types
545*
546 nerrs = 0
547*
548 DO 310 jsize = 1, nsizes
549 m = mm( jsize )
550 n = nn( jsize )
551 mnmin = min( m, n )
552*
553 IF( nsizes.NE.1 ) THEN
554 mtypes = min( maxtyp, ntypes )
555 ELSE
556 mtypes = min( maxtyp+1, ntypes )
557 END IF
558*
559 DO 300 jtype = 1, mtypes
560 IF( .NOT.dotype( jtype ) )
561 $ GO TO 300
562 ntest = 0
563*
564 DO 20 j = 1, 4
565 ioldsd( j ) = iseed( j )
566 20 CONTINUE
567*
568* Compute "A"
569*
570 IF( mtypes.GT.maxtyp )
571 $ GO TO 50
572*
573 IF( jtype.EQ.1 ) THEN
574*
575* Zero matrix
576*
577 CALL claset( 'Full', m, n, czero, czero, a, lda )
578 DO 30 i = 1, min( m, n )
579 s( i ) = zero
580 30 CONTINUE
581*
582 ELSE IF( jtype.EQ.2 ) THEN
583*
584* Identity matrix
585*
586 CALL claset( 'Full', m, n, czero, cone, a, lda )
587 DO 40 i = 1, min( m, n )
588 s( i ) = one
589 40 CONTINUE
590*
591 ELSE
592*
593* (Scaled) random matrix
594*
595 IF( jtype.EQ.3 )
596 $ anorm = one
597 IF( jtype.EQ.4 )
598 $ anorm = unfl / ulp
599 IF( jtype.EQ.5 )
600 $ anorm = ovfl*ulp
601 CALL clatms( m, n, 'U', iseed, 'N', s, 4, real( mnmin ),
602 $ anorm, m-1, n-1, 'n', A, LDA, WORK, IINFO )
603.NE. IF( IINFO0 ) THEN
604 WRITE( NOUNIT, FMT = 9996 )'generator', IINFO, M, N,
605 $ JTYPE, IOLDSD
606 INFO = ABS( IINFO )
607 RETURN
608 END IF
609 END IF
610*
611 50 CONTINUE
612 CALL CLACPY( 'f', M, N, A, LDA, ASAV, LDA )
613*
614* Do for minimal and adequate (for blocking) workspace
615*
616 DO 290 IWSPC = 1, 4
617*
618* Test for CGESVD
619*
620 IWTMP = 2*MIN( M, N )+MAX( M, N )
621 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
622 LSWORK = MIN( LSWORK, LWORK )
623 LSWORK = MAX( LSWORK, 1 )
624.EQ. IF( IWSPC4 )
625 $ LSWORK = LWORK
626*
627 DO 60 J = 1, 35
628 RESULT( J ) = -ONE
629 60 CONTINUE
630*
631* Factorize A
632*
633.GT. IF( IWSPC1 )
634 $ CALL CLACPY( 'f', M, N, ASAV, LDA, A, LDA )
635 SRNAMT = 'cgesvd'
636 CALL CGESVD( 'a', 'a', M, N, A, LDA, SSAV, USAV, LDU,
637 $ VTSAV, LDVT, WORK, LSWORK, RWORK, IINFO )
638.NE. IF( IINFO0 ) THEN
639 WRITE( NOUNIT, FMT = 9995 )'gesvd', IINFO, M, N,
640 $ JTYPE, LSWORK, IOLDSD
641 INFO = ABS( IINFO )
642 RETURN
643 END IF
644*
645* Do tests 1--4
646*
647 CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
648 $ VTSAV, LDVT, WORK, RWORK, RESULT( 1 ) )
649.NE..AND..NE. IF( M0 N0 ) THEN
650 CALL CUNT01( 'columns', MNMIN, M, USAV, LDU, WORK,
651 $ LWORK, RWORK, RESULT( 2 ) )
652 CALL CUNT01( 'rows', MNMIN, N, VTSAV, LDVT, WORK,
653 $ LWORK, RWORK, RESULT( 3 ) )
654 END IF
655 RESULT( 4 ) = 0
656 DO 70 I = 1, MNMIN - 1
657.LT. IF( SSAV( I )SSAV( I+1 ) )
658 $ RESULT( 4 ) = ULPINV
659.LT. IF( SSAV( I )ZERO )
660 $ RESULT( 4 ) = ULPINV
661 70 CONTINUE
662.GE. IF( MNMIN1 ) THEN
663.LT. IF( SSAV( MNMIN )ZERO )
664 $ RESULT( 4 ) = ULPINV
665 END IF
666*
667* Do partial SVDs, comparing to SSAV, USAV, and VTSAV
668*
669 RESULT( 5 ) = ZERO
670 RESULT( 6 ) = ZERO
671 RESULT( 7 ) = ZERO
672 DO 100 IJU = 0, 3
673 DO 90 IJVT = 0, 3
674.EQ..AND..EQ..OR. IF( ( IJU3 IJVT3 )
675.EQ..AND..EQ. $ ( IJU1 IJVT1 ) )GO TO 90
676 JOBU = CJOB( IJU+1 )
677 JOBVT = CJOB( IJVT+1 )
678 CALL CLACPY( 'f', m, n, asav, lda, a, lda )
679 srnamt = 'CGESVD'
680 CALL cgesvd( jobu, jobvt, m, n, a, lda, s, u, ldu,
681 $ vt, ldvt, work, lswork, rwork, iinfo )
682*
683* Compare U
684*
685 dif = zero
686 IF( m.GT.0 .AND. n.GT.0 ) THEN
687 IF( iju.EQ.1 ) THEN
688 CALL cunt03( 'C', m, mnmin, m, mnmin, usav,
689 $ ldu, a, lda, work, lwork, rwork,
690 $ dif, iinfo )
691 ELSE IF( iju.EQ.2 ) THEN
692 CALL cunt03( 'C', m, mnmin, m, mnmin, usav,
693 $ ldu, u, ldu, work, lwork, rwork,
694 $ dif, iinfo )
695 ELSE IF( iju.EQ.3 ) THEN
696 CALL cunt03( 'C', m, m, m, mnmin, usav, ldu,
697 $ u, ldu, work, lwork, rwork, dif,
698 $ iinfo )
699 END IF
700 END IF
701 result( 5 ) = max( result( 5 ), dif )
702*
703* Compare VT
704*
705 dif = zero
706 IF( m.GT.0 .AND. n.GT.0 ) THEN
707 IF( ijvt.EQ.1 ) THEN
708 CALL cunt03( 'R', n, mnmin, n, mnmin, vtsav,
709 $ ldvt, a, lda, work, lwork,
710 $ rwork, dif, iinfo )
711 ELSE IF( ijvt.EQ.2 ) THEN
712 CALL cunt03( 'R', n, mnmin, n, mnmin, vtsav,
713 $ ldvt, vt, ldvt, work, lwork,
714 $ rwork, dif, iinfo )
715 ELSE IF( ijvt.EQ.3 ) THEN
716 CALL cunt03( 'R', n, n, n, mnmin, vtsav,
717 $ ldvt, vt, ldvt, work, lwork,
718 $ rwork, dif, iinfo )
719 END IF
720 END IF
721 result( 6 ) = max( result( 6 ), dif )
722*
723* Compare S
724*
725 dif = zero
726 div = max( real( mnmin )*ulp*s( 1 ),
727 $ slamch( 'Safe minimum' ) )
728 DO 80 i = 1, mnmin - 1
729 IF( ssav( i ).LT.ssav( i+1 ) )
730 $ dif = ulpinv
731 IF( ssav( i ).LT.zero )
732 $ dif = ulpinv
733 dif = max( dif, abs( ssav( i )-s( i ) ) / div )
734 80 CONTINUE
735 result( 7 ) = max( result( 7 ), dif )
736 90 CONTINUE
737 100 CONTINUE
738*
739* Test for CGESDD
740*
741 iwtmp = 2*mnmin*mnmin + 2*mnmin + max( m, n )
742 lswork = iwtmp + ( iwspc-1 )*( lwork-iwtmp ) / 3
743 lswork = min( lswork, lwork )
744 lswork = max( lswork, 1 )
745 IF( iwspc.EQ.4 )
746 $ lswork = lwork
747*
748* Factorize A
749*
750 CALL clacpy( 'F', m, n, asav, lda, a, lda )
751 srnamt = 'CGESDD'
752 CALL cgesdd( 'A', m, n, a, lda, ssav, usav, ldu, vtsav,
753 $ ldvt, work, lswork, rwork, iwork, iinfo )
754 IF( iinfo.NE.0 ) THEN
755 WRITE( nounit, fmt = 9995 )'GESDD', iinfo, m, n,
756 $ jtype, lswork, ioldsd
757 info = abs( iinfo )
758 RETURN
759 END IF
760*
761* Do tests 1--4
762*
763 CALL cbdt01( m, n, 0, asav, lda, usav, ldu, ssav, e,
764 $ vtsav, ldvt, work, rwork, result( 8 ) )
765 IF( m.NE.0 .AND. n.NE.0 ) THEN
766 CALL cunt01( 'Columns', mnmin, m, usav, ldu, work,
767 $ lwork, rwork, result( 9 ) )
768 CALL cunt01( 'Rows', mnmin, n, vtsav, ldvt, work,
769 $ lwork, rwork, result( 10 ) )
770 END IF
771 result( 11 ) = 0
772 DO 110 i = 1, mnmin - 1
773 IF( ssav( i ).LT.ssav( i+1 ) )
774 $ result( 11 ) = ulpinv
775 IF( ssav( i ).LT.zero )
776 $ result( 11 ) = ulpinv
777 110 CONTINUE
778 IF( mnmin.GE.1 ) THEN
779 IF( ssav( mnmin ).LT.zero )
780 $ result( 11 ) = ulpinv
781 END IF
782*
783* Do partial SVDs, comparing to SSAV, USAV, and VTSAV
784*
785 result( 12 ) = zero
786 result( 13 ) = zero
787 result( 14 ) = zero
788 DO 130 ijq = 0, 2
789 jobq = cjob( ijq+1 )
790 CALL clacpy( 'F', m, n, asav, lda, a, lda )
791 srnamt = 'CGESDD'
792 CALL cgesdd( jobq, m, n, a, lda, s, u, ldu, vt, ldvt,
793 $ work, lswork, rwork, iwork, iinfo )
794*
795* Compare U
796*
797 dif = zero
798 IF( m.GT.0 .AND. n.GT.0 ) THEN
799 IF( ijq.EQ.1 ) THEN
800 IF( m.GE.n ) THEN
801 CALL cunt03( 'C', m, mnmin, m, mnmin, usav,
802 $ ldu, a, lda, work, lwork, rwork,
803 $ dif, iinfo )
804 ELSE
805 CALL cunt03( 'C', m, mnmin, m, mnmin, usav,
806 $ ldu, u, ldu, work, lwork, rwork,
807 $ dif, iinfo )
808 END IF
809 ELSE IF( ijq.EQ.2 ) THEN
810 CALL cunt03( 'C', m, mnmin, m, mnmin, usav, ldu,
811 $ u, ldu, work, lwork, rwork, dif,
812 $ iinfo )
813 END IF
814 END IF
815 result( 12 ) = max( result( 12 ), dif )
816*
817* Compare VT
818*
819 dif = zero
820 IF( m.GT.0 .AND. n.GT.0 ) THEN
821 IF( ijq.EQ.1 ) THEN
822 IF( m.GE.n ) THEN
823 CALL cunt03( 'r', N, MNMIN, N, MNMIN, VTSAV,
824 $ LDVT, VT, LDVT, WORK, LWORK,
825 $ RWORK, DIF, IINFO )
826 ELSE
827 CALL CUNT03( 'r', N, MNMIN, N, MNMIN, VTSAV,
828 $ LDVT, A, LDA, WORK, LWORK,
829 $ RWORK, DIF, IINFO )
830 END IF
831.EQ. ELSE IF( IJQ2 ) THEN
832 CALL CUNT03( 'r', N, MNMIN, N, MNMIN, VTSAV,
833 $ LDVT, VT, LDVT, WORK, LWORK, RWORK,
834 $ DIF, IINFO )
835 END IF
836 END IF
837 RESULT( 13 ) = MAX( RESULT( 13 ), DIF )
838*
839* Compare S
840*
841 DIF = ZERO
842 DIV = MAX( REAL( MNMIN )*ULP*S( 1 ),
843 $ SLAMCH( 'safe minimum' ) )
844 DO 120 I = 1, MNMIN - 1
845.LT. IF( SSAV( I )SSAV( I+1 ) )
846 $ DIF = ULPINV
847.LT. IF( SSAV( I )ZERO )
848 $ DIF = ULPINV
849 DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
850 120 CONTINUE
851 RESULT( 14 ) = MAX( RESULT( 14 ), DIF )
852 130 CONTINUE
853
854*
855* Test CGESVDQ
856* Note: CGESVDQ only works for M >= N
857*
858 RESULT( 36 ) = ZERO
859 RESULT( 37 ) = ZERO
860 RESULT( 38 ) = ZERO
861 RESULT( 39 ) = ZERO
862*
863.GE. IF( MN ) THEN
864 IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
865 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
866 LSWORK = MIN( LSWORK, LWORK )
867 LSWORK = MAX( LSWORK, 1 )
868.EQ. IF( IWSPC4 )
869 $ LSWORK = LWORK
870*
871 CALL CLACPY( 'f', M, N, ASAV, LDA, A, LDA )
872 SRNAMT = 'cgesvdq'
873*
874 LRWORK = MAX(2, M, 5*N)
875 LIWORK = MAX( N, 1 )
876 CALL CGESVDQ( 'h', 'n', 'n', 'a', 'a',
877 $ M, N, A, LDA, SSAV, USAV, LDU,
878 $ VTSAV, LDVT, NUMRANK, IWORK, LIWORK,
879 $ WORK, LWORK, RWORK, LRWORK, IINFO )
880*
881.NE. IF( IINFO0 ) THEN
882 WRITE( NOUNIT, FMT = 9995 )'cgesvdq', IINFO, M, N,
883 $ JTYPE, LSWORK, IOLDSD
884 INFO = ABS( IINFO )
885 RETURN
886 END IF
887*
888* Do tests 36--39
889*
890 CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
891 $ VTSAV, LDVT, WORK, RWORK, RESULT( 36 ) )
892.NE..AND..NE. IF( M0 N0 ) THEN
893 CALL CUNT01( 'columns', M, M, USAV, LDU, WORK,
894 $ LWORK, RWORK, RESULT( 37 ) )
895 CALL CUNT01( 'rows', N, N, VTSAV, LDVT, WORK,
896 $ LWORK, RWORK, RESULT( 38 ) )
897 END IF
898 RESULT( 39 ) = ZERO
899 DO 199 I = 1, MNMIN - 1
900.LT. IF( SSAV( I )SSAV( I+1 ) )
901 $ RESULT( 39 ) = ULPINV
902.LT. IF( SSAV( I )ZERO )
903 $ RESULT( 39 ) = ULPINV
904 199 CONTINUE
905.GE. IF( MNMIN1 ) THEN
906.LT. IF( SSAV( MNMIN )ZERO )
907 $ RESULT( 39 ) = ULPINV
908 END IF
909 END IF
910*
911* Test CGESVJ
912* Note: CGESVJ only works for M >= N
913*
914 RESULT( 15 ) = ZERO
915 RESULT( 16 ) = ZERO
916 RESULT( 17 ) = ZERO
917 RESULT( 18 ) = ZERO
918*
919.GE. IF( MN ) THEN
920 IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
921 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
922 LSWORK = MIN( LSWORK, LWORK )
923 LSWORK = MAX( LSWORK, 1 )
924 LRWORK = MAX(6,N)
925.EQ. IF( IWSPC4 )
926 $ LSWORK = LWORK
927*
928 CALL CLACPY( 'f', M, N, ASAV, LDA, USAV, LDA )
929 SRNAMT = 'cgesvj'
930 CALL CGESVJ( 'g', 'u', 'v', M, N, USAV, LDA, SSAV,
931 & 0, A, LDVT, WORK, LWORK, RWORK,
932 & LRWORK, IINFO )
933*
934* CGESVJ returns V not VH
935*
936 DO J=1,N
937 DO I=1,N
938 VTSAV(J,I) = CONJG (A(I,J))
939 END DO
940 END DO
941*
942.NE. IF( IINFO0 ) THEN
943 WRITE( NOUNIT, FMT = 9995 )'gesvj', IINFO, M, N,
944 $ JTYPE, LSWORK, IOLDSD
945 INFO = ABS( IINFO )
946 RETURN
947 END IF
948*
949* Do tests 15--18
950*
951 CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
952 $ VTSAV, LDVT, WORK, RWORK, RESULT( 15 ) )
953.NE..AND..NE. IF( M0 N0 ) THEN
954 CALL CUNT01( 'columns', M, M, USAV, LDU, WORK,
955 $ LWORK, RWORK, RESULT( 16 ) )
956 CALL CUNT01( 'rows', N, N, VTSAV, LDVT, WORK,
957 $ LWORK, RWORK, RESULT( 17 ) )
958 END IF
959 RESULT( 18 ) = ZERO
960 DO 131 I = 1, MNMIN - 1
961.LT. IF( SSAV( I )SSAV( I+1 ) )
962 $ RESULT( 18 ) = ULPINV
963.LT. IF( SSAV( I )ZERO )
964 $ RESULT( 18 ) = ULPINV
965 131 CONTINUE
966.GE. IF( MNMIN1 ) THEN
967.LT. IF( SSAV( MNMIN )ZERO )
968 $ RESULT( 18 ) = ULPINV
969 END IF
970 END IF
971*
972* Test CGEJSV
973* Note: CGEJSV only works for M >= N
974*
975 RESULT( 19 ) = ZERO
976 RESULT( 20 ) = ZERO
977 RESULT( 21 ) = ZERO
978 RESULT( 22 ) = ZERO
979.GE. IF( MN ) THEN
980 IWTMP = 2*MNMIN*MNMIN + 2*MNMIN + MAX( M, N )
981 LSWORK = IWTMP + ( IWSPC-1 )*( LWORK-IWTMP ) / 3
982 LSWORK = MIN( LSWORK, LWORK )
983 LSWORK = MAX( LSWORK, 1 )
984.EQ. IF( IWSPC4 )
985 $ LSWORK = LWORK
986 LRWORK = MAX( 7, N + 2*M)
987*
988 CALL CLACPY( 'f', M, N, ASAV, LDA, VTSAV, LDA )
989 SRNAMT = 'cgejsv'
990 CALL CGEJSV( 'g', 'u', 'v', 'r', 'n', 'n',
991 & M, N, VTSAV, LDA, SSAV, USAV, LDU, A, LDVT,
992 & WORK, LWORK, RWORK,
993 & LRWORK, IWORK, IINFO )
994*
995* CGEJSV returns V not VH
996*
997 DO 133 J=1,N
998 DO 132 I=1,N
999 VTSAV(J,I) = CONJG (A(I,J))
1000 132 END DO
1001 133 END DO
1002*
1003.NE. IF( IINFO0 ) THEN
1004 WRITE( NOUNIT, FMT = 9995 )'gejsv', IINFO, M, N,
1005 $ JTYPE, LSWORK, IOLDSD
1006 INFO = ABS( IINFO )
1007 RETURN
1008 END IF
1009*
1010* Do tests 19--22
1011*
1012 CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
1013 $ VTSAV, LDVT, WORK, RWORK, RESULT( 19 ) )
1014.NE..AND..NE. IF( M0 N0 ) THEN
1015 CALL CUNT01( 'columns', M, M, USAV, LDU, WORK,
1016 $ LWORK, RWORK, RESULT( 20 ) )
1017 CALL CUNT01( 'rows', N, N, VTSAV, LDVT, WORK,
1018 $ LWORK, RWORK, RESULT( 21 ) )
1019 END IF
1020 RESULT( 22 ) = ZERO
1021 DO 134 I = 1, MNMIN - 1
1022.LT. IF( SSAV( I )SSAV( I+1 ) )
1023 $ RESULT( 22 ) = ULPINV
1024.LT. IF( SSAV( I )ZERO )
1025 $ RESULT( 22 ) = ULPINV
1026 134 CONTINUE
1027.GE. IF( MNMIN1 ) THEN
1028.LT. IF( SSAV( MNMIN )ZERO )
1029 $ RESULT( 22 ) = ULPINV
1030 END IF
1031 END IF
1032*
1033* Test CGESVDX
1034*
1035* Factorize A
1036*
1037 CALL CLACPY( 'f', M, N, ASAV, LDA, A, LDA )
1038 SRNAMT = 'cgesvdx'
1039 CALL CGESVDX( 'v', 'v', 'a', M, N, A, LDA,
1040 $ VL, VU, IL, IU, NS, SSAV, USAV, LDU,
1041 $ VTSAV, LDVT, WORK, LWORK, RWORK,
1042 $ IWORK, IINFO )
1043.NE. IF( IINFO0 ) THEN
1044 WRITE( NOUNIT, FMT = 9995 )'gesvdx', IINFO, M, N,
1045 $ JTYPE, LSWORK, IOLDSD
1046 INFO = ABS( IINFO )
1047 RETURN
1048 END IF
1049*
1050* Do tests 1--4
1051*
1052 RESULT( 23 ) = ZERO
1053 RESULT( 24 ) = ZERO
1054 RESULT( 25 ) = ZERO
1055 CALL CBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
1056 $ VTSAV, LDVT, WORK, RWORK, RESULT( 23 ) )
1057.NE..AND..NE. IF( M0 N0 ) THEN
1058 CALL CUNT01( 'columns', MNMIN, M, USAV, LDU, WORK,
1059 $ LWORK, RWORK, RESULT( 24 ) )
1060 CALL CUNT01( 'rows', MNMIN, N, VTSAV, LDVT, WORK,
1061 $ LWORK, RWORK, RESULT( 25 ) )
1062 END IF
1063 RESULT( 26 ) = ZERO
1064 DO 140 I = 1, MNMIN - 1
1065.LT. IF( SSAV( I )SSAV( I+1 ) )
1066 $ RESULT( 26 ) = ULPINV
1067.LT. IF( SSAV( I )ZERO )
1068 $ RESULT( 26 ) = ULPINV
1069 140 CONTINUE
1070.GE. IF( MNMIN1 ) THEN
1071.LT. IF( SSAV( MNMIN )ZERO )
1072 $ RESULT( 26 ) = ULPINV
1073 END IF
1074*
1075* Do partial SVDs, comparing to SSAV, USAV, and VTSAV
1076*
1077 RESULT( 27 ) = ZERO
1078 RESULT( 28 ) = ZERO
1079 RESULT( 29 ) = ZERO
1080 DO 170 IJU = 0, 1
1081 DO 160 IJVT = 0, 1
1082.EQ..AND..EQ..OR. IF( ( IJU0 IJVT0 )
1083.EQ..AND..EQ. $ ( IJU1 IJVT1 ) ) GO TO 160
1084 JOBU = CJOBV( IJU+1 )
1085 JOBVT = CJOBV( IJVT+1 )
1086 RANGE = CJOBR( 1 )
1087 CALL CLACPY( 'f', M, N, ASAV, LDA, A, LDA )
1088 SRNAMT = 'cgesvdx'
1089 CALL CGESVDX( JOBU, JOBVT, 'a', M, N, A, LDA,
1090 $ VL, VU, IL, IU, NS, SSAV, U, LDU,
1091 $ VT, LDVT, WORK, LWORK, RWORK,
1092 $ IWORK, IINFO )
1093*
1094* Compare U
1095*
1096 DIF = ZERO
1097.GT..AND..GT. IF( M0 N0 ) THEN
1098.EQ. IF( IJU1 ) THEN
1099 CALL CUNT03( 'c', M, MNMIN, M, MNMIN, USAV,
1100 $ LDU, U, LDU, WORK, LWORK, RWORK,
1101 $ DIF, IINFO )
1102 END IF
1103 END IF
1104 RESULT( 27 ) = MAX( RESULT( 27 ), DIF )
1105*
1106* Compare VT
1107*
1108 DIF = ZERO
1109.GT..AND..GT. IF( M0 N0 ) THEN
1110.EQ. IF( IJVT1 ) THEN
1111 CALL CUNT03( 'r', N, MNMIN, N, MNMIN, VTSAV,
1112 $ LDVT, VT, LDVT, WORK, LWORK,
1113 $ RWORK, DIF, IINFO )
1114 END IF
1115 END IF
1116 RESULT( 28 ) = MAX( RESULT( 28 ), DIF )
1117*
1118* Compare S
1119*
1120 DIF = ZERO
1121 DIV = MAX( REAL( MNMIN )*ULP*S( 1 ),
1122 $ SLAMCH( 'safe minimum' ) )
1123 DO 150 I = 1, MNMIN - 1
1124.LT. IF( SSAV( I )SSAV( I+1 ) )
1125 $ DIF = ULPINV
1126.LT. IF( SSAV( I )ZERO )
1127 $ DIF = ULPINV
1128 DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
1129 150 CONTINUE
1130 RESULT( 29) = MAX( RESULT( 29 ), DIF )
1131 160 CONTINUE
1132 170 CONTINUE
1133*
1134* Do tests 8--10
1135*
1136 DO 180 I = 1, 4
1137 ISEED2( I ) = ISEED( I )
1138 180 CONTINUE
1139.LE. IF( MNMIN1 ) THEN
1140 IL = 1
1141 IU = MAX( 1, MNMIN )
1142 ELSE
1143 IL = 1 + INT( ( MNMIN-1 )*SLARND( 1, ISEED2 ) )
1144 IU = 1 + INT( ( MNMIN-1 )*SLARND( 1, ISEED2 ) )
1145.LT. IF( IUIL ) THEN
1146 ITEMP = IU
1147 IU = IL
1148 IL = ITEMP
1149 END IF
1150 END IF
1151 CALL CLACPY( 'f', M, N, ASAV, LDA, A, LDA )
1152 SRNAMT = 'cgesvdx'
1153 CALL CGESVDX( 'v', 'v', 'i', M, N, A, LDA,
1154 $ VL, VU, IL, IU, NSI, S, U, LDU,
1155 $ VT, LDVT, WORK, LWORK, RWORK,
1156 $ IWORK, IINFO )
1157.NE. IF( IINFO0 ) THEN
1158 WRITE( NOUNIT, FMT = 9995 )'gesvdx', IINFO, M, N,
1159 $ JTYPE, LSWORK, IOLDSD
1160 INFO = ABS( IINFO )
1161 RETURN
1162 END IF
1163*
1164 RESULT( 30 ) = ZERO
1165 RESULT( 31 ) = ZERO
1166 RESULT( 32 ) = ZERO
1167 CALL CBDT05( M, N, ASAV, LDA, S, NSI, U, LDU,
1168 $ VT, LDVT, WORK, RESULT( 30 ) )
1169.NE..AND..NE. IF( M0 N0 ) THEN
1170 CALL CUNT01( 'columns', M, NSI, U, LDU, WORK,
1171 $ LWORK, RWORK, RESULT( 31 ) )
1172 CALL CUNT01( 'rows', nsi, n, vt, ldvt, work,
1173 $ lwork, rwork, result( 32 ) )
1174 END IF
1175*
1176* Do tests 11--13
1177*
1178 IF( mnmin.GT.0 .AND. nsi.GT.1 ) THEN
1179 IF( il.NE.1 ) THEN
1180 vu = ssav( il ) +
1181 $ max( half*abs( ssav( il )-ssav( il-1 ) ),
1182 $ ulp*anorm, two*rtunfl )
1183 ELSE
1184 vu = ssav( 1 ) +
1185 $ max( half*abs( ssav( ns )-ssav( 1 ) ),
1186 $ ulp*anorm, two*rtunfl )
1187 END IF
1188 IF( iu.NE.ns ) THEN
1189 vl = ssav( iu ) - max( ulp*anorm, two*rtunfl,
1190 $ half*abs( ssav( iu+1 )-ssav( iu ) ) )
1191 ELSE
1192 vl = ssav( ns ) - max( ulp*anorm, two*rtunfl,
1193 $ half*abs( ssav( ns )-ssav( 1 ) ) )
1194 END IF
1195 vl = max( vl,zero )
1196 vu = max( vu,zero )
1197 IF( vl.GE.vu ) vu = max( vu*2, vu+vl+half )
1198 ELSE
1199 vl = zero
1200 vu = one
1201 END IF
1202 CALL clacpy( 'F', m, n, asav, lda, a, lda )
1203 srnamt = 'CGESVDX'
1204 CALL cgesvdx( 'V', 'V', 'V', m, n, a, lda,
1205 $ vl, vu, il, iu, nsv, s, u, ldu,
1206 $ vt, ldvt, work, lwork, rwork,
1207 $ iwork, iinfo )
1208 IF( iinfo.NE.0 ) THEN
1209 WRITE( nounit, fmt = 9995 )'GESVDX', iinfo, m, n,
1210 $ jtype, lswork, ioldsd
1211 info = abs( iinfo )
1212 RETURN
1213 END IF
1214*
1215 result( 33 ) = zero
1216 result( 34 ) = zero
1217 result( 35 ) = zero
1218 CALL cbdt05( m, n, asav, lda, s, nsv, u, ldu,
1219 $ vt, ldvt, work, result( 33 ) )
1220 IF( m.NE.0 .AND. n.NE.0 ) THEN
1221 CALL cunt01( 'Columns', m, nsv, u, ldu, work,
1222 $ lwork, rwork, result( 34 ) )
1223 CALL cunt01( 'Rows', nsv, n, vt, ldvt, work,
1224 $ lwork, rwork, result( 35 ) )
1225 END IF
1226*
1227* End of Loop -- Check for RESULT(j) > THRESH
1228*
1229 ntest = 0
1230 nfail = 0
1231 DO 190 j = 1, 39
1232 IF( result( j ).GE.zero )
1233 $ ntest = ntest + 1
1234 IF( result( j ).GE.thresh )
1235 $ nfail = nfail + 1
1236 190 CONTINUE
1237*
1238 IF( nfail.GT.0 )
1239 $ ntestf = ntestf + 1
1240 IF( ntestf.EQ.1 ) THEN
1241 WRITE( nounit, fmt = 9999 )
1242 WRITE( nounit, fmt = 9998 )thresh
1243 ntestf = 2
1244 END IF
1245*
1246 DO 200 j = 1, 39
1247 IF( result( j ).GE.thresh ) THEN
1248 WRITE( nounit, fmt = 9997 )m, n, jtype, iwspc,
1249 $ ioldsd, j, result( j )
1250 END IF
1251 200 CONTINUE
1252*
1253 nerrs = nerrs + nfail
1254 ntestt = ntestt + ntest
1255*
1256 290 CONTINUE
1257*
1258 300 CONTINUE
1259 310 CONTINUE
1260*
1261* Summary
1262*
1263 CALL alasvm( 'CBD', nounit, nerrs, ntestt, 0 )
1264*
1265 9999 FORMAT( ' SVD -- Complex Singular Value Decomposition Driver ',
1266 $ / ' Matrix types (see CDRVBD for details):',
1267 $ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
1268 $ / ' 3 = Evenly spaced singular values near 1',
1269 $ / ' 4 = Evenly spaced singular values near underflow',
1270 $ / ' 5 = Evenly spaced singular values near overflow',
1271 $ / / ' Tests performed: ( A is dense, U and V are unitary,',
1272 $ / 19x, ' S is an array, and Upartial, VTpartial, and',
1273 $ / 19x, ' Spartial are partially computed U, VT and S),', / )
1274 9998 FORMAT( ' Tests performed with Test Threshold = ', f8.2,
1275 $ / ' CGESVD: ', /
1276 $ ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1277 $ / ' 2 = | I - U**T U | / ( M ulp ) ',
1278 $ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
1279 $ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
1280 $ ' decreasing order, else 1/ulp',
1281 $ / ' 5 = | U - Upartial | / ( M ulp )',
1282 $ / ' 6 = | VT - VTpartial | / ( N ulp )',
1283 $ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
1284 $ / ' CGESDD: ', /
1285 $ ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1286 $ / ' 9 = | I - U**T U | / ( M ulp ) ',
1287 $ / '10 = | I - VT VT**T | / ( N ulp ) ',
1288 $ / '11 = 0 if S contains min(M,N) nonnegative values in',
1289 $ ' decreasing order, else 1/ulp',
1290 $ / '12 = | U - Upartial | / ( M ulp )',
1291 $ / '13 = | VT - VTpartial | / ( N ulp )',
1292 $ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )',
1293 $ / ' CGESVJ: ', /
1294 $ / '15 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1295 $ / '16 = | I - U**T U | / ( M ulp ) ',
1296 $ / '17 = | I - VT VT**T | / ( N ulp ) ',
1297 $ / '18 = 0 if S contains min(M,N) nonnegative values in',
1298 $ ' decreasing order, else 1/ulp',
1299 $ / ' CGESJV: ', /
1300 $ / '19 = | A - U diag(S) VT | / ( |A| max(M,N) ulp )',
1301 $ / '20 = | I - U**T U | / ( M ulp ) ',
1302 $ / '21 = | I - VT VT**T | / ( N ulp ) ',
1303 $ / '22 = 0 if S contains min(M,N) nonnegative values in',
1304 $ ' decreasing order, else 1/ulp',
1305 $ / ' CGESVDX(V,V,A): ', /
1306 $ '23 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1307 $ / '24 = | I - U**T U | / ( M ulp ) ',
1308 $ / '25 = | I - VT VT**T | / ( N ulp ) ',
1309 $ / '26 = 0 if S contains min(M,N) nonnegative values in',
1310 $ ' decreasing order, else 1/ulp',
1311 $ / '27 = | U - Upartial | / ( M ulp )',
1312 $ / '28 = | VT - VTpartial | / ( N ulp )',
1313 $ / '29 = | S - Spartial | / ( min(M,N) ulp |S| )',
1314 $ / ' CGESVDX(V,V,I): ',
1315 $ / '30 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp )',
1316 $ / '31 = | I - U**T U | / ( M ulp ) ',
1317 $ / '32 = | I - VT VT**T | / ( N ulp ) ',
1318 $ / ' CGESVDX(V,V,V) ',
1319 $ / '33 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp )',
1320 $ / '34 = | I - U**T U | / ( M ulp ) ',
1321 $ / '35 = | I - VT VT**T | / ( N ulp ) ',
1322 $ ' CGESVDQ(H,N,N,A,A',
1323 $ / '36 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1324 $ / '37 = | I - U**T U | / ( M ulp ) ',
1325 $ / '38 = | I - VT VT**T | / ( N ulp ) ',
1326 $ / '39 = 0 if S contains min(M,N) nonnegative values in',
1327 $ ' decreasing order, else 1/ulp',
1328 $ / / )
1329 9997 FORMAT( ' M=', i5, ', N=', i5, ', type ', i1, ', IWS=', i1,
1330 $ ', seed=', 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1331 9996 FORMAT( ' CDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1332 $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1333 $ i5, ')' )
1334 9995 FORMAT( ' CDRVBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1335 $ i6, ', N=', i6, ', JTYPE=', i6, ', LSWORK=', i6, / 9x,
1336 $ 'ISEED=(', 3( i5, ',' ), i5, ')' )
1337*
1338 RETURN
1339*
1340* End of CDRVBD
1341*
1342 END
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
alasvm
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
cgesvj
subroutine cgesvj(joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv, cwork, lwork, rwork, lrwork, info)
CGESVJ
Definition cgesvj.f:351
cgesvdx
subroutine cgesvdx(jobu, jobvt, range, m, n, a, lda, vl, vu, il, iu, ns, s, u, ldu, vt, ldvt, work, lwork, rwork, iwork, info)
CGESVDX computes the singular value decomposition (SVD) for GE matrices
Definition cgesvdx.f:270
cgesvd
subroutine cgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, info)
CGESVD computes the singular value decomposition (SVD) for GE matrices
Definition cgesvd.f:214
cgesdd
subroutine cgesdd(jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, iwork, info)
CGESDD
Definition cgesdd.f:227
cgejsv
subroutine cgejsv(joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, cwork, lwork, rwork, lrwork, iwork, info)
CGEJSV
Definition cgejsv.f:568
cgesvdq
subroutine cgesvdq(joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, cwork, lcwork, rwork, lrwork, info)
CGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE...
Definition cgesvdq.f:413
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
cbdt01
subroutine cbdt01(m, n, kd, a, lda, q, ldq, d, e, pt, ldpt, work, rwork, resid)
CBDT01
Definition cbdt01.f:147
cunt03
subroutine cunt03(rc, mu, mv, n, k, u, ldu, v, ldv, work, lwork, rwork, result, info)
CUNT03
Definition cunt03.f:162
cunt01
subroutine cunt01(rowcol, m, n, u, ldu, work, lwork, rwork, resid)
CUNT01
Definition cunt01.f:126
cdrvbd
subroutine cdrvbd(nsizes, mm, nn, ntypes, dotype, iseed, thresh, a, lda, u, ldu, vt, ldvt, asav, usav, vtsav, s, ssav, e, work, lwork, rwork, iwork, nounit, info)
CDRVBD
Definition cdrvbd.f:401
clatms
subroutine clatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
CLATMS
Definition clatms.f:332
cbdt05
subroutine cbdt05(m, n, a, lda, s, ns, u, ldu, vt, ldvt, work, resid)
CBDT05
Definition cbdt05.f:125
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21