1*> \brief \b DDRVBD
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 DDRVBD( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
12*                          A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S,
13*                          SSAV, E, WORK, LWORK, IWORK, NOUT, INFO )
14*
15*       .. Scalar Arguments ..
16*       INTEGER            INFO, LDA, LDU, LDVT, LWORK, NOUT, NSIZES,
17*      $                   NTYPES
18*       DOUBLE PRECISION   THRESH
19*       ..
20*       .. Array Arguments ..
21*       LOGICAL            DOTYPE( * )
22*       INTEGER            ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
23*       DOUBLE PRECISION   A( LDA, * ), ASAV( LDA, * ), E( * ), S( * ),
24*      $                   SSAV( * ), U( LDU, * ), USAV( LDU, * ),
25*      $                   VT( LDVT, * ), VTSAV( LDVT, * ), WORK( * )
26*       ..
27*
28*
29*> \par Purpose:
30*  =============
31*>
32*> \verbatim
33*>
34*> DDRVBD checks the singular value decomposition (SVD) drivers
35*> DGESVD, DGESDD, DGESVDQ, DGESVJ, DGEJSV, and DGESVDX.
36*>
37*> Both DGESVD and DGESDD factor A = U diag(S) VT, where U and VT are
38*> orthogonal and diag(S) is diagonal with the entries of the array S
39*> on its diagonal. The entries of S are the singular values,
40*> nonnegative and stored in decreasing order.  U and VT can be
41*> optionally not computed, overwritten on A, or computed partially.
42*>
43*> A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN.
44*> U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N.
45*>
46*> When DDRVBD is called, a number of matrix "sizes" (M's and N's)
47*> and a number of matrix "types" are specified.  For each size (M,N)
48*> and each type of matrix, and for the minimal workspace as well as
49*> workspace adequate to permit blocking, an  M x N  matrix "A" will be
50*> generated and used to test the SVD routines.  For each matrix, A will
51*> be factored as A = U diag(S) VT and the following 12 tests computed:
52*>
53*> Test for DGESVD:
54*>
55*> (1)    | A - U diag(S) VT | / ( |A| max(M,N) ulp )
56*>
57*> (2)    | I - U'U | / ( M ulp )
58*>
59*> (3)    | I - VT VT' | / ( N ulp )
60*>
61*> (4)    S contains MNMIN nonnegative values in decreasing order.
62*>        (Return 0 if true, 1/ULP if false.)
63*>
64*> (5)    | U - Upartial | / ( M ulp ) where Upartial is a partially
65*>        computed U.
66*>
67*> (6)    | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
68*>        computed VT.
69*>
70*> (7)    | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
71*>        vector of singular values from the partial SVD
72*>
73*> Test for DGESDD:
74*>
75*> (8)    | A - U diag(S) VT | / ( |A| max(M,N) ulp )
76*>
77*> (9)    | I - U'U | / ( M ulp )
78*>
79*> (10)   | I - VT VT' | / ( N ulp )
80*>
81*> (11)   S contains MNMIN nonnegative values in decreasing order.
82*>        (Return 0 if true, 1/ULP if false.)
83*>
84*> (12)   | U - Upartial | / ( M ulp ) where Upartial is a partially
85*>        computed U.
86*>
87*> (13)   | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
88*>        computed VT.
89*>
90*> (14)   | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
91*>        vector of singular values from the partial SVD
92*>
93*> Test for DGESVDQ:
94*>
95*> (36)   | A - U diag(S) VT | / ( |A| max(M,N) ulp )
96*>
97*> (37)   | I - U'U | / ( M ulp )
98*>
99*> (38)   | I - VT VT' | / ( N ulp )
100*>
101*> (39)   S contains MNMIN nonnegative values in decreasing order.
102*>        (Return 0 if true, 1/ULP if false.)
103*>
104*> Test for DGESVJ:
105*>
106*> (15)   | A - U diag(S) VT | / ( |A| max(M,N) ulp )
107*>
108*> (16)   | I - U'U | / ( M ulp )
109*>
110*> (17)   | I - VT VT' | / ( N ulp )
111*>
112*> (18)   S contains MNMIN nonnegative values in decreasing order.
113*>        (Return 0 if true, 1/ULP if false.)
114*>
115*> Test for DGEJSV:
116*>
117*> (19)   | A - U diag(S) VT | / ( |A| max(M,N) ulp )
118*>
119*> (20)   | I - U'U | / ( M ulp )
120*>
121*> (21)   | I - VT VT' | / ( N ulp )
122*>
123*> (22)   S contains MNMIN nonnegative values in decreasing order.
124*>        (Return 0 if true, 1/ULP if false.)
125*>
126*> Test for DGESVDX( 'V', 'V', 'A' )/DGESVDX( 'N', 'N', 'A' )
127*>
128*> (23)   | A - U diag(S) VT | / ( |A| max(M,N) ulp )
129*>
130*> (24)   | I - U'U | / ( M ulp )
131*>
132*> (25)   | I - VT VT' | / ( N ulp )
133*>
134*> (26)   S contains MNMIN nonnegative values in decreasing order.
135*>        (Return 0 if true, 1/ULP if false.)
136*>
137*> (27)   | U - Upartial | / ( M ulp ) where Upartial is a partially
138*>        computed U.
139*>
140*> (28)   | VT - VTpartial | / ( N ulp ) where VTpartial is a partially
141*>        computed VT.
142*>
143*> (29)   | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the
144*>        vector of singular values from the partial SVD
145*>
146*> Test for DGESVDX( 'V', 'V', 'I' )
147*>
148*> (30)   | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
149*>
150*> (31)   | I - U'U | / ( M ulp )
151*>
152*> (32)   | I - VT VT' | / ( N ulp )
153*>
154*> Test for DGESVDX( 'V', 'V', 'V' )
155*>
156*> (33)   | U' A VT''' - diag(S) | / ( |A| max(M,N) ulp )
157*>
158*> (34)   | I - U'U | / ( M ulp )
159*>
160*> (35)   | I - VT VT' | / ( N ulp )
161*>
162*> The "sizes" are specified by the arrays MM(1:NSIZES) and
163*> NN(1:NSIZES); the value of each element pair (MM(j),NN(j))
164*> specifies one size.  The "types" are specified by a logical array
165*> DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j"
166*> will be generated.
167*> Currently, the list of possible types is:
168*>
169*> (1)  The zero matrix.
170*> (2)  The identity matrix.
171*> (3)  A matrix of the form  U D V, where U and V are orthogonal and
172*>      D has evenly spaced entries 1, ..., ULP with random signs
173*>      on the diagonal.
174*> (4)  Same as (3), but multiplied by the underflow-threshold / ULP.
175*> (5)  Same as (3), but multiplied by the overflow-threshold * ULP.
176*> \endverbatim
177*
178*  Arguments:
179*  ==========
180*
181*> \param[in] NSIZES
182*> \verbatim
183*>          NSIZES is INTEGER
184*>          The number of matrix sizes (M,N) contained in the vectors
185*>          MM and NN.
186*> \endverbatim
187*>
188*> \param[in] MM
189*> \verbatim
190*>          MM is INTEGER array, dimension (NSIZES)
191*>          The values of the matrix row dimension M.
192*> \endverbatim
193*>
194*> \param[in] NN
195*> \verbatim
196*>          NN is INTEGER array, dimension (NSIZES)
197*>          The values of the matrix column dimension N.
198*> \endverbatim
199*>
200*> \param[in] NTYPES
201*> \verbatim
202*>          NTYPES is INTEGER
203*>          The number of elements in DOTYPE.   If it is zero, DDRVBD
204*>          does nothing.  It must be at least zero.  If it is MAXTYP+1
205*>          and NSIZES is 1, then an additional type, MAXTYP+1 is
206*>          defined, which is to use whatever matrices are in A and B.
207*>          This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
208*>          DOTYPE(MAXTYP+1) is .TRUE. .
209*> \endverbatim
210*>
211*> \param[in] DOTYPE
212*> \verbatim
213*>          DOTYPE is LOGICAL array, dimension (NTYPES)
214*>          If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
215*>          of type j will be generated.  If NTYPES is smaller than the
216*>          maximum number of types defined (PARAMETER MAXTYP), then
217*>          types NTYPES+1 through MAXTYP will not be generated.  If
218*>          NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
219*>          DOTYPE(NTYPES) will be ignored.
220*> \endverbatim
221*>
222*> \param[in,out] ISEED
223*> \verbatim
224*>          ISEED is INTEGER array, dimension (4)
225*>          On entry, the seed of the random number generator.  The array
226*>          elements should be between 0 and 4095; if not they will be
227*>          reduced mod 4096.  Also, ISEED(4) must be odd.
228*>          On exit, ISEED is changed and can be used in the next call to
229*>          DDRVBD to continue the same random number sequence.
230*> \endverbatim
231*>
232*> \param[in] THRESH
233*> \verbatim
234*>          THRESH is DOUBLE PRECISION
235*>          The threshold value for the test ratios.  A result is
236*>          included in the output file if RESULT >= THRESH.  The test
237*>          ratios are scaled to be O(1), so THRESH should be a small
238*>          multiple of 1, e.g., 10 or 100.  To have every test ratio
239*>          printed, use THRESH = 0.
240*> \endverbatim
241*>
242*> \param[out] A
243*> \verbatim
244*>          A is DOUBLE PRECISION array, dimension (LDA,NMAX)
245*>          where NMAX is the maximum value of N in NN.
246*> \endverbatim
247*>
248*> \param[in] LDA
249*> \verbatim
250*>          LDA is INTEGER
251*>          The leading dimension of the array A.  LDA >= max(1,MMAX),
252*>          where MMAX is the maximum value of M in MM.
253*> \endverbatim
254*>
255*> \param[out] U
256*> \verbatim
257*>          U is DOUBLE PRECISION array, dimension (LDU,MMAX)
258*> \endverbatim
259*>
260*> \param[in] LDU
261*> \verbatim
262*>          LDU is INTEGER
263*>          The leading dimension of the array U.  LDU >= max(1,MMAX).
264*> \endverbatim
265*>
266*> \param[out] VT
267*> \verbatim
268*>          VT is DOUBLE PRECISION array, dimension (LDVT,NMAX)
269*> \endverbatim
270*>
271*> \param[in] LDVT
272*> \verbatim
273*>          LDVT is INTEGER
274*>          The leading dimension of the array VT.  LDVT >= max(1,NMAX).
275*> \endverbatim
276*>
277*> \param[out] ASAV
278*> \verbatim
279*>          ASAV is DOUBLE PRECISION array, dimension (LDA,NMAX)
280*> \endverbatim
281*>
282*> \param[out] USAV
283*> \verbatim
284*>          USAV is DOUBLE PRECISION array, dimension (LDU,MMAX)
285*> \endverbatim
286*>
287*> \param[out] VTSAV
288*> \verbatim
289*>          VTSAV is DOUBLE PRECISION array, dimension (LDVT,NMAX)
290*> \endverbatim
291*>
292*> \param[out] S
293*> \verbatim
294*>          S is DOUBLE PRECISION array, dimension
295*>                      (max(min(MM,NN)))
296*> \endverbatim
297*>
298*> \param[out] SSAV
299*> \verbatim
300*>          SSAV is DOUBLE PRECISION array, dimension
301*>                      (max(min(MM,NN)))
302*> \endverbatim
303*>
304*> \param[out] E
305*> \verbatim
306*>          E is DOUBLE PRECISION array, dimension
307*>                      (max(min(MM,NN)))
308*> \endverbatim
309*>
310*> \param[out] WORK
311*> \verbatim
312*>          WORK is DOUBLE PRECISION array, dimension (LWORK)
313*> \endverbatim
314*>
315*> \param[in] LWORK
316*> \verbatim
317*>          LWORK is INTEGER
318*>          The number of entries in WORK.  This must be at least
319*>          max(3*MN+MX,5*MN-4)+2*MN**2 for all pairs
320*>          pairs  (MN,MX)=( min(MM(j),NN(j), max(MM(j),NN(j)) )
321*> \endverbatim
322*>
323*> \param[out] IWORK
324*> \verbatim
325*>          IWORK is INTEGER array, dimension at least 8*min(M,N)
326*> \endverbatim
327*>
328*> \param[in] NOUT
329*> \verbatim
330*>          NOUT is INTEGER
331*>          The FORTRAN unit number for printing out error messages
332*>          (e.g., if a routine returns IINFO not equal to 0.)
333*> \endverbatim
334*>
335*> \param[out] INFO
336*> \verbatim
337*>          INFO is INTEGER
338*>          If 0, then everything ran OK.
339*>           -1: NSIZES < 0
340*>           -2: Some MM(j) < 0
341*>           -3: Some NN(j) < 0
342*>           -4: NTYPES < 0
343*>           -7: THRESH < 0
344*>          -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
345*>          -12: LDU < 1 or LDU < MMAX.
346*>          -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ).
347*>          -21: LWORK too small.
348*>          If  DLATMS, or DGESVD returns an error code, the
349*>              absolute value of it is returned.
350*> \endverbatim
351*
352*  Authors:
353*  ========
354*
355*> \author Univ. of Tennessee
356*> \author Univ. of California Berkeley
357*> \author Univ. of Colorado Denver
358*> \author NAG Ltd.
359*
360*> \date June 2016
361*
362*> \ingroup double_eig
363*
364*  =====================================================================
365      SUBROUTINE DDRVBD( NSIZES, MM, NN, NTYPES, DOTYPE, ISEED, THRESH,
366     $                   A, LDA, U, LDU, VT, LDVT, ASAV, USAV, VTSAV, S,
367     $                   SSAV, E, WORK, LWORK, IWORK, NOUT, INFO )
368*
369      IMPLICIT NONE
370*
371*  -- LAPACK test routine (version 3.7.0) --
372*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
373*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
374*     June 2016
375*
376*     .. Scalar Arguments ..
377      INTEGER            INFO, LDA, LDU, LDVT, LWORK, NOUT, NSIZES,
378     $                   NTYPES
379      DOUBLE PRECISION   THRESH
380*     ..
381*     .. Array Arguments ..
382      LOGICAL            DOTYPE( * )
383      INTEGER            ISEED( 4 ), IWORK( * ), MM( * ), NN( * )
384      DOUBLE PRECISION   A( LDA, * ), ASAV( LDA, * ), E( * ), S( * ),
385     $                   SSAV( * ), U( LDU, * ), USAV( LDU, * ),
386     $                   VT( LDVT, * ), VTSAV( LDVT, * ), WORK( * )
387*     ..
388*
389*  =====================================================================
390*
391*     .. Parameters ..
392      DOUBLE PRECISION  ZERO, ONE, TWO, HALF
393      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
394     $                   HALF = 0.5D0 )
395      INTEGER            MAXTYP
396      PARAMETER          ( MAXTYP = 5 )
397*     ..
398*     .. Local Scalars ..
399      LOGICAL            BADMM, BADNN
400      CHARACTER          JOBQ, JOBU, JOBVT, RANGE
401      CHARACTER*3        PATH
402      INTEGER            I, IINFO, IJQ, IJU, IJVT, IL,IU, IWS, IWTMP,
403     $                   ITEMP, J, JSIZE, JTYPE, LSWORK, M, MINWRK,
404     $                   MMAX, MNMAX, MNMIN, MTYPES, N, NFAIL,
405     $                   NMAX, NS, NSI, NSV, NTEST
406      DOUBLE PRECISION   ANORM, DIF, DIV, OVFL, RTUNFL, ULP,
407     $                   ULPINV, UNFL, VL, VU
408*     ..
409*     .. Local Scalars for DGESVDQ ..
410      INTEGER            LIWORK, LRWORK, NUMRANK
411*     ..
412*     .. Local Arrays for DGESVDQ ..
413      DOUBLE PRECISION   RWORK( 2 )
414*     ..
415*     .. Local Arrays ..
416      CHARACTER          CJOB( 4 ), CJOBR( 3 ), CJOBV( 2 )
417      INTEGER            IOLDSD( 4 ), ISEED2( 4 )
418      DOUBLE PRECISION   RESULT( 39 )
419*     ..
420*     .. External Functions ..
421      DOUBLE PRECISION   DLAMCH, DLARND
422      EXTERNAL           DLAMCH, DLARND
423*     ..
424*     .. External Subroutines ..
425      EXTERNAL           ALASVM, DBDT01, DGEJSV, DGESDD, DGESVD,
426     $                   DGESVDQ, DGESVDX, DGESVJ, DLABAD, DLACPY,
427     $                   DLASET, DLATMS, DORT01, DORT03, XERBLA
428*     ..
429*     .. Intrinsic Functions ..
430      INTRINSIC          ABS, DBLE, INT, MAX, MIN
431*     ..
432*     .. Scalars in Common ..
433      LOGICAL            LERR, OK
434      CHARACTER*32       SRNAMT
435      INTEGER            INFOT, NUNIT
436*     ..
437*     .. Common blocks ..
438      COMMON             / INFOC / INFOT, NUNIT, OK, LERR
439      COMMON             / SRNAMC / SRNAMT
440*     ..
441*     .. Data statements ..
442      DATA               CJOB / 'N', 'O', 'S', 'A' /
443      DATA               CJOBR / 'A', 'V', 'I' /
444      DATA               CJOBV / 'N', 'V' /
445*     ..
446*     .. Executable Statements ..
447*
448*     Check for errors
449*
450      INFO = 0
451      BADMM = .FALSE.
452      BADNN = .FALSE.
453      MMAX = 1
454      NMAX = 1
455      MNMAX = 1
456      MINWRK = 1
457      DO 10 J = 1, NSIZES
458         MMAX = MAX( MMAX, MM( J ) )
459         IF( MM( J ).LT.0 )
460     $      BADMM = .TRUE.
461         NMAX = MAX( NMAX, NN( J ) )
462         IF( NN( J ).LT.0 )
463     $      BADNN = .TRUE.
464         MNMAX = MAX( MNMAX, MIN( MM( J ), NN( J ) ) )
465         MINWRK = MAX( MINWRK, MAX( 3*MIN( MM( J ),
466     $            NN( J ) )+MAX( MM( J ), NN( J ) ), 5*MIN( MM( J ),
467     $            NN( J )-4 ) )+2*MIN( MM( J ), NN( J ) )**2 )
468   10 CONTINUE
469*
470*     Check for errors
471*
472      IF( NSIZES.LT.0 ) THEN
473         INFO = -1
474      ELSE IF( BADMM ) THEN
475         INFO = -2
476      ELSE IF( BADNN ) THEN
477         INFO = -3
478      ELSE IF( NTYPES.LT.0 ) THEN
479         INFO = -4
480      ELSE IF( LDA.LT.MAX( 1, MMAX ) ) THEN
481         INFO = -10
482      ELSE IF( LDU.LT.MAX( 1, MMAX ) ) THEN
483         INFO = -12
484      ELSE IF( LDVT.LT.MAX( 1, NMAX ) ) THEN
485         INFO = -14
486      ELSE IF( MINWRK.GT.LWORK ) THEN
487         INFO = -21
488      END IF
489*
490      IF( INFO.NE.0 ) THEN
491         CALL XERBLA( 'DDRVBD', -INFO )
492         RETURN
493      END IF
494*
495*     Initialize constants
496*
497      PATH( 1: 1 ) = 'Double precision'
498      PATH( 2: 3 ) = 'BD'
499      NFAIL = 0
500      NTEST = 0
501      UNFL = DLAMCH( 'Safe minimum' )
502      OVFL = ONE / UNFL
503      CALL DLABAD( UNFL, OVFL )
504      ULP = DLAMCH( 'Precision' )
505      RTUNFL = SQRT( UNFL )
506      ULPINV = ONE / ULP
507      INFOT = 0
508*
509*     Loop over sizes, types
510*
511      DO 240 JSIZE = 1, NSIZES
512         M = MM( JSIZE )
513         N = NN( JSIZE )
514         MNMIN = MIN( M, N )
515*
516         IF( NSIZES.NE.1 ) THEN
517            MTYPES = MIN( MAXTYP, NTYPES )
518         ELSE
519            MTYPES = MIN( MAXTYP+1, NTYPES )
520         END IF
521*
522         DO 230 JTYPE = 1, MTYPES
523            IF( .NOT.DOTYPE( JTYPE ) )
524     $         GO TO 230
525*
526            DO 20 J = 1, 4
527               IOLDSD( J ) = ISEED( J )
528   20       CONTINUE
529*
530*           Compute "A"
531*
532            IF( MTYPES.GT.MAXTYP )
533     $         GO TO 30
534*
535            IF( JTYPE.EQ.1 ) THEN
536*
537*              Zero matrix
538*
539               CALL DLASET( 'Full', M, N, ZERO, ZERO, A, LDA )
540*
541            ELSE IF( JTYPE.EQ.2 ) THEN
542*
543*              Identity matrix
544*
545               CALL DLASET( 'Full', M, N, ZERO, ONE, A, LDA )
546*
547            ELSE
548*
549*              (Scaled) random matrix
550*
551               IF( JTYPE.EQ.3 )
552     $            ANORM = ONE
553               IF( JTYPE.EQ.4 )
554     $            ANORM = UNFL / ULP
555               IF( JTYPE.EQ.5 )
556     $            ANORM = OVFL*ULP
557               CALL DLATMS( M, N, 'U', ISEED, 'N', S, 4, DBLE( MNMIN ),
558     $                      ANORM, M-1, N-1, 'N', A, LDA, WORK, IINFO )
559               IF( IINFO.NE.0 ) THEN
560                  WRITE( NOUT, FMT = 9996 )'Generator', IINFO, M, N,
561     $               JTYPE, IOLDSD
562                  INFO = ABS( IINFO )
563                  RETURN
564               END IF
565            END IF
566*
567   30       CONTINUE
568            CALL DLACPY( 'F', M, N, A, LDA, ASAV, LDA )
569*
570*           Do for minimal and adequate (for blocking) workspace
571*
572            DO 220 IWS = 1, 4
573*
574               DO 40 J = 1, 32
575                  RESULT( J ) = -ONE
576   40          CONTINUE
577*
578*              Test DGESVD: Factorize A
579*
580               IWTMP = MAX( 3*MIN( M, N )+MAX( M, N ), 5*MIN( M, N ) )
581               LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
582               LSWORK = MIN( LSWORK, LWORK )
583               LSWORK = MAX( LSWORK, 1 )
584               IF( IWS.EQ.4 )
585     $            LSWORK = LWORK
586*
587               IF( IWS.GT.1 )
588     $            CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
589               SRNAMT = 'DGESVD'
590               CALL DGESVD( 'A', 'A', M, N, A, LDA, SSAV, USAV, LDU,
591     $                      VTSAV, LDVT, WORK, LSWORK, IINFO )
592               IF( IINFO.NE.0 ) THEN
593                  WRITE( NOUT, FMT = 9995 )'GESVD', IINFO, M, N, JTYPE,
594     $               LSWORK, IOLDSD
595                  INFO = ABS( IINFO )
596                  RETURN
597               END IF
598*
599*              Do tests 1--4
600*
601               CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
602     $                      VTSAV, LDVT, WORK, RESULT( 1 ) )
603               IF( M.NE.0 .AND. N.NE.0 ) THEN
604                  CALL DORT01( 'Columns', M, M, USAV, LDU, WORK, LWORK,
605     $                         RESULT( 2 ) )
606                  CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK, LWORK,
607     $                         RESULT( 3 ) )
608               END IF
609               RESULT( 4 ) = ZERO
610               DO 50 I = 1, MNMIN - 1
611                  IF( SSAV( I ).LT.SSAV( I+1 ) )
612     $               RESULT( 4 ) = ULPINV
613                  IF( SSAV( I ).LT.ZERO )
614     $               RESULT( 4 ) = ULPINV
615   50          CONTINUE
616               IF( MNMIN.GE.1 ) THEN
617                  IF( SSAV( MNMIN ).LT.ZERO )
618     $               RESULT( 4 ) = ULPINV
619               END IF
620*
621*              Do partial SVDs, comparing to SSAV, USAV, and VTSAV
622*
623               RESULT( 5 ) = ZERO
624               RESULT( 6 ) = ZERO
625               RESULT( 7 ) = ZERO
626               DO 80 IJU = 0, 3
627                  DO 70 IJVT = 0, 3
628                     IF( ( IJU.EQ.3 .AND. IJVT.EQ.3 ) .OR.
629     $                   ( IJU.EQ.1 .AND. IJVT.EQ.1 ) )GO TO 70
630                     JOBU = CJOB( IJU+1 )
631                     JOBVT = CJOB( IJVT+1 )
632                     CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
633                     SRNAMT = 'DGESVD'
634                     CALL DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
635     $                            VT, LDVT, WORK, LSWORK, IINFO )
636*
637*                    Compare U
638*
639                     DIF = ZERO
640                     IF( M.GT.0 .AND. N.GT.0 ) THEN
641                        IF( IJU.EQ.1 ) THEN
642                           CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
643     $                                  LDU, A, LDA, WORK, LWORK, DIF,
644     $                                  IINFO )
645                        ELSE IF( IJU.EQ.2 ) THEN
646                           CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
647     $                                  LDU, U, LDU, WORK, LWORK, DIF,
648     $                                  IINFO )
649                        ELSE IF( IJU.EQ.3 ) THEN
650                           CALL DORT03( 'C', M, M, M, MNMIN, USAV, LDU,
651     $                                  U, LDU, WORK, LWORK, DIF,
652     $                                  IINFO )
653                        END IF
654                     END IF
655                     RESULT( 5 ) = MAX( RESULT( 5 ), DIF )
656*
657*                    Compare VT
658*
659                     DIF = ZERO
660                     IF( M.GT.0 .AND. N.GT.0 ) THEN
661                        IF( IJVT.EQ.1 ) THEN
662                           CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
663     $                                  LDVT, A, LDA, WORK, LWORK, DIF,
664     $                                  IINFO )
665                        ELSE IF( IJVT.EQ.2 ) THEN
666                           CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
667     $                                  LDVT, VT, LDVT, WORK, LWORK,
668     $                                  DIF, IINFO )
669                        ELSE IF( IJVT.EQ.3 ) THEN
670                           CALL DORT03( 'R', N, N, N, MNMIN, VTSAV,
671     $                                  LDVT, VT, LDVT, WORK, LWORK,
672     $                                  DIF, IINFO )
673                        END IF
674                     END IF
675                     RESULT( 6 ) = MAX( RESULT( 6 ), DIF )
676*
677*                    Compare S
678*
679                     DIF = ZERO
680                     DIV = MAX( MNMIN*ULP*S( 1 ), UNFL )
681                     DO 60 I = 1, MNMIN - 1
682                        IF( SSAV( I ).LT.SSAV( I+1 ) )
683     $                     DIF = ULPINV
684                        IF( SSAV( I ).LT.ZERO )
685     $                     DIF = ULPINV
686                        DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
687   60                CONTINUE
688                     RESULT( 7 ) = MAX( RESULT( 7 ), DIF )
689   70             CONTINUE
690   80          CONTINUE
691*
692*              Test DGESDD: Factorize A
693*
694               IWTMP = 5*MNMIN*MNMIN + 9*MNMIN + MAX( M, N )
695               LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
696               LSWORK = MIN( LSWORK, LWORK )
697               LSWORK = MAX( LSWORK, 1 )
698               IF( IWS.EQ.4 )
699     $            LSWORK = LWORK
700*
701               CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
702               SRNAMT = 'DGESDD'
703               CALL DGESDD( 'A', M, N, A, LDA, SSAV, USAV, LDU, VTSAV,
704     $                      LDVT, WORK, LSWORK, IWORK, IINFO )
705               IF( IINFO.NE.0 ) THEN
706                  WRITE( NOUT, FMT = 9995 )'GESDD', IINFO, M, N, JTYPE,
707     $               LSWORK, IOLDSD
708                  INFO = ABS( IINFO )
709                  RETURN
710               END IF
711*
712*              Do tests 8--11
713*
714               CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
715     $                      VTSAV, LDVT, WORK, RESULT( 8 ) )
716               IF( M.NE.0 .AND. N.NE.0 ) THEN
717                  CALL DORT01( 'Columns', M, M, USAV, LDU, WORK, LWORK,
718     $                         RESULT( 9 ) )
719                  CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK, LWORK,
720     $                         RESULT( 10 ) )
721               END IF
722               RESULT( 11 ) = ZERO
723               DO 90 I = 1, MNMIN - 1
724                  IF( SSAV( I ).LT.SSAV( I+1 ) )
725     $               RESULT( 11 ) = ULPINV
726                  IF( SSAV( I ).LT.ZERO )
727     $               RESULT( 11 ) = ULPINV
728   90          CONTINUE
729               IF( MNMIN.GE.1 ) THEN
730                  IF( SSAV( MNMIN ).LT.ZERO )
731     $               RESULT( 11 ) = ULPINV
732               END IF
733*
734*              Do partial SVDs, comparing to SSAV, USAV, and VTSAV
735*
736               RESULT( 12 ) = ZERO
737               RESULT( 13 ) = ZERO
738               RESULT( 14 ) = ZERO
739               DO 110 IJQ = 0, 2
740                  JOBQ = CJOB( IJQ+1 )
741                  CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
742                  SRNAMT = 'DGESDD'
743                  CALL DGESDD( JOBQ, M, N, A, LDA, S, U, LDU, VT, LDVT,
744     $                         WORK, LSWORK, IWORK, IINFO )
745*
746*                 Compare U
747*
748                  DIF = ZERO
749                  IF( M.GT.0 .AND. N.GT.0 ) THEN
750                     IF( IJQ.EQ.1 ) THEN
751                        IF( M.GE.N ) THEN
752                           CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
753     $                                  LDU, A, LDA, WORK, LWORK, DIF,
754     $                                  INFO )
755                        ELSE
756                           CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
757     $                                  LDU, U, LDU, WORK, LWORK, DIF,
758     $                                  INFO )
759                        END IF
760                     ELSE IF( IJQ.EQ.2 ) THEN
761                        CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV, LDU,
762     $                               U, LDU, WORK, LWORK, DIF, INFO )
763                     END IF
764                  END IF
765                  RESULT( 12 ) = MAX( RESULT( 12 ), DIF )
766*
767*                 Compare VT
768*
769                  DIF = ZERO
770                  IF( M.GT.0 .AND. N.GT.0 ) THEN
771                     IF( IJQ.EQ.1 ) THEN
772                        IF( M.GE.N ) THEN
773                           CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
774     $                                  LDVT, VT, LDVT, WORK, LWORK,
775     $                                  DIF, INFO )
776                        ELSE
777                           CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
778     $                                  LDVT, A, LDA, WORK, LWORK, DIF,
779     $                                  INFO )
780                        END IF
781                     ELSE IF( IJQ.EQ.2 ) THEN
782                        CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
783     $                               LDVT, VT, LDVT, WORK, LWORK, DIF,
784     $                               INFO )
785                     END IF
786                  END IF
787                  RESULT( 13 ) = MAX( RESULT( 13 ), DIF )
788*
789*                 Compare S
790*
791                  DIF = ZERO
792                  DIV = MAX( MNMIN*ULP*S( 1 ), UNFL )
793                  DO 100 I = 1, MNMIN - 1
794                     IF( SSAV( I ).LT.SSAV( I+1 ) )
795     $                  DIF = ULPINV
796                     IF( SSAV( I ).LT.ZERO )
797     $                  DIF = ULPINV
798                     DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
799  100             CONTINUE
800                  RESULT( 14 ) = MAX( RESULT( 14 ), DIF )
801  110          CONTINUE
802*
803*              Test DGESVDQ
804*              Note: DGESVDQ only works for M >= N
805*
806               RESULT( 36 ) = ZERO
807               RESULT( 37 ) = ZERO
808               RESULT( 38 ) = ZERO
809               RESULT( 39 ) = ZERO
810*
811               IF( M.GE.N ) THEN
812                  IWTMP = 5*MNMIN*MNMIN + 9*MNMIN + MAX( M, N )
813                  LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
814                  LSWORK = MIN( LSWORK, LWORK )
815                  LSWORK = MAX( LSWORK, 1 )
816                  IF( IWS.EQ.4 )
817     $               LSWORK = LWORK
818*
819                  CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
820                  SRNAMT = 'DGESVDQ'
821*
822                  LRWORK = 2
823                  LIWORK = MAX( N, 1 )
824                  CALL DGESVDQ( 'H', 'N', 'N', 'A', 'A',
825     $                          M, N, A, LDA, SSAV, USAV, LDU,
826     $                          VTSAV, LDVT, NUMRANK, IWORK, LIWORK,
827     $                          WORK, LWORK, RWORK, LRWORK, IINFO )
828*
829                  IF( IINFO.NE.0 ) THEN
830                     WRITE( NOUT, FMT = 9995 )'DGESVDQ', IINFO, M, N,
831     $               JTYPE, LSWORK, IOLDSD
832                     INFO = ABS( IINFO )
833                     RETURN
834                  END IF
835*
836*                 Do tests 36--39
837*
838                  CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
839     $                         VTSAV, LDVT, WORK, RESULT( 36 ) )
840                  IF( M.NE.0 .AND. N.NE.0 ) THEN
841                     CALL DORT01( 'Columns', M, M, USAV, LDU, WORK,
842     $                            LWORK, RESULT( 37 ) )
843                     CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK,
844     $                            LWORK, RESULT( 38 ) )
845                  END IF
846                  RESULT( 39 ) = ZERO
847                  DO 199 I = 1, MNMIN - 1
848                     IF( SSAV( I ).LT.SSAV( I+1 ) )
849     $                  RESULT( 39 ) = ULPINV
850                     IF( SSAV( I ).LT.ZERO )
851     $                  RESULT( 39 ) = ULPINV
852  199             CONTINUE
853                  IF( MNMIN.GE.1 ) THEN
854                     IF( SSAV( MNMIN ).LT.ZERO )
855     $                  RESULT( 39 ) = ULPINV
856                  END IF
857               END IF
858*
859*              Test DGESVJ
860*              Note: DGESVJ only works for M >= N
861*
862               RESULT( 15 ) = ZERO
863               RESULT( 16 ) = ZERO
864               RESULT( 17 ) = ZERO
865               RESULT( 18 ) = ZERO
866*
867               IF( M.GE.N ) THEN
868                  IWTMP = 5*MNMIN*MNMIN + 9*MNMIN + MAX( M, N )
869                  LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
870                  LSWORK = MIN( LSWORK, LWORK )
871                  LSWORK = MAX( LSWORK, 1 )
872                  IF( IWS.EQ.4 )
873     $               LSWORK = LWORK
874*
875                  CALL DLACPY( 'F', M, N, ASAV, LDA, USAV, LDA )
876                  SRNAMT = 'DGESVJ'
877                  CALL DGESVJ( 'G', 'U', 'V', M, N, USAV, LDA, SSAV,
878     &                        0, A, LDVT, WORK, LWORK, INFO )
879*
880*                 DGESVJ returns V not VT
881*
882                  DO J=1,N
883                     DO I=1,N
884                        VTSAV(J,I) = A(I,J)
885                     END DO
886                  END DO
887*
888                  IF( IINFO.NE.0 ) THEN
889                     WRITE( NOUT, FMT = 9995 )'GESVJ', IINFO, M, N,
890     $               JTYPE, LSWORK, IOLDSD
891                     INFO = ABS( IINFO )
892                     RETURN
893                  END IF
894*
895*                 Do tests 15--18
896*
897                  CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
898     $                         VTSAV, LDVT, WORK, RESULT( 15 ) )
899                  IF( M.NE.0 .AND. N.NE.0 ) THEN
900                     CALL DORT01( 'Columns', M, M, USAV, LDU, WORK,
901     $                            LWORK, RESULT( 16 ) )
902                     CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK,
903     $                            LWORK, RESULT( 17 ) )
904                  END IF
905                  RESULT( 18 ) = ZERO
906                  DO 120 I = 1, MNMIN - 1
907                     IF( SSAV( I ).LT.SSAV( I+1 ) )
908     $                  RESULT( 18 ) = ULPINV
909                     IF( SSAV( I ).LT.ZERO )
910     $                  RESULT( 18 ) = ULPINV
911  120             CONTINUE
912                  IF( MNMIN.GE.1 ) THEN
913                     IF( SSAV( MNMIN ).LT.ZERO )
914     $                  RESULT( 18 ) = ULPINV
915                  END IF
916               END IF
917*
918*              Test DGEJSV
919*              Note: DGEJSV only works for M >= N
920*
921               RESULT( 19 ) = ZERO
922               RESULT( 20 ) = ZERO
923               RESULT( 21 ) = ZERO
924               RESULT( 22 ) = ZERO
925               IF( M.GE.N ) THEN
926                  IWTMP = 5*MNMIN*MNMIN + 9*MNMIN + MAX( M, N )
927                  LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
928                  LSWORK = MIN( LSWORK, LWORK )
929                  LSWORK = MAX( LSWORK, 1 )
930                  IF( IWS.EQ.4 )
931     $               LSWORK = LWORK
932*
933                  CALL DLACPY( 'F', M, N, ASAV, LDA, VTSAV, LDA )
934                  SRNAMT = 'DGEJSV'
935                  CALL DGEJSV( 'G', 'U', 'V', 'R', 'N', 'N',
936     &                   M, N, VTSAV, LDA, SSAV, USAV, LDU, A, LDVT,
937     &                   WORK, LWORK, IWORK, INFO )
938*
939*                 DGEJSV returns V not VT
940*
941                  DO 140 J=1,N
942                     DO 130 I=1,N
943                        VTSAV(J,I) = A(I,J)
944  130                END DO
945  140             END DO
946*
947                  IF( IINFO.NE.0 ) THEN
948                     WRITE( NOUT, FMT = 9995 )'GEJSV', IINFO, M, N,
949     $               JTYPE, LSWORK, IOLDSD
950                     INFO = ABS( IINFO )
951                     RETURN
952                  END IF
953*
954*                 Do tests 19--22
955*
956                  CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
957     $                         VTSAV, LDVT, WORK, RESULT( 19 ) )
958                  IF( M.NE.0 .AND. N.NE.0 ) THEN
959                     CALL DORT01( 'Columns', M, M, USAV, LDU, WORK,
960     $                            LWORK, RESULT( 20 ) )
961                     CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK,
962     $                            LWORK, RESULT( 21 ) )
963                  END IF
964                  RESULT( 22 ) = ZERO
965                  DO 150 I = 1, MNMIN - 1
966                     IF( SSAV( I ).LT.SSAV( I+1 ) )
967     $                  RESULT( 22 ) = ULPINV
968                     IF( SSAV( I ).LT.ZERO )
969     $                  RESULT( 22 ) = ULPINV
970  150             CONTINUE
971                  IF( MNMIN.GE.1 ) THEN
972                     IF( SSAV( MNMIN ).LT.ZERO )
973     $                  RESULT( 22 ) = ULPINV
974                  END IF
975               END IF
976*
977*              Test DGESVDX
978*
979               CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
980               CALL DGESVDX( 'V', 'V', 'A', M, N, A, LDA,
981     $                       VL, VU, IL, IU, NS, SSAV, USAV, LDU,
982     $                       VTSAV, LDVT, WORK, LWORK, IWORK,
983     $                       IINFO )
984               IF( IINFO.NE.0 ) THEN
985                  WRITE( NOUT, FMT = 9995 )'GESVDX', IINFO, M, N,
986     $               JTYPE, LSWORK, IOLDSD
987                  INFO = ABS( IINFO )
988                  RETURN
989               END IF
990*
991*              Do tests 23--29
992*
993               RESULT( 23 ) = ZERO
994               RESULT( 24 ) = ZERO
995               RESULT( 25 ) = ZERO
996               CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
997     $                      VTSAV, LDVT, WORK, RESULT( 23 ) )
998               IF( M.NE.0 .AND. N.NE.0 ) THEN
999                  CALL DORT01( 'Columns', M, M, USAV, LDU, WORK, LWORK,
1000     $                         RESULT( 24 ) )
1001                  CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK, LWORK,
1002     $                         RESULT( 25 ) )
1003               END IF
1004               RESULT( 26 ) = ZERO
1005               DO 160 I = 1, MNMIN - 1
1006                  IF( SSAV( I ).LT.SSAV( I+1 ) )
1007     $               RESULT( 26 ) = ULPINV
1008                  IF( SSAV( I ).LT.ZERO )
1009     $               RESULT( 26 ) = ULPINV
1010  160          CONTINUE
1011               IF( MNMIN.GE.1 ) THEN
1012                  IF( SSAV( MNMIN ).LT.ZERO )
1013     $               RESULT( 26 ) = ULPINV
1014               END IF
1015*
1016*              Do partial SVDs, comparing to SSAV, USAV, and VTSAV
1017*
1018               RESULT( 27 ) = ZERO
1019               RESULT( 28 ) = ZERO
1020               RESULT( 29 ) = ZERO
1021               DO 180 IJU = 0, 1
1022                  DO 170 IJVT = 0, 1
1023                     IF( ( IJU.EQ.0 .AND. IJVT.EQ.0 ) .OR.
1024     $                   ( IJU.EQ.1 .AND. IJVT.EQ.1 ) )GO TO 170
1025                     JOBU = CJOBV( IJU+1 )
1026                     JOBVT = CJOBV( IJVT+1 )
1027                     RANGE = CJOBR( 1 )
1028                     CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
1029                     CALL DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA,
1030     $                             VL, VU, IL, IU, NS, S, U, LDU,
1031     $                             VT, LDVT, WORK, LWORK, IWORK,
1032     $                             IINFO )
1033*
1034*                    Compare U
1035*
1036                     DIF = ZERO
1037                     IF( M.GT.0 .AND. N.GT.0 ) THEN
1038                        IF( IJU.EQ.1 ) THEN
1039                           CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
1040     $                                  LDU, U, LDU, WORK, LWORK, DIF,
1041     $                                  IINFO )
1042                        END IF
1043                     END IF
1044                     RESULT( 27 ) = MAX( RESULT( 27 ), DIF )
1045*
1046*                    Compare VT
1047*
1048                     DIF = ZERO
1049                     IF( M.GT.0 .AND. N.GT.0 ) THEN
1050                        IF( IJVT.EQ.1 ) THEN
1051                           CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
1052     $                                  LDVT, VT, LDVT, WORK, LWORK,
1053     $                                  DIF, IINFO )
1054                        END IF
1055                     END IF
1056                     RESULT( 28 ) = MAX( RESULT( 28 ), DIF )
1057*
1058*                    Compare S
1059*
1060                     DIF = ZERO
1061                     DIV = MAX( MNMIN*ULP*S( 1 ), UNFL )
1062                     DO 190 I = 1, MNMIN - 1
1063                        IF( SSAV( I ).LT.SSAV( I+1 ) )
1064     $                     DIF = ULPINV
1065                        IF( SSAV( I ).LT.ZERO )
1066     $                     DIF = ULPINV
1067                        DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
1068  190                CONTINUE
1069                     RESULT( 29 ) = MAX( RESULT( 29 ), DIF )
1070  170             CONTINUE
1071  180          CONTINUE
1072*
1073*              Do tests 30--32: DGESVDX( 'V', 'V', 'I' )
1074*
1075               DO 200 I = 1, 4
1076                  ISEED2( I ) = ISEED( I )
1077  200          CONTINUE
1078               IF( MNMIN.LE.1 ) THEN
1079                  IL = 1
1080                  IU = MAX( 1, MNMIN )
1081               ELSE
1082                  IL = 1 + INT( ( MNMIN-1 )*DLARND( 1, ISEED2 ) )
1083                  IU = 1 + INT( ( MNMIN-1 )*DLARND( 1, ISEED2 ) )
1084                  IF( IU.LT.IL ) THEN
1085                     ITEMP = IU
1086                     IU = IL
1087                     IL = ITEMP
1088                  END IF
1089               END IF
1090               CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
1091               CALL DGESVDX( 'V', 'V', 'I', M, N, A, LDA,
1092     $                       VL, VU, IL, IU, NSI, S, U, LDU,
1093     $                       VT, LDVT, WORK, LWORK, IWORK,
1094     $                       IINFO )
1095               IF( IINFO.NE.0 ) THEN
1096                  WRITE( NOUT, FMT = 9995 )'GESVDX', IINFO, M, N,
1097     $               JTYPE, LSWORK, IOLDSD
1098                  INFO = ABS( IINFO )
1099                  RETURN
1100               END IF
1101*
1102               RESULT( 30 ) = ZERO
1103               RESULT( 31 ) = ZERO
1104               RESULT( 32 ) = ZERO
1105               CALL DBDT05( M, N, ASAV, LDA, S, NSI, U, LDU,
1106     $                      VT, LDVT, WORK, RESULT( 30 ) )
1107               CALL DORT01( 'Columns', M, NSI, U, LDU, WORK, LWORK,
1108     $                      RESULT( 31 ) )
1109               CALL DORT01( 'Rows', NSI, N, VT, LDVT, WORK, LWORK,
1110     $                      RESULT( 32 ) )
1111*
1112*              Do tests 33--35: DGESVDX( 'V', 'V', 'V' )
1113*
1114               IF( MNMIN.GT.0 .AND. NSI.GT.1 ) THEN
1115                  IF( IL.NE.1 ) THEN
1116                     VU = SSAV( IL ) +
1117     $                    MAX( HALF*ABS( SSAV( IL )-SSAV( IL-1 ) ),
1118     $                    ULP*ANORM, TWO*RTUNFL )
1119                  ELSE
1120                     VU = SSAV( 1 ) +
1121     $                    MAX( HALF*ABS( SSAV( NS )-SSAV( 1 ) ),
1122     $                    ULP*ANORM, TWO*RTUNFL )
1123                  END IF
1124                  IF( IU.NE.NS ) THEN
1125                     VL = SSAV( IU ) - MAX( ULP*ANORM, TWO*RTUNFL,
1126     $                    HALF*ABS( SSAV( IU+1 )-SSAV( IU ) ) )
1127                  ELSE
1128                     VL = SSAV( NS ) - MAX( ULP*ANORM, TWO*RTUNFL,
1129     $                    HALF*ABS( SSAV( NS )-SSAV( 1 ) ) )
1130                  END IF
1131                  VL = MAX( VL,ZERO )
1132                  VU = MAX( VU,ZERO )
1133                  IF( VL.GE.VU ) VU = MAX( VU*2, VU+VL+HALF )
1134               ELSE
1135                  VL = ZERO
1136                  VU = ONE
1137               END IF
1138               CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
1139               CALL DGESVDX( 'V', 'V', 'V', M, N, A, LDA,
1140     $                       VL, VU, IL, IU, NSV, S, U, LDU,
1141     $                       VT, LDVT, WORK, LWORK, IWORK,
1142     $                       IINFO )
1143               IF( IINFO.NE.0 ) THEN
1144                  WRITE( NOUT, FMT = 9995 )'GESVDX', IINFO, M, N,
1145     $               JTYPE, LSWORK, IOLDSD
1146                  INFO = ABS( IINFO )
1147                  RETURN
1148               END IF
1149*
1150               RESULT( 33 ) = ZERO
1151               RESULT( 34 ) = ZERO
1152               RESULT( 35 ) = ZERO
1153               CALL DBDT05( M, N, ASAV, LDA, S, NSV, U, LDU,
1154     $                      VT, LDVT, WORK, RESULT( 33 ) )
1155               CALL DORT01( 'Columns', M, NSV, U, LDU, WORK, LWORK,
1156     $                      RESULT( 34 ) )
1157               CALL DORT01( 'Rows', NSV, N, VT, LDVT, WORK, LWORK,
1158     $                      RESULT( 35 ) )
1159*
1160*              End of Loop -- Check for RESULT(j) > THRESH
1161*
1162               DO 210 J = 1, 39
1163                  IF( RESULT( J ).GE.THRESH ) THEN
1164                     IF( NFAIL.EQ.0 ) THEN
1165                        WRITE( NOUT, FMT = 9999 )
1166                        WRITE( NOUT, FMT = 9998 )
1167                     END IF
1168                     WRITE( NOUT, FMT = 9997 )M, N, JTYPE, IWS, IOLDSD,
1169     $                  J, RESULT( J )
1170                     NFAIL = NFAIL + 1
1171                  END IF
1172  210          CONTINUE
1173               NTEST = NTEST + 39
1174  220       CONTINUE
1175  230    CONTINUE
1176  240 CONTINUE
1177*
1178*     Summary
1179*
1180      CALL ALASVM( PATH, NOUT, NFAIL, NTEST, 0 )
1181*
1182 9999 FORMAT( ' SVD -- Real Singular Value Decomposition Driver ',
1183     $      / ' Matrix types (see DDRVBD for details):',
1184     $      / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
1185     $      / ' 3 = Evenly spaced singular values near 1',
1186     $      / ' 4 = Evenly spaced singular values near underflow',
1187     $      / ' 5 = Evenly spaced singular values near overflow', / /
1188     $      ' Tests performed: ( A is dense, U and V are orthogonal,',
1189     $      / 19X, ' S is an array, and Upartial, VTpartial, and',
1190     $      / 19X, ' Spartial are partially computed U, VT and S),', / )
1191 9998 FORMAT( ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1192     $      / ' 2 = | I - U**T U | / ( M ulp ) ',
1193     $      / ' 3 = | I - VT VT**T | / ( N ulp ) ',
1194     $      / ' 4 = 0 if S contains min(M,N) nonnegative values in',
1195     $      ' decreasing order, else 1/ulp',
1196     $      / ' 5 = | U - Upartial | / ( M ulp )',
1197     $      / ' 6 = | VT - VTpartial | / ( N ulp )',
1198     $      / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
1199     $      / ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1200     $      / ' 9 = | I - U**T U | / ( M ulp ) ',
1201     $      / '10 = | I - VT VT**T | / ( N ulp ) ',
1202     $      / '11 = 0 if S contains min(M,N) nonnegative values in',
1203     $      ' decreasing order, else 1/ulp',
1204     $      / '12 = | U - Upartial | / ( M ulp )',
1205     $      / '13 = | VT - VTpartial | / ( N ulp )',
1206     $      / '14 = | S - Spartial | / ( min(M,N) ulp |S| )',
1207     $      / '15 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1208     $      / '16 = | I - U**T U | / ( M ulp ) ',
1209     $      / '17 = | I - VT VT**T | / ( N ulp ) ',
1210     $      / '18 = 0 if S contains min(M,N) nonnegative values in',
1211     $      ' decreasing order, else 1/ulp',
1212     $      / '19 = | U - Upartial | / ( M ulp )',
1213     $      / '20 = | VT - VTpartial | / ( N ulp )',
1214     $      / '21 = | S - Spartial | / ( min(M,N) ulp |S| )',
1215     $      / '22 = 0 if S contains min(M,N) nonnegative values in',
1216     $      ' decreasing order, else 1/ulp',
1217     $      / '23 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ),',
1218     $      ' DGESVDX(V,V,A) ',
1219     $      / '24 = | I - U**T U | / ( M ulp ) ',
1220     $      / '25 = | I - VT VT**T | / ( N ulp ) ',
1221     $      / '26 = 0 if S contains min(M,N) nonnegative values in',
1222     $      ' decreasing order, else 1/ulp',
1223     $      / '27 = | U - Upartial | / ( M ulp )',
1224     $      / '28 = | VT - VTpartial | / ( N ulp )',
1225     $      / '29 = | S - Spartial | / ( min(M,N) ulp |S| )',
1226     $      / '30 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp ),',
1227     $      ' DGESVDX(V,V,I) ',
1228     $      / '31 = | I - U**T U | / ( M ulp ) ',
1229     $      / '32 = | I - VT VT**T | / ( N ulp ) ',
1230     $      / '33 = | U**T A VT**T - diag(S) | / ( |A| max(M,N) ulp ),',
1231     $      ' DGESVDX(V,V,V) ',
1232     $      / '34 = | I - U**T U | / ( M ulp ) ',
1233     $      / '35 = | I - VT VT**T | / ( N ulp ) ',
1234     $      ' DGESVDQ(H,N,N,A,A',
1235     $      / '36 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
1236     $      / '37 = | I - U**T U | / ( M ulp ) ',
1237     $      / '38 = | I - VT VT**T | / ( N ulp ) ',
1238     $      / '39 = 0 if S contains min(M,N) nonnegative values in',
1239     $      ' decreasing order, else 1/ulp',
1240     $      / / )
1241 9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1,
1242     $      ', seed=', 4( I4, ',' ), ' test(', I2, ')=', G11.4 )
1243 9996 FORMAT( ' DDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
1244     $      I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
1245     $      I5, ')' )
1246 9995 FORMAT( ' DDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
1247     $      I6, ', N=', I6, ', JTYPE=', I6, ', LSWORK=', I6, / 9X,
1248     $      'ISEED=(', 3( I5, ',' ), I5, ')' )
1249*
1250      RETURN
1251*
1252*     End of DDRVBD
1253*
1254      END
1255