1*> \brief \b DGESDD
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGESDD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesdd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesdd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesdd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
22*                          WORK, LWORK, IWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          JOBZ
26*       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
27*       ..
28*       .. Array Arguments ..
29*       INTEGER            IWORK( * )
30*       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
31*      $                   VT( LDVT, * ), WORK( * )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> DGESDD computes the singular value decomposition (SVD) of a real
41*> M-by-N matrix A, optionally computing the left and right singular
42*> vectors.  If singular vectors are desired, it uses a
43*> divide-and-conquer algorithm.
44*>
45*> The SVD is written
46*>
47*>      A = U * SIGMA * transpose(V)
48*>
49*> where SIGMA is an M-by-N matrix which is zero except for its
50*> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
51*> V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
52*> are the singular values of A; they are real and non-negative, and
53*> are returned in descending order.  The first min(m,n) columns of
54*> U and V are the left and right singular vectors of A.
55*>
56*> Note that the routine returns VT = V**T, not V.
57*>
58*> The divide and conquer algorithm makes very mild assumptions about
59*> floating point arithmetic. It will work on machines with a guard
60*> digit in add/subtract, or on those binary machines without guard
61*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
62*> Cray-2. It could conceivably fail on hexadecimal or decimal machines
63*> without guard digits, but we know of none.
64*> \endverbatim
65*
66*  Arguments:
67*  ==========
68*
69*> \param[in] JOBZ
70*> \verbatim
71*>          JOBZ is CHARACTER*1
72*>          Specifies options for computing all or part of the matrix U:
73*>          = 'A':  all M columns of U and all N rows of V**T are
74*>                  returned in the arrays U and VT;
75*>          = 'S':  the first min(M,N) columns of U and the first
76*>                  min(M,N) rows of V**T are returned in the arrays U
77*>                  and VT;
78*>          = 'O':  If M >= N, the first N columns of U are overwritten
79*>                  on the array A and all rows of V**T are returned in
80*>                  the array VT;
81*>                  otherwise, all columns of U are returned in the
82*>                  array U and the first M rows of V**T are overwritten
83*>                  in the array A;
84*>          = 'N':  no columns of U or rows of V**T are computed.
85*> \endverbatim
86*>
87*> \param[in] M
88*> \verbatim
89*>          M is INTEGER
90*>          The number of rows of the input matrix A.  M >= 0.
91*> \endverbatim
92*>
93*> \param[in] N
94*> \verbatim
95*>          N is INTEGER
96*>          The number of columns of the input matrix A.  N >= 0.
97*> \endverbatim
98*>
99*> \param[in,out] A
100*> \verbatim
101*>          A is DOUBLE PRECISION array, dimension (LDA,N)
102*>          On entry, the M-by-N matrix A.
103*>          On exit,
104*>          if JOBZ = 'O',  A is overwritten with the first N columns
105*>                          of U (the left singular vectors, stored
106*>                          columnwise) if M >= N;
107*>                          A is overwritten with the first M rows
108*>                          of V**T (the right singular vectors, stored
109*>                          rowwise) otherwise.
110*>          if JOBZ .ne. 'O', the contents of A are destroyed.
111*> \endverbatim
112*>
113*> \param[in] LDA
114*> \verbatim
115*>          LDA is INTEGER
116*>          The leading dimension of the array A.  LDA >= max(1,M).
117*> \endverbatim
118*>
119*> \param[out] S
120*> \verbatim
121*>          S is DOUBLE PRECISION array, dimension (min(M,N))
122*>          The singular values of A, sorted so that S(i) >= S(i+1).
123*> \endverbatim
124*>
125*> \param[out] U
126*> \verbatim
127*>          U is DOUBLE PRECISION array, dimension (LDU,UCOL)
128*>          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
129*>          UCOL = min(M,N) if JOBZ = 'S'.
130*>          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
131*>          orthogonal matrix U;
132*>          if JOBZ = 'S', U contains the first min(M,N) columns of U
133*>          (the left singular vectors, stored columnwise);
134*>          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
135*> \endverbatim
136*>
137*> \param[in] LDU
138*> \verbatim
139*>          LDU is INTEGER
140*>          The leading dimension of the array U.  LDU >= 1; if
141*>          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
142*> \endverbatim
143*>
144*> \param[out] VT
145*> \verbatim
146*>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
147*>          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
148*>          N-by-N orthogonal matrix V**T;
149*>          if JOBZ = 'S', VT contains the first min(M,N) rows of
150*>          V**T (the right singular vectors, stored rowwise);
151*>          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
152*> \endverbatim
153*>
154*> \param[in] LDVT
155*> \verbatim
156*>          LDVT is INTEGER
157*>          The leading dimension of the array VT.  LDVT >= 1;
158*>          if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
159*>          if JOBZ = 'S', LDVT >= min(M,N).
160*> \endverbatim
161*>
162*> \param[out] WORK
163*> \verbatim
164*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
165*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
166*> \endverbatim
167*>
168*> \param[in] LWORK
169*> \verbatim
170*>          LWORK is INTEGER
171*>          The dimension of the array WORK. LWORK >= 1.
172*>          If LWORK = -1, a workspace query is assumed.  The optimal
173*>          size for the WORK array is calculated and stored in WORK(1),
174*>          and no other work except argument checking is performed.
175*>
176*>          Let mx = max(M,N) and mn = min(M,N).
177*>          If JOBZ = 'N', LWORK >= 3*mn + max( mx, 7*mn ).
178*>          If JOBZ = 'O', LWORK >= 3*mn + max( mx, 5*mn*mn + 4*mn ).
179*>          If JOBZ = 'S', LWORK >= 4*mn*mn + 7*mn.
180*>          If JOBZ = 'A', LWORK >= 4*mn*mn + 6*mn + mx.
181*>          These are not tight minimums in all cases; see comments inside code.
182*>          For good performance, LWORK should generally be larger;
183*>          a query is recommended.
184*> \endverbatim
185*>
186*> \param[out] IWORK
187*> \verbatim
188*>          IWORK is INTEGER array, dimension (8*min(M,N))
189*> \endverbatim
190*>
191*> \param[out] INFO
192*> \verbatim
193*>          INFO is INTEGER
194*>          = 0:  successful exit.
195*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
196*>          > 0:  DBDSDC did not converge, updating process failed.
197*> \endverbatim
198*
199*  Authors:
200*  ========
201*
202*> \author Univ. of Tennessee
203*> \author Univ. of California Berkeley
204*> \author Univ. of Colorado Denver
205*> \author NAG Ltd.
206*
207*> \date June 2016
208*
209*> \ingroup doubleGEsing
210*
211*> \par Contributors:
212*  ==================
213*>
214*>     Ming Gu and Huan Ren, Computer Science Division, University of
215*>     California at Berkeley, USA
216*>
217*  =====================================================================
218      SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
219     $                   WORK, LWORK, IWORK, INFO )
220      implicit none
221*
222*  -- LAPACK driver routine (version 3.7.0) --
223*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
224*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
225*     June 2016
226*
227*     .. Scalar Arguments ..
228      CHARACTER          JOBZ
229      INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
230*     ..
231*     .. Array Arguments ..
232      INTEGER            IWORK( * )
233      DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
234     $                   VT( LDVT, * ), WORK( * )
235*     ..
236*
237*  =====================================================================
238*
239*     .. Parameters ..
240      DOUBLE PRECISION   ZERO, ONE
241      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
242*     ..
243*     .. Local Scalars ..
244      LOGICAL            LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
245      INTEGER            BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
246     $                   IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
247     $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
248     $                   MNTHR, NWORK, WRKBL
249      INTEGER            LWORK_DGEBRD_MN, LWORK_DGEBRD_MM,
250     $                   LWORK_DGEBRD_NN, LWORK_DGELQF_MN,
251     $                   LWORK_DGEQRF_MN,
252     $                   LWORK_DORGBR_P_MM, LWORK_DORGBR_Q_NN,
253     $                   LWORK_DORGLQ_MN, LWORK_DORGLQ_NN,
254     $                   LWORK_DORGQR_MM, LWORK_DORGQR_MN,
255     $                   LWORK_DORMBR_PRT_MM, LWORK_DORMBR_QLN_MM,
256     $                   LWORK_DORMBR_PRT_MN, LWORK_DORMBR_QLN_MN,
257     $                   LWORK_DORMBR_PRT_NN, LWORK_DORMBR_QLN_NN
258      DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM
259*     ..
260*     .. Local Arrays ..
261      INTEGER            IDUM( 1 )
262      DOUBLE PRECISION   DUM( 1 )
263*     ..
264*     .. External Subroutines ..
265      EXTERNAL           DBDSDC, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY,
266     $                   DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR,
267     $                   XERBLA
268*     ..
269*     .. External Functions ..
270      LOGICAL            LSAME, DISNAN
271      DOUBLE PRECISION   DLAMCH, DLANGE
272      EXTERNAL           DLAMCH, DLANGE, LSAME, DISNAN
273*     ..
274*     .. Intrinsic Functions ..
275      INTRINSIC          INT, MAX, MIN, SQRT
276*     ..
277*     .. Executable Statements ..
278*
279*     Test the input arguments
280*
281      INFO   = 0
282      MINMN  = MIN( M, N )
283      WNTQA  = LSAME( JOBZ, 'A' )
284      WNTQS  = LSAME( JOBZ, 'S' )
285      WNTQAS = WNTQA .OR. WNTQS
286      WNTQO  = LSAME( JOBZ, 'O' )
287      WNTQN  = LSAME( JOBZ, 'N' )
288      LQUERY = ( LWORK.EQ.-1 )
289*
290      IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
291         INFO = -1
292      ELSE IF( M.LT.0 ) THEN
293         INFO = -2
294      ELSE IF( N.LT.0 ) THEN
295         INFO = -3
296      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
297         INFO = -5
298      ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
299     $         ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
300         INFO = -8
301      ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
302     $         ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
303     $         ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
304         INFO = -10
305      END IF
306*
307*     Compute workspace
308*       Note: Comments in the code beginning "Workspace:" describe the
309*       minimal amount of workspace allocated at that point in the code,
310*       as well as the preferred amount for good performance.
311*       NB refers to the optimal block size for the immediately
312*       following subroutine, as returned by ILAENV.
313*
314      IF( INFO.EQ.0 ) THEN
315         MINWRK = 1
316         MAXWRK = 1
317         BDSPAC = 0
318         MNTHR  = INT( MINMN*11.0D0 / 6.0D0 )
319         IF( M.GE.N .AND. MINMN.GT.0 ) THEN
320*
321*           Compute space needed for DBDSDC
322*
323            IF( WNTQN ) THEN
324*              dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
325*              keep 7*N for backwards compatibility.
326               BDSPAC = 7*N
327            ELSE
328               BDSPAC = 3*N*N + 4*N
329            END IF
330*
331*           Compute space preferred for each routine
332            CALL DGEBRD( M, N, DUM(1), M, DUM(1), DUM(1), DUM(1),
333     $                   DUM(1), DUM(1), -1, IERR )
334            LWORK_DGEBRD_MN = INT( DUM(1) )
335*
336            CALL DGEBRD( N, N, DUM(1), N, DUM(1), DUM(1), DUM(1),
337     $                   DUM(1), DUM(1), -1, IERR )
338            LWORK_DGEBRD_NN = INT( DUM(1) )
339*
340            CALL DGEQRF( M, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
341            LWORK_DGEQRF_MN = INT( DUM(1) )
342*
343            CALL DORGBR( 'Q', N, N, N, DUM(1), N, DUM(1), DUM(1), -1,
344     $                   IERR )
345            LWORK_DORGBR_Q_NN = INT( DUM(1) )
346*
347            CALL DORGQR( M, M, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
348            LWORK_DORGQR_MM = INT( DUM(1) )
349*
350            CALL DORGQR( M, N, N, DUM(1), M, DUM(1), DUM(1), -1, IERR )
351            LWORK_DORGQR_MN = INT( DUM(1) )
352*
353            CALL DORMBR( 'P', 'R', 'T', N, N, N, DUM(1), N,
354     $                   DUM(1), DUM(1), N, DUM(1), -1, IERR )
355            LWORK_DORMBR_PRT_NN = INT( DUM(1) )
356*
357            CALL DORMBR( 'Q', 'L', 'N', N, N, N, DUM(1), N,
358     $                   DUM(1), DUM(1), N, DUM(1), -1, IERR )
359            LWORK_DORMBR_QLN_NN = INT( DUM(1) )
360*
361            CALL DORMBR( 'Q', 'L', 'N', M, N, N, DUM(1), M,
362     $                   DUM(1), DUM(1), M, DUM(1), -1, IERR )
363            LWORK_DORMBR_QLN_MN = INT( DUM(1) )
364*
365            CALL DORMBR( 'Q', 'L', 'N', M, M, N, DUM(1), M,
366     $                   DUM(1), DUM(1), M, DUM(1), -1, IERR )
367            LWORK_DORMBR_QLN_MM = INT( DUM(1) )
368*
369            IF( M.GE.MNTHR ) THEN
370               IF( WNTQN ) THEN
371*
372*                 Path 1 (M >> N, JOBZ='N')
373*
374                  WRKBL = N + LWORK_DGEQRF_MN
375                  WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
376                  MAXWRK = MAX( WRKBL, BDSPAC + N )
377                  MINWRK = BDSPAC + N
378               ELSE IF( WNTQO ) THEN
379*
380*                 Path 2 (M >> N, JOBZ='O')
381*
382                  WRKBL = N + LWORK_DGEQRF_MN
383                  WRKBL = MAX( WRKBL,   N + LWORK_DORGQR_MN )
384                  WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
385                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
386                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
387                  WRKBL = MAX( WRKBL, 3*N + BDSPAC )
388                  MAXWRK = WRKBL + 2*N*N
389                  MINWRK = BDSPAC + 2*N*N + 3*N
390               ELSE IF( WNTQS ) THEN
391*
392*                 Path 3 (M >> N, JOBZ='S')
393*
394                  WRKBL = N + LWORK_DGEQRF_MN
395                  WRKBL = MAX( WRKBL,   N + LWORK_DORGQR_MN )
396                  WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
397                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
398                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
399                  WRKBL = MAX( WRKBL, 3*N + BDSPAC )
400                  MAXWRK = WRKBL + N*N
401                  MINWRK = BDSPAC + N*N + 3*N
402               ELSE IF( WNTQA ) THEN
403*
404*                 Path 4 (M >> N, JOBZ='A')
405*
406                  WRKBL = N + LWORK_DGEQRF_MN
407                  WRKBL = MAX( WRKBL,   N + LWORK_DORGQR_MM )
408                  WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
409                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN )
410                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
411                  WRKBL = MAX( WRKBL, 3*N + BDSPAC )
412                  MAXWRK = WRKBL + N*N
413                  MINWRK = N*N + MAX( 3*N + BDSPAC, N + M )
414               END IF
415            ELSE
416*
417*              Path 5 (M >= N, but not much larger)
418*
419               WRKBL = 3*N + LWORK_DGEBRD_MN
420               IF( WNTQN ) THEN
421*                 Path 5n (M >= N, jobz='N')
422                  MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
423                  MINWRK = 3*N + MAX( M, BDSPAC )
424               ELSE IF( WNTQO ) THEN
425*                 Path 5o (M >= N, jobz='O')
426                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
427                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MN )
428                  WRKBL = MAX( WRKBL, 3*N + BDSPAC )
429                  MAXWRK = WRKBL + M*N
430                  MINWRK = 3*N + MAX( M, N*N + BDSPAC )
431               ELSE IF( WNTQS ) THEN
432*                 Path 5s (M >= N, jobz='S')
433                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MN )
434                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
435                  MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
436                  MINWRK = 3*N + MAX( M, BDSPAC )
437               ELSE IF( WNTQA ) THEN
438*                 Path 5a (M >= N, jobz='A')
439                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MM )
440                  WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN )
441                  MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
442                  MINWRK = 3*N + MAX( M, BDSPAC )
443               END IF
444            END IF
445         ELSE IF( MINMN.GT.0 ) THEN
446*
447*           Compute space needed for DBDSDC
448*
449            IF( WNTQN ) THEN
450*              dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
451*              keep 7*N for backwards compatibility.
452               BDSPAC = 7*M
453            ELSE
454               BDSPAC = 3*M*M + 4*M
455            END IF
456*
457*           Compute space preferred for each routine
458            CALL DGEBRD( M, N, DUM(1), M, DUM(1), DUM(1), DUM(1),
459     $                   DUM(1), DUM(1), -1, IERR )
460            LWORK_DGEBRD_MN = INT( DUM(1) )
461*
462            CALL DGEBRD( M, M, A, M, S, DUM(1), DUM(1),
463     $                   DUM(1), DUM(1), -1, IERR )
464            LWORK_DGEBRD_MM = INT( DUM(1) )
465*
466            CALL DGELQF( M, N, A, M, DUM(1), DUM(1), -1, IERR )
467            LWORK_DGELQF_MN = INT( DUM(1) )
468*
469            CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR )
470            LWORK_DORGLQ_NN = INT( DUM(1) )
471*
472            CALL DORGLQ( M, N, M, A, M, DUM(1), DUM(1), -1, IERR )
473            LWORK_DORGLQ_MN = INT( DUM(1) )
474*
475            CALL DORGBR( 'P', M, M, M, A, N, DUM(1), DUM(1), -1, IERR )
476            LWORK_DORGBR_P_MM = INT( DUM(1) )
477*
478            CALL DORMBR( 'P', 'R', 'T', M, M, M, DUM(1), M,
479     $                   DUM(1), DUM(1), M, DUM(1), -1, IERR )
480            LWORK_DORMBR_PRT_MM = INT( DUM(1) )
481*
482            CALL DORMBR( 'P', 'R', 'T', M, N, M, DUM(1), M,
483     $                   DUM(1), DUM(1), M, DUM(1), -1, IERR )
484            LWORK_DORMBR_PRT_MN = INT( DUM(1) )
485*
486            CALL DORMBR( 'P', 'R', 'T', N, N, M, DUM(1), N,
487     $                   DUM(1), DUM(1), N, DUM(1), -1, IERR )
488            LWORK_DORMBR_PRT_NN = INT( DUM(1) )
489*
490            CALL DORMBR( 'Q', 'L', 'N', M, M, M, DUM(1), M,
491     $                   DUM(1), DUM(1), M, DUM(1), -1, IERR )
492            LWORK_DORMBR_QLN_MM = INT( DUM(1) )
493*
494            IF( N.GE.MNTHR ) THEN
495               IF( WNTQN ) THEN
496*
497*                 Path 1t (N >> M, JOBZ='N')
498*
499                  WRKBL = M + LWORK_DGELQF_MN
500                  WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
501                  MAXWRK = MAX( WRKBL, BDSPAC + M )
502                  MINWRK = BDSPAC + M
503               ELSE IF( WNTQO ) THEN
504*
505*                 Path 2t (N >> M, JOBZ='O')
506*
507                  WRKBL = M + LWORK_DGELQF_MN
508                  WRKBL = MAX( WRKBL,   M + LWORK_DORGLQ_MN )
509                  WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
510                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
511                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
512                  WRKBL = MAX( WRKBL, 3*M + BDSPAC )
513                  MAXWRK = WRKBL + 2*M*M
514                  MINWRK = BDSPAC + 2*M*M + 3*M
515               ELSE IF( WNTQS ) THEN
516*
517*                 Path 3t (N >> M, JOBZ='S')
518*
519                  WRKBL = M + LWORK_DGELQF_MN
520                  WRKBL = MAX( WRKBL,   M + LWORK_DORGLQ_MN )
521                  WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
522                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
523                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
524                  WRKBL = MAX( WRKBL, 3*M + BDSPAC )
525                  MAXWRK = WRKBL + M*M
526                  MINWRK = BDSPAC + M*M + 3*M
527               ELSE IF( WNTQA ) THEN
528*
529*                 Path 4t (N >> M, JOBZ='A')
530*
531                  WRKBL = M + LWORK_DGELQF_MN
532                  WRKBL = MAX( WRKBL,   M + LWORK_DORGLQ_NN )
533                  WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
534                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
535                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM )
536                  WRKBL = MAX( WRKBL, 3*M + BDSPAC )
537                  MAXWRK = WRKBL + M*M
538                  MINWRK = M*M + MAX( 3*M + BDSPAC, M + N )
539               END IF
540            ELSE
541*
542*              Path 5t (N > M, but not much larger)
543*
544               WRKBL = 3*M + LWORK_DGEBRD_MN
545               IF( WNTQN ) THEN
546*                 Path 5tn (N > M, jobz='N')
547                  MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
548                  MINWRK = 3*M + MAX( N, BDSPAC )
549               ELSE IF( WNTQO ) THEN
550*                 Path 5to (N > M, jobz='O')
551                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
552                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MN )
553                  WRKBL = MAX( WRKBL, 3*M + BDSPAC )
554                  MAXWRK = WRKBL + M*N
555                  MINWRK = 3*M + MAX( N, M*M + BDSPAC )
556               ELSE IF( WNTQS ) THEN
557*                 Path 5ts (N > M, jobz='S')
558                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
559                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MN )
560                  MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
561                  MINWRK = 3*M + MAX( N, BDSPAC )
562               ELSE IF( WNTQA ) THEN
563*                 Path 5ta (N > M, jobz='A')
564                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM )
565                  WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_NN )
566                  MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
567                  MINWRK = 3*M + MAX( N, BDSPAC )
568               END IF
569            END IF
570         END IF
571
572         MAXWRK = MAX( MAXWRK, MINWRK )
573         WORK( 1 ) = MAXWRK
574*
575         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
576            INFO = -12
577         END IF
578      END IF
579*
580      IF( INFO.NE.0 ) THEN
581         CALL XERBLA( 'DGESDD', -INFO )
582         RETURN
583      ELSE IF( LQUERY ) THEN
584         RETURN
585      END IF
586*
587*     Quick return if possible
588*
589      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
590         RETURN
591      END IF
592*
593*     Get machine constants
594*
595      EPS = DLAMCH( 'P' )
596      SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
597      BIGNUM = ONE / SMLNUM
598*
599*     Scale A if max element outside range [SMLNUM,BIGNUM]
600*
601      ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
602      IF( DISNAN( ANRM ) ) THEN
603          INFO = -4
604          RETURN
605      END IF
606      ISCL = 0
607      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
608         ISCL = 1
609         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
610      ELSE IF( ANRM.GT.BIGNUM ) THEN
611         ISCL = 1
612         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
613      END IF
614*
615      IF( M.GE.N ) THEN
616*
617*        A has at least as many rows as columns. If A has sufficiently
618*        more rows than columns, first reduce using the QR
619*        decomposition (if sufficient workspace available)
620*
621         IF( M.GE.MNTHR ) THEN
622*
623            IF( WNTQN ) THEN
624*
625*              Path 1 (M >> N, JOBZ='N')
626*              No singular vectors to be computed
627*
628               ITAU = 1
629               NWORK = ITAU + N
630*
631*              Compute A=Q*R
632*              Workspace: need   N [tau] + N    [work]
633*              Workspace: prefer N [tau] + N*NB [work]
634*
635               CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
636     $                      LWORK - NWORK + 1, IERR )
637*
638*              Zero out below R
639*
640               CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
641               IE = 1
642               ITAUQ = IE + N
643               ITAUP = ITAUQ + N
644               NWORK = ITAUP + N
645*
646*              Bidiagonalize R in A
647*              Workspace: need   3*N [e, tauq, taup] + N      [work]
648*              Workspace: prefer 3*N [e, tauq, taup] + 2*N*NB [work]
649*
650               CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
651     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
652     $                      IERR )
653               NWORK = IE + N
654*
655*              Perform bidiagonal SVD, computing singular values only
656*              Workspace: need   N [e] + BDSPAC
657*
658               CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
659     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
660*
661            ELSE IF( WNTQO ) THEN
662*
663*              Path 2 (M >> N, JOBZ = 'O')
664*              N left singular vectors to be overwritten on A and
665*              N right singular vectors to be computed in VT
666*
667               IR = 1
668*
669*              WORK(IR) is LDWRKR by N
670*
671               IF( LWORK .GE. LDA*N + N*N + 3*N + BDSPAC ) THEN
672                  LDWRKR = LDA
673               ELSE
674                  LDWRKR = ( LWORK - N*N - 3*N - BDSPAC ) / N
675               END IF
676               ITAU = IR + LDWRKR*N
677               NWORK = ITAU + N
678*
679*              Compute A=Q*R
680*              Workspace: need   N*N [R] + N [tau] + N    [work]
681*              Workspace: prefer N*N [R] + N [tau] + N*NB [work]
682*
683               CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
684     $                      LWORK - NWORK + 1, IERR )
685*
686*              Copy R to WORK(IR), zeroing out below it
687*
688               CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
689               CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
690     $                      LDWRKR )
691*
692*              Generate Q in A
693*              Workspace: need   N*N [R] + N [tau] + N    [work]
694*              Workspace: prefer N*N [R] + N [tau] + N*NB [work]
695*
696               CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
697     $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
698               IE = ITAU
699               ITAUQ = IE + N
700               ITAUP = ITAUQ + N
701               NWORK = ITAUP + N
702*
703*              Bidiagonalize R in WORK(IR)
704*              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N      [work]
705*              Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work]
706*
707               CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
708     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
709     $                      LWORK - NWORK + 1, IERR )
710*
711*              WORK(IU) is N by N
712*
713               IU = NWORK
714               NWORK = IU + N*N
715*
716*              Perform bidiagonal SVD, computing left singular vectors
717*              of bidiagonal matrix in WORK(IU) and computing right
718*              singular vectors of bidiagonal matrix in VT
719*              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N*N [U] + BDSPAC
720*
721               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
722     $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
723     $                      INFO )
724*
725*              Overwrite WORK(IU) by left singular vectors of R
726*              and VT by right singular vectors of R
727*              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N    [work]
728*              Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
729*
730               CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
731     $                      WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
732     $                      LWORK - NWORK + 1, IERR )
733               CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
734     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
735     $                      LWORK - NWORK + 1, IERR )
736*
737*              Multiply Q in A by left singular vectors of R in
738*              WORK(IU), storing result in WORK(IR) and copying to A
739*              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N*N [U]
740*              Workspace: prefer M*N [R] + 3*N [e, tauq, taup] + N*N [U]
741*
742               DO 10 I = 1, M, LDWRKR
743                  CHUNK = MIN( M - I + 1, LDWRKR )
744                  CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
745     $                        LDA, WORK( IU ), N, ZERO, WORK( IR ),
746     $                        LDWRKR )
747                  CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
748     $                         A( I, 1 ), LDA )
749   10          CONTINUE
750*
751            ELSE IF( WNTQS ) THEN
752*
753*              Path 3 (M >> N, JOBZ='S')
754*              N left singular vectors to be computed in U and
755*              N right singular vectors to be computed in VT
756*
757               IR = 1
758*
759*              WORK(IR) is N by N
760*
761               LDWRKR = N
762               ITAU = IR + LDWRKR*N
763               NWORK = ITAU + N
764*
765*              Compute A=Q*R
766*              Workspace: need   N*N [R] + N [tau] + N    [work]
767*              Workspace: prefer N*N [R] + N [tau] + N*NB [work]
768*
769               CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
770     $                      LWORK - NWORK + 1, IERR )
771*
772*              Copy R to WORK(IR), zeroing out below it
773*
774               CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
775               CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
776     $                      LDWRKR )
777*
778*              Generate Q in A
779*              Workspace: need   N*N [R] + N [tau] + N    [work]
780*              Workspace: prefer N*N [R] + N [tau] + N*NB [work]
781*
782               CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
783     $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
784               IE = ITAU
785               ITAUQ = IE + N
786               ITAUP = ITAUQ + N
787               NWORK = ITAUP + N
788*
789*              Bidiagonalize R in WORK(IR)
790*              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N      [work]
791*              Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work]
792*
793               CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
794     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
795     $                      LWORK - NWORK + 1, IERR )
796*
797*              Perform bidiagonal SVD, computing left singular vectors
798*              of bidiagoal matrix in U and computing right singular
799*              vectors of bidiagonal matrix in VT
800*              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + BDSPAC
801*
802               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
803     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
804     $                      INFO )
805*
806*              Overwrite U by left singular vectors of R and VT
807*              by right singular vectors of R
808*              Workspace: need   N*N [R] + 3*N [e, tauq, taup] + N    [work]
809*              Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*NB [work]
810*
811               CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
812     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
813     $                      LWORK - NWORK + 1, IERR )
814*
815               CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
816     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
817     $                      LWORK - NWORK + 1, IERR )
818*
819*              Multiply Q in A by left singular vectors of R in
820*              WORK(IR), storing result in U
821*              Workspace: need   N*N [R]
822*
823               CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
824               CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),
825     $                     LDWRKR, ZERO, U, LDU )
826*
827            ELSE IF( WNTQA ) THEN
828*
829*              Path 4 (M >> N, JOBZ='A')
830*              M left singular vectors to be computed in U and
831*              N right singular vectors to be computed in VT
832*
833               IU = 1
834*
835*              WORK(IU) is N by N
836*
837               LDWRKU = N
838               ITAU = IU + LDWRKU*N
839               NWORK = ITAU + N
840*
841*              Compute A=Q*R, copying result to U
842*              Workspace: need   N*N [U] + N [tau] + N    [work]
843*              Workspace: prefer N*N [U] + N [tau] + N*NB [work]
844*
845               CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
846     $                      LWORK - NWORK + 1, IERR )
847               CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
848*
849*              Generate Q in U
850*              Workspace: need   N*N [U] + N [tau] + M    [work]
851*              Workspace: prefer N*N [U] + N [tau] + M*NB [work]
852               CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
853     $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
854*
855*              Produce R in A, zeroing out other entries
856*
857               CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
858               IE = ITAU
859               ITAUQ = IE + N
860               ITAUP = ITAUQ + N
861               NWORK = ITAUP + N
862*
863*              Bidiagonalize R in A
864*              Workspace: need   N*N [U] + 3*N [e, tauq, taup] + N      [work]
865*              Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + 2*N*NB [work]
866*
867               CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
868     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
869     $                      IERR )
870*
871*              Perform bidiagonal SVD, computing left singular vectors
872*              of bidiagonal matrix in WORK(IU) and computing right
873*              singular vectors of bidiagonal matrix in VT
874*              Workspace: need   N*N [U] + 3*N [e, tauq, taup] + BDSPAC
875*
876               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
877     $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
878     $                      INFO )
879*
880*              Overwrite WORK(IU) by left singular vectors of R and VT
881*              by right singular vectors of R
882*              Workspace: need   N*N [U] + 3*N [e, tauq, taup] + N    [work]
883*              Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + N*NB [work]
884*
885               CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
886     $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
887     $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
888               CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
889     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
890     $                      LWORK - NWORK + 1, IERR )
891*
892*              Multiply Q in U by left singular vectors of R in
893*              WORK(IU), storing result in A
894*              Workspace: need   N*N [U]
895*
896               CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),
897     $                     LDWRKU, ZERO, A, LDA )
898*
899*              Copy left singular vectors of A from A to U
900*
901               CALL DLACPY( 'F', M, N, A, LDA, U, LDU )
902*
903            END IF
904*
905         ELSE
906*
907*           M .LT. MNTHR
908*
909*           Path 5 (M >= N, but not much larger)
910*           Reduce to bidiagonal form without QR decomposition
911*
912            IE = 1
913            ITAUQ = IE + N
914            ITAUP = ITAUQ + N
915            NWORK = ITAUP + N
916*
917*           Bidiagonalize A
918*           Workspace: need   3*N [e, tauq, taup] + M        [work]
919*           Workspace: prefer 3*N [e, tauq, taup] + (M+N)*NB [work]
920*
921            CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
922     $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
923     $                   IERR )
924            IF( WNTQN ) THEN
925*
926*              Path 5n (M >= N, JOBZ='N')
927*              Perform bidiagonal SVD, only computing singular values
928*              Workspace: need   3*N [e, tauq, taup] + BDSPAC
929*
930               CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
931     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
932            ELSE IF( WNTQO ) THEN
933*              Path 5o (M >= N, JOBZ='O')
934               IU = NWORK
935               IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN
936*
937*                 WORK( IU ) is M by N
938*
939                  LDWRKU = M
940                  NWORK = IU + LDWRKU*N
941                  CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),
942     $                         LDWRKU )
943*                 IR is unused; silence compile warnings
944                  IR = -1
945               ELSE
946*
947*                 WORK( IU ) is N by N
948*
949                  LDWRKU = N
950                  NWORK = IU + LDWRKU*N
951*
952*                 WORK(IR) is LDWRKR by N
953*
954                  IR = NWORK
955                  LDWRKR = ( LWORK - N*N - 3*N ) / N
956               END IF
957               NWORK = IU + LDWRKU*N
958*
959*              Perform bidiagonal SVD, computing left singular vectors
960*              of bidiagonal matrix in WORK(IU) and computing right
961*              singular vectors of bidiagonal matrix in VT
962*              Workspace: need   3*N [e, tauq, taup] + N*N [U] + BDSPAC
963*
964               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
965     $                      LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
966     $                      IWORK, INFO )
967*
968*              Overwrite VT by right singular vectors of A
969*              Workspace: need   3*N [e, tauq, taup] + N*N [U] + N    [work]
970*              Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
971*
972               CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
973     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
974     $                      LWORK - NWORK + 1, IERR )
975*
976               IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN
977*
978*                 Path 5o-fast
979*                 Overwrite WORK(IU) by left singular vectors of A
980*                 Workspace: need   3*N [e, tauq, taup] + M*N [U] + N    [work]
981*                 Workspace: prefer 3*N [e, tauq, taup] + M*N [U] + N*NB [work]
982*
983                  CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
984     $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,
985     $                         WORK( NWORK ), LWORK - NWORK + 1, IERR )
986*
987*                 Copy left singular vectors of A from WORK(IU) to A
988*
989                  CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
990               ELSE
991*
992*                 Path 5o-slow
993*                 Generate Q in A
994*                 Workspace: need   3*N [e, tauq, taup] + N*N [U] + N    [work]
995*                 Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work]
996*
997                  CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
998     $                         WORK( NWORK ), LWORK - NWORK + 1, IERR )
999*
1000*                 Multiply Q in A by left singular vectors of
1001*                 bidiagonal matrix in WORK(IU), storing result in
1002*                 WORK(IR) and copying to A
1003*                 Workspace: need   3*N [e, tauq, taup] + N*N [U] + NB*N [R]
1004*                 Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + M*N  [R]
1005*
1006                  DO 20 I = 1, M, LDWRKR
1007                     CHUNK = MIN( M - I + 1, LDWRKR )
1008                     CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
1009     $                           LDA, WORK( IU ), LDWRKU, ZERO,
1010     $                           WORK( IR ), LDWRKR )
1011                     CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
1012     $                            A( I, 1 ), LDA )
1013   20             CONTINUE
1014               END IF
1015*
1016            ELSE IF( WNTQS ) THEN
1017*
1018*              Path 5s (M >= N, JOBZ='S')
1019*              Perform bidiagonal SVD, computing left singular vectors
1020*              of bidiagonal matrix in U and computing right singular
1021*              vectors of bidiagonal matrix in VT
1022*              Workspace: need   3*N [e, tauq, taup] + BDSPAC
1023*
1024               CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU )
1025               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
1026     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1027     $                      INFO )
1028*
1029*              Overwrite U by left singular vectors of A and VT
1030*              by right singular vectors of A
1031*              Workspace: need   3*N [e, tauq, taup] + N    [work]
1032*              Workspace: prefer 3*N [e, tauq, taup] + N*NB [work]
1033*
1034               CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1035     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1036     $                      LWORK - NWORK + 1, IERR )
1037               CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
1038     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1039     $                      LWORK - NWORK + 1, IERR )
1040            ELSE IF( WNTQA ) THEN
1041*
1042*              Path 5a (M >= N, JOBZ='A')
1043*              Perform bidiagonal SVD, computing left singular vectors
1044*              of bidiagonal matrix in U and computing right singular
1045*              vectors of bidiagonal matrix in VT
1046*              Workspace: need   3*N [e, tauq, taup] + BDSPAC
1047*
1048               CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU )
1049               CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
1050     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1051     $                      INFO )
1052*
1053*              Set the right corner of U to identity matrix
1054*
1055               IF( M.GT.N ) THEN
1056                  CALL DLASET( 'F', M - N, M - N, ZERO, ONE, U(N+1,N+1),
1057     $                         LDU )
1058               END IF
1059*
1060*              Overwrite U by left singular vectors of A and VT
1061*              by right singular vectors of A
1062*              Workspace: need   3*N [e, tauq, taup] + M    [work]
1063*              Workspace: prefer 3*N [e, tauq, taup] + M*NB [work]
1064*
1065               CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1066     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1067     $                      LWORK - NWORK + 1, IERR )
1068               CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1069     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1070     $                      LWORK - NWORK + 1, IERR )
1071            END IF
1072*
1073         END IF
1074*
1075      ELSE
1076*
1077*        A has more columns than rows. If A has sufficiently more
1078*        columns than rows, first reduce using the LQ decomposition (if
1079*        sufficient workspace available)
1080*
1081         IF( N.GE.MNTHR ) THEN
1082*
1083            IF( WNTQN ) THEN
1084*
1085*              Path 1t (N >> M, JOBZ='N')
1086*              No singular vectors to be computed
1087*
1088               ITAU = 1
1089               NWORK = ITAU + M
1090*
1091*              Compute A=L*Q
1092*              Workspace: need   M [tau] + M [work]
1093*              Workspace: prefer M [tau] + M*NB [work]
1094*
1095               CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1096     $                      LWORK - NWORK + 1, IERR )
1097*
1098*              Zero out above L
1099*
1100               CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1101               IE = 1
1102               ITAUQ = IE + M
1103               ITAUP = ITAUQ + M
1104               NWORK = ITAUP + M
1105*
1106*              Bidiagonalize L in A
1107*              Workspace: need   3*M [e, tauq, taup] + M      [work]
1108*              Workspace: prefer 3*M [e, tauq, taup] + 2*M*NB [work]
1109*
1110               CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1111     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1112     $                      IERR )
1113               NWORK = IE + M
1114*
1115*              Perform bidiagonal SVD, computing singular values only
1116*              Workspace: need   M [e] + BDSPAC
1117*
1118               CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1119     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1120*
1121            ELSE IF( WNTQO ) THEN
1122*
1123*              Path 2t (N >> M, JOBZ='O')
1124*              M right singular vectors to be overwritten on A and
1125*              M left singular vectors to be computed in U
1126*
1127               IVT = 1
1128*
1129*              WORK(IVT) is M by M
1130*              WORK(IL)  is M by M; it is later resized to M by chunk for gemm
1131*
1132               IL = IVT + M*M
1133               IF( LWORK .GE. M*N + M*M + 3*M + BDSPAC ) THEN
1134                  LDWRKL = M
1135                  CHUNK = N
1136               ELSE
1137                  LDWRKL = M
1138                  CHUNK = ( LWORK - M*M ) / M
1139               END IF
1140               ITAU = IL + LDWRKL*M
1141               NWORK = ITAU + M
1142*
1143*              Compute A=L*Q
1144*              Workspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work]
1145*              Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1146*
1147               CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1148     $                      LWORK - NWORK + 1, IERR )
1149*
1150*              Copy L to WORK(IL), zeroing about above it
1151*
1152               CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1153               CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
1154     $                      WORK( IL + LDWRKL ), LDWRKL )
1155*
1156*              Generate Q in A
1157*              Workspace: need   M*M [VT] + M*M [L] + M [tau] + M    [work]
1158*              Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1159*
1160               CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1161     $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
1162               IE = ITAU
1163               ITAUQ = IE + M
1164               ITAUP = ITAUQ + M
1165               NWORK = ITAUP + M
1166*
1167*              Bidiagonalize L in WORK(IL)
1168*              Workspace: need   M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M      [work]
1169*              Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work]
1170*
1171               CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1172     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1173     $                      LWORK - NWORK + 1, IERR )
1174*
1175*              Perform bidiagonal SVD, computing left singular vectors
1176*              of bidiagonal matrix in U, and computing right singular
1177*              vectors of bidiagonal matrix in WORK(IVT)
1178*              Workspace: need   M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1179*
1180               CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1181     $                      WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
1182     $                      IWORK, INFO )
1183*
1184*              Overwrite U by left singular vectors of L and WORK(IVT)
1185*              by right singular vectors of L
1186*              Workspace: need   M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M    [work]
1187*              Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M*NB [work]
1188*
1189               CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1190     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1191     $                      LWORK - NWORK + 1, IERR )
1192               CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1193     $                      WORK( ITAUP ), WORK( IVT ), M,
1194     $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
1195*
1196*              Multiply right singular vectors of L in WORK(IVT) by Q
1197*              in A, storing result in WORK(IL) and copying to A
1198*              Workspace: need   M*M [VT] + M*M [L]
1199*              Workspace: prefer M*M [VT] + M*N [L]
1200*              At this point, L is resized as M by chunk.
1201*
1202               DO 30 I = 1, N, CHUNK
1203                  BLK = MIN( N - I + 1, CHUNK )
1204                  CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,
1205     $                        A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
1206                  CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
1207     $                         A( 1, I ), LDA )
1208   30          CONTINUE
1209*
1210            ELSE IF( WNTQS ) THEN
1211*
1212*              Path 3t (N >> M, JOBZ='S')
1213*              M right singular vectors to be computed in VT and
1214*              M left singular vectors to be computed in U
1215*
1216               IL = 1
1217*
1218*              WORK(IL) is M by M
1219*
1220               LDWRKL = M
1221               ITAU = IL + LDWRKL*M
1222               NWORK = ITAU + M
1223*
1224*              Compute A=L*Q
1225*              Workspace: need   M*M [L] + M [tau] + M    [work]
1226*              Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1227*
1228               CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1229     $                      LWORK - NWORK + 1, IERR )
1230*
1231*              Copy L to WORK(IL), zeroing out above it
1232*
1233               CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1234               CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
1235     $                      WORK( IL + LDWRKL ), LDWRKL )
1236*
1237*              Generate Q in A
1238*              Workspace: need   M*M [L] + M [tau] + M    [work]
1239*              Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1240*
1241               CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1242     $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
1243               IE = ITAU
1244               ITAUQ = IE + M
1245               ITAUP = ITAUQ + M
1246               NWORK = ITAUP + M
1247*
1248*              Bidiagonalize L in WORK(IU).
1249*              Workspace: need   M*M [L] + 3*M [e, tauq, taup] + M      [work]
1250*              Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work]
1251*
1252               CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1253     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1254     $                      LWORK - NWORK + 1, IERR )
1255*
1256*              Perform bidiagonal SVD, computing left singular vectors
1257*              of bidiagonal matrix in U and computing right singular
1258*              vectors of bidiagonal matrix in VT
1259*              Workspace: need   M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1260*
1261               CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
1262     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1263     $                      INFO )
1264*
1265*              Overwrite U by left singular vectors of L and VT
1266*              by right singular vectors of L
1267*              Workspace: need   M*M [L] + 3*M [e, tauq, taup] + M    [work]
1268*              Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + M*NB [work]
1269*
1270               CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1271     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1272     $                      LWORK - NWORK + 1, IERR )
1273               CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1274     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1275     $                      LWORK - NWORK + 1, IERR )
1276*
1277*              Multiply right singular vectors of L in WORK(IL) by
1278*              Q in A, storing result in VT
1279*              Workspace: need   M*M [L]
1280*
1281               CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1282               CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,
1283     $                     A, LDA, ZERO, VT, LDVT )
1284*
1285            ELSE IF( WNTQA ) THEN
1286*
1287*              Path 4t (N >> M, JOBZ='A')
1288*              N right singular vectors to be computed in VT and
1289*              M left singular vectors to be computed in U
1290*
1291               IVT = 1
1292*
1293*              WORK(IVT) is M by M
1294*
1295               LDWKVT = M
1296               ITAU = IVT + LDWKVT*M
1297               NWORK = ITAU + M
1298*
1299*              Compute A=L*Q, copying result to VT
1300*              Workspace: need   M*M [VT] + M [tau] + M    [work]
1301*              Workspace: prefer M*M [VT] + M [tau] + M*NB [work]
1302*
1303               CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1304     $                      LWORK - NWORK + 1, IERR )
1305               CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
1306*
1307*              Generate Q in VT
1308*              Workspace: need   M*M [VT] + M [tau] + N    [work]
1309*              Workspace: prefer M*M [VT] + M [tau] + N*NB [work]
1310*
1311               CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1312     $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
1313*
1314*              Produce L in A, zeroing out other entries
1315*
1316               CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1317               IE = ITAU
1318               ITAUQ = IE + M
1319               ITAUP = ITAUQ + M
1320               NWORK = ITAUP + M
1321*
1322*              Bidiagonalize L in A
1323*              Workspace: need   M*M [VT] + 3*M [e, tauq, taup] + M      [work]
1324*              Workspace: prefer M*M [VT] + 3*M [e, tauq, taup] + 2*M*NB [work]
1325*
1326               CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1327     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1328     $                      IERR )
1329*
1330*              Perform bidiagonal SVD, computing left singular vectors
1331*              of bidiagonal matrix in U and computing right singular
1332*              vectors of bidiagonal matrix in WORK(IVT)
1333*              Workspace: need   M*M [VT] + 3*M [e, tauq, taup] + BDSPAC
1334*
1335               CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1336     $                      WORK( IVT ), LDWKVT, DUM, IDUM,
1337     $                      WORK( NWORK ), IWORK, INFO )
1338*
1339*              Overwrite U by left singular vectors of L and WORK(IVT)
1340*              by right singular vectors of L
1341*              Workspace: need   M*M [VT] + 3*M [e, tauq, taup]+ M    [work]
1342*              Workspace: prefer M*M [VT] + 3*M [e, tauq, taup]+ M*NB [work]
1343*
1344               CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1345     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1346     $                      LWORK - NWORK + 1, IERR )
1347               CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
1348     $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
1349     $                      WORK( NWORK ), LWORK - NWORK + 1, IERR )
1350*
1351*              Multiply right singular vectors of L in WORK(IVT) by
1352*              Q in VT, storing result in A
1353*              Workspace: need   M*M [VT]
1354*
1355               CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,
1356     $                     VT, LDVT, ZERO, A, LDA )
1357*
1358*              Copy right singular vectors of A from A to VT
1359*
1360               CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT )
1361*
1362            END IF
1363*
1364         ELSE
1365*
1366*           N .LT. MNTHR
1367*
1368*           Path 5t (N > M, but not much larger)
1369*           Reduce to bidiagonal form without LQ decomposition
1370*
1371            IE = 1
1372            ITAUQ = IE + M
1373            ITAUP = ITAUQ + M
1374            NWORK = ITAUP + M
1375*
1376*           Bidiagonalize A
1377*           Workspace: need   3*M [e, tauq, taup] + N        [work]
1378*           Workspace: prefer 3*M [e, tauq, taup] + (M+N)*NB [work]
1379*
1380            CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1381     $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1382     $                   IERR )
1383            IF( WNTQN ) THEN
1384*
1385*              Path 5tn (N > M, JOBZ='N')
1386*              Perform bidiagonal SVD, only computing singular values
1387*              Workspace: need   3*M [e, tauq, taup] + BDSPAC
1388*
1389               CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1390     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1391            ELSE IF( WNTQO ) THEN
1392*              Path 5to (N > M, JOBZ='O')
1393               LDWKVT = M
1394               IVT = NWORK
1395               IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN
1396*
1397*                 WORK( IVT ) is M by N
1398*
1399                  CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),
1400     $                         LDWKVT )
1401                  NWORK = IVT + LDWKVT*N
1402*                 IL is unused; silence compile warnings
1403                  IL = -1
1404               ELSE
1405*
1406*                 WORK( IVT ) is M by M
1407*
1408                  NWORK = IVT + LDWKVT*M
1409                  IL = NWORK
1410*
1411*                 WORK(IL) is M by CHUNK
1412*
1413                  CHUNK = ( LWORK - M*M - 3*M ) / M
1414               END IF
1415*
1416*              Perform bidiagonal SVD, computing left singular vectors
1417*              of bidiagonal matrix in U and computing right singular
1418*              vectors of bidiagonal matrix in WORK(IVT)
1419*              Workspace: need   3*M [e, tauq, taup] + M*M [VT] + BDSPAC
1420*
1421               CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,
1422     $                      WORK( IVT ), LDWKVT, DUM, IDUM,
1423     $                      WORK( NWORK ), IWORK, INFO )
1424*
1425*              Overwrite U by left singular vectors of A
1426*              Workspace: need   3*M [e, tauq, taup] + M*M [VT] + M    [work]
1427*              Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work]
1428*
1429               CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1430     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1431     $                      LWORK - NWORK + 1, IERR )
1432*
1433               IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN
1434*
1435*                 Path 5to-fast
1436*                 Overwrite WORK(IVT) by left singular vectors of A
1437*                 Workspace: need   3*M [e, tauq, taup] + M*N [VT] + M    [work]
1438*                 Workspace: prefer 3*M [e, tauq, taup] + M*N [VT] + M*NB [work]
1439*
1440                  CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1441     $                         WORK( ITAUP ), WORK( IVT ), LDWKVT,
1442     $                         WORK( NWORK ), LWORK - NWORK + 1, IERR )
1443*
1444*                 Copy right singular vectors of A from WORK(IVT) to A
1445*
1446                  CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
1447               ELSE
1448*
1449*                 Path 5to-slow
1450*                 Generate P**T in A
1451*                 Workspace: need   3*M [e, tauq, taup] + M*M [VT] + M    [work]
1452*                 Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work]
1453*
1454                  CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1455     $                         WORK( NWORK ), LWORK - NWORK + 1, IERR )
1456*
1457*                 Multiply Q in A by right singular vectors of
1458*                 bidiagonal matrix in WORK(IVT), storing result in
1459*                 WORK(IL) and copying to A
1460*                 Workspace: need   3*M [e, tauq, taup] + M*M [VT] + M*NB [L]
1461*                 Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*N  [L]
1462*
1463                  DO 40 I = 1, N, CHUNK
1464                     BLK = MIN( N - I + 1, CHUNK )
1465                     CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),
1466     $                           LDWKVT, A( 1, I ), LDA, ZERO,
1467     $                           WORK( IL ), M )
1468                     CALL DLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ),
1469     $                            LDA )
1470   40             CONTINUE
1471               END IF
1472            ELSE IF( WNTQS ) THEN
1473*
1474*              Path 5ts (N > M, JOBZ='S')
1475*              Perform bidiagonal SVD, computing left singular vectors
1476*              of bidiagonal matrix in U and computing right singular
1477*              vectors of bidiagonal matrix in VT
1478*              Workspace: need   3*M [e, tauq, taup] + BDSPAC
1479*
1480               CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
1481               CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1482     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1483     $                      INFO )
1484*
1485*              Overwrite U by left singular vectors of A and VT
1486*              by right singular vectors of A
1487*              Workspace: need   3*M [e, tauq, taup] + M    [work]
1488*              Workspace: prefer 3*M [e, tauq, taup] + M*NB [work]
1489*
1490               CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1491     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1492     $                      LWORK - NWORK + 1, IERR )
1493               CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1494     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1495     $                      LWORK - NWORK + 1, IERR )
1496            ELSE IF( WNTQA ) THEN
1497*
1498*              Path 5ta (N > M, JOBZ='A')
1499*              Perform bidiagonal SVD, computing left singular vectors
1500*              of bidiagonal matrix in U and computing right singular
1501*              vectors of bidiagonal matrix in VT
1502*              Workspace: need   3*M [e, tauq, taup] + BDSPAC
1503*
1504               CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
1505               CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1506     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1507     $                      INFO )
1508*
1509*              Set the right corner of VT to identity matrix
1510*
1511               IF( N.GT.M ) THEN
1512                  CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT(M+1,M+1),
1513     $                         LDVT )
1514               END IF
1515*
1516*              Overwrite U by left singular vectors of A and VT
1517*              by right singular vectors of A
1518*              Workspace: need   3*M [e, tauq, taup] + N    [work]
1519*              Workspace: prefer 3*M [e, tauq, taup] + N*NB [work]
1520*
1521               CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1522     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1523     $                      LWORK - NWORK + 1, IERR )
1524               CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1525     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1526     $                      LWORK - NWORK + 1, IERR )
1527            END IF
1528*
1529         END IF
1530*
1531      END IF
1532*
1533*     Undo scaling if necessary
1534*
1535      IF( ISCL.EQ.1 ) THEN
1536         IF( ANRM.GT.BIGNUM )
1537     $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
1538     $                   IERR )
1539         IF( ANRM.LT.SMLNUM )
1540     $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
1541     $                   IERR )
1542      END IF
1543*
1544*     Return optimal workspace in WORK(1)
1545*
1546      WORK( 1 ) = MAXWRK
1547*
1548      RETURN
1549*
1550*     End of DGESDD
1551*
1552      END
1553