1*> \brief \b SGESDD
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SGESDD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgesdd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgesdd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgesdd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
22*                          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*       REAL               A( LDA, * ), S( * ), U( LDU, * ),
31*      $                   VT( LDVT, * ), WORK( * )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> SGESDD 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 REAL 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 REAL 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 REAL 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 REAL 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; if
158*>          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 REAL 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 JOBZ = 'N',
173*>            LWORK >= 3*min(M,N) + max(max(M,N),6*min(M,N)).
174*>          If JOBZ = 'O',
175*>            LWORK >= 3*min(M,N) +
176*>                     max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
177*>          If JOBZ = 'S' or 'A'
178*>            LWORK >= min(M,N)*(7+4*min(M,N))
179*>          For good performance, LWORK should generally be larger.
180*>          If LWORK = -1 but other input arguments are legal, WORK(1)
181*>          returns the optimal LWORK.
182*> \endverbatim
183*>
184*> \param[out] IWORK
185*> \verbatim
186*>          IWORK is INTEGER array, dimension (8*min(M,N))
187*> \endverbatim
188*>
189*> \param[out] INFO
190*> \verbatim
191*>          INFO is INTEGER
192*>          = 0:  successful exit.
193*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
194*>          > 0:  SBDSDC did not converge, updating process failed.
195*> \endverbatim
196*
197*  Authors:
198*  ========
199*
200*> \author Univ. of Tennessee
201*> \author Univ. of California Berkeley
202*> \author Univ. of Colorado Denver
203*> \author NAG Ltd.
204*
205*> \date November 2015
206*
207*> \ingroup realGEsing
208*
209*> \par Contributors:
210*  ==================
211*>
212*>     Ming Gu and Huan Ren, Computer Science Division, University of
213*>     California at Berkeley, USA
214*>
215*  =====================================================================
216      SUBROUTINE SGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
217     $                   LWORK, IWORK, INFO )
218*
219*  -- LAPACK driver routine (version 3.6.0) --
220*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
221*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
222*     November 2015
223*
224*     .. Scalar Arguments ..
225      CHARACTER          JOBZ
226      INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
227*     ..
228*     .. Array Arguments ..
229      INTEGER            IWORK( * )
230      REAL               A( LDA, * ), S( * ), U( LDU, * ),
231     $                   VT( LDVT, * ), WORK( * )
232*     ..
233*
234*  =====================================================================
235*
236*     .. Parameters ..
237      REAL               ZERO, ONE
238      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
239*     ..
240*     .. Local Scalars ..
241      LOGICAL            LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
242      INTEGER            BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
243     $                   IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
244     $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
245     $                   MNTHR, NWORK, WRKBL
246      REAL               ANRM, BIGNUM, EPS, SMLNUM
247*     ..
248*     .. Local Arrays ..
249      INTEGER            IDUM( 1 )
250      REAL               DUM( 1 )
251*     ..
252*     .. External Subroutines ..
253      EXTERNAL           SBDSDC, SGEBRD, SGELQF, SGEMM, SGEQRF, SLACPY,
254     $                   SLASCL, SLASET, SORGBR, SORGLQ, SORGQR, SORMBR,
255     $                   XERBLA
256*     ..
257*     .. External Functions ..
258      LOGICAL            LSAME
259      INTEGER            ILAENV
260      REAL               SLAMCH, SLANGE
261      EXTERNAL           ILAENV, LSAME, SLAMCH, SLANGE
262*     ..
263*     .. Intrinsic Functions ..
264      INTRINSIC          INT, MAX, MIN, SQRT
265*     ..
266*     .. Executable Statements ..
267*
268*     Test the input arguments
269*
270      INFO = 0
271      MINMN = MIN( M, N )
272      WNTQA = LSAME( JOBZ, 'A' )
273      WNTQS = LSAME( JOBZ, 'S' )
274      WNTQAS = WNTQA .OR. WNTQS
275      WNTQO = LSAME( JOBZ, 'O' )
276      WNTQN = LSAME( JOBZ, 'N' )
277      LQUERY = ( LWORK.EQ.-1 )
278*
279      IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
280         INFO = -1
281      ELSE IF( M.LT.0 ) THEN
282         INFO = -2
283      ELSE IF( N.LT.0 ) THEN
284         INFO = -3
285      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
286         INFO = -5
287      ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
288     $         ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
289         INFO = -8
290      ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
291     $         ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
292     $         ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
293         INFO = -10
294      END IF
295*
296*     Compute workspace
297*      (Note: Comments in the code beginning "Workspace:" describe the
298*       minimal amount of workspace needed at that point in the code,
299*       as well as the preferred amount for good performance.
300*       NB refers to the optimal block size for the immediately
301*       following subroutine, as returned by ILAENV.)
302*
303      IF( INFO.EQ.0 ) THEN
304         MINWRK = 1
305         MAXWRK = 1
306         IF( M.GE.N .AND. MINMN.GT.0 ) THEN
307*
308*           Compute space needed for SBDSDC
309*
310            MNTHR = INT( MINMN*11.0E0 / 6.0E0 )
311            IF( WNTQN ) THEN
312               BDSPAC = 7*N
313            ELSE
314               BDSPAC = 3*N*N + 4*N
315            END IF
316            IF( M.GE.MNTHR ) THEN
317               IF( WNTQN ) THEN
318*
319*                 Path 1 (M much larger than N, JOBZ='N')
320*
321                  WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1,
322     $                    -1 )
323                  WRKBL = MAX( WRKBL, 3*N+2*N*
324     $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
325                  MAXWRK = MAX( WRKBL, BDSPAC+N )
326                  MINWRK = BDSPAC + N
327               ELSE IF( WNTQO ) THEN
328*
329*                 Path 2 (M much larger than N, JOBZ='O')
330*
331                  WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
332                  WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M,
333     $                    N, N, -1 ) )
334                  WRKBL = MAX( WRKBL, 3*N+2*N*
335     $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
336                  WRKBL = MAX( WRKBL, 3*N+N*
337     $                    ILAENV( 1, 'SORMBR', 'QLN', N, N, N, -1 ) )
338                  WRKBL = MAX( WRKBL, 3*N+N*
339     $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
340                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
341                  MAXWRK = WRKBL + 2*N*N
342                  MINWRK = BDSPAC + 2*N*N + 3*N
343               ELSE IF( WNTQS ) THEN
344*
345*                 Path 3 (M much larger than N, JOBZ='S')
346*
347                  WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
348                  WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'SORGQR', ' ', M,
349     $                    N, N, -1 ) )
350                  WRKBL = MAX( WRKBL, 3*N+2*N*
351     $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
352                  WRKBL = MAX( WRKBL, 3*N+N*
353     $                    ILAENV( 1, 'SORMBR', 'QLN', N, N, N, -1 ) )
354                  WRKBL = MAX( WRKBL, 3*N+N*
355     $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
356                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
357                  MAXWRK = WRKBL + N*N
358                  MINWRK = BDSPAC + N*N + 3*N
359               ELSE IF( WNTQA ) THEN
360*
361*                 Path 4 (M much larger than N, JOBZ='A')
362*
363                  WRKBL = N + N*ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
364                  WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'SORGQR', ' ', M,
365     $                    M, N, -1 ) )
366                  WRKBL = MAX( WRKBL, 3*N+2*N*
367     $                    ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
368                  WRKBL = MAX( WRKBL, 3*N+N*
369     $                    ILAENV( 1, 'SORMBR', 'QLN', N, N, N, -1 ) )
370                  WRKBL = MAX( WRKBL, 3*N+N*
371     $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
372                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
373                  MAXWRK = WRKBL + N*N
374                  MINWRK = BDSPAC + N*N + 2*N + M
375               END IF
376            ELSE
377*
378*              Path 5 (M at least N, but not much larger)
379*
380               WRKBL = 3*N + ( M+N )*ILAENV( 1, 'SGEBRD', ' ', M, N, -1,
381     $                 -1 )
382               IF( WNTQN ) THEN
383                  MAXWRK = MAX( WRKBL, BDSPAC+3*N )
384                  MINWRK = 3*N + MAX( M, BDSPAC )
385               ELSE IF( WNTQO ) THEN
386                  WRKBL = MAX( WRKBL, 3*N+N*
387     $                    ILAENV( 1, 'SORMBR', 'QLN', M, N, N, -1 ) )
388                  WRKBL = MAX( WRKBL, 3*N+N*
389     $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
390                  WRKBL = MAX( WRKBL, BDSPAC+3*N )
391                  MAXWRK = WRKBL + M*N
392                  MINWRK = 3*N + MAX( M, N*N+BDSPAC )
393               ELSE IF( WNTQS ) THEN
394                  WRKBL = MAX( WRKBL, 3*N+N*
395     $                    ILAENV( 1, 'SORMBR', 'QLN', M, N, N, -1 ) )
396                  WRKBL = MAX( WRKBL, 3*N+N*
397     $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
398                  MAXWRK = MAX( WRKBL, BDSPAC+3*N )
399                  MINWRK = 3*N + MAX( M, BDSPAC )
400               ELSE IF( WNTQA ) THEN
401                  WRKBL = MAX( WRKBL, 3*N+M*
402     $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, N, -1 ) )
403                  WRKBL = MAX( WRKBL, 3*N+N*
404     $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, N, -1 ) )
405                  MAXWRK = MAX( MAXWRK, BDSPAC+3*N )
406                  MINWRK = 3*N + MAX( M, BDSPAC )
407               END IF
408            END IF
409         ELSE IF ( MINMN.GT.0 ) THEN
410*
411*           Compute space needed for SBDSDC
412*
413            MNTHR = INT( MINMN*11.0E0 / 6.0E0 )
414            IF( WNTQN ) THEN
415               BDSPAC = 7*M
416            ELSE
417               BDSPAC = 3*M*M + 4*M
418            END IF
419            IF( N.GE.MNTHR ) THEN
420               IF( WNTQN ) THEN
421*
422*                 Path 1t (N much larger than M, JOBZ='N')
423*
424                  WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1,
425     $                    -1 )
426                  WRKBL = MAX( WRKBL, 3*M+2*M*
427     $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
428                  MAXWRK = MAX( WRKBL, BDSPAC+M )
429                  MINWRK = BDSPAC + M
430               ELSE IF( WNTQO ) THEN
431*
432*                 Path 2t (N much larger than M, JOBZ='O')
433*
434                  WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
435                  WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M,
436     $                    N, M, -1 ) )
437                  WRKBL = MAX( WRKBL, 3*M+2*M*
438     $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
439                  WRKBL = MAX( WRKBL, 3*M+M*
440     $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, M, -1 ) )
441                  WRKBL = MAX( WRKBL, 3*M+M*
442     $                    ILAENV( 1, 'SORMBR', 'PRT', M, M, M, -1 ) )
443                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
444                  MAXWRK = WRKBL + 2*M*M
445                  MINWRK = BDSPAC + 2*M*M + 3*M
446               ELSE IF( WNTQS ) THEN
447*
448*                 Path 3t (N much larger than M, JOBZ='S')
449*
450                  WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
451                  WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'SORGLQ', ' ', M,
452     $                    N, M, -1 ) )
453                  WRKBL = MAX( WRKBL, 3*M+2*M*
454     $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
455                  WRKBL = MAX( WRKBL, 3*M+M*
456     $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, M, -1 ) )
457                  WRKBL = MAX( WRKBL, 3*M+M*
458     $                    ILAENV( 1, 'SORMBR', 'PRT', M, M, M, -1 ) )
459                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
460                  MAXWRK = WRKBL + M*M
461                  MINWRK = BDSPAC + M*M + 3*M
462               ELSE IF( WNTQA ) THEN
463*
464*                 Path 4t (N much larger than M, JOBZ='A')
465*
466                  WRKBL = M + M*ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
467                  WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'SORGLQ', ' ', N,
468     $                    N, M, -1 ) )
469                  WRKBL = MAX( WRKBL, 3*M+2*M*
470     $                    ILAENV( 1, 'SGEBRD', ' ', M, M, -1, -1 ) )
471                  WRKBL = MAX( WRKBL, 3*M+M*
472     $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, M, -1 ) )
473                  WRKBL = MAX( WRKBL, 3*M+M*
474     $                    ILAENV( 1, 'SORMBR', 'PRT', M, M, M, -1 ) )
475                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
476                  MAXWRK = WRKBL + M*M
477                  MINWRK = BDSPAC + M*M + 3*M
478               END IF
479            ELSE
480*
481*              Path 5t (N greater than M, but not much larger)
482*
483               WRKBL = 3*M + ( M+N )*ILAENV( 1, 'SGEBRD', ' ', M, N, -1,
484     $                 -1 )
485               IF( WNTQN ) THEN
486                  MAXWRK = MAX( WRKBL, BDSPAC+3*M )
487                  MINWRK = 3*M + MAX( N, BDSPAC )
488               ELSE IF( WNTQO ) THEN
489                  WRKBL = MAX( WRKBL, 3*M+M*
490     $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, N, -1 ) )
491                  WRKBL = MAX( WRKBL, 3*M+M*
492     $                    ILAENV( 1, 'SORMBR', 'PRT', M, N, M, -1 ) )
493                  WRKBL = MAX( WRKBL, BDSPAC+3*M )
494                  MAXWRK = WRKBL + M*N
495                  MINWRK = 3*M + MAX( N, M*M+BDSPAC )
496               ELSE IF( WNTQS ) THEN
497                  WRKBL = MAX( WRKBL, 3*M+M*
498     $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, N, -1 ) )
499                  WRKBL = MAX( WRKBL, 3*M+M*
500     $                    ILAENV( 1, 'SORMBR', 'PRT', M, N, M, -1 ) )
501                  MAXWRK = MAX( WRKBL, BDSPAC+3*M )
502                  MINWRK = 3*M + MAX( N, BDSPAC )
503               ELSE IF( WNTQA ) THEN
504                  WRKBL = MAX( WRKBL, 3*M+M*
505     $                    ILAENV( 1, 'SORMBR', 'QLN', M, M, N, -1 ) )
506                  WRKBL = MAX( WRKBL, 3*M+M*
507     $                    ILAENV( 1, 'SORMBR', 'PRT', N, N, M, -1 ) )
508                  MAXWRK = MAX( WRKBL, BDSPAC+3*M )
509                  MINWRK = 3*M + MAX( N, BDSPAC )
510               END IF
511            END IF
512         END IF
513         MAXWRK = MAX( MAXWRK, MINWRK )
514         WORK( 1 ) = MAXWRK
515*
516         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
517            INFO = -12
518         END IF
519      END IF
520*
521      IF( INFO.NE.0 ) THEN
522         CALL XERBLA( 'SGESDD', -INFO )
523         RETURN
524      ELSE IF( LQUERY ) THEN
525         RETURN
526      END IF
527*
528*     Quick return if possible
529*
530      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
531         RETURN
532      END IF
533*
534*     Get machine constants
535*
536      EPS = SLAMCH( 'P' )
537      SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS
538      BIGNUM = ONE / SMLNUM
539*
540*     Scale A if max element outside range [SMLNUM,BIGNUM]
541*
542      ANRM = SLANGE( 'M', M, N, A, LDA, DUM )
543      ISCL = 0
544      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
545         ISCL = 1
546         CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
547      ELSE IF( ANRM.GT.BIGNUM ) THEN
548         ISCL = 1
549         CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
550      END IF
551*
552      IF( M.GE.N ) THEN
553*
554*        A has at least as many rows as columns. If A has sufficiently
555*        more rows than columns, first reduce using the QR
556*        decomposition (if sufficient workspace available)
557*
558         IF( M.GE.MNTHR ) THEN
559*
560            IF( WNTQN ) THEN
561*
562*              Path 1 (M much larger than N, JOBZ='N')
563*              No singular vectors to be computed
564*
565               ITAU = 1
566               NWORK = ITAU + N
567*
568*              Compute A=Q*R
569*              (Workspace: need 2*N, prefer N+N*NB)
570*
571               CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
572     $                      LWORK-NWORK+1, IERR )
573*
574*              Zero out below R
575*
576               CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
577               IE = 1
578               ITAUQ = IE + N
579               ITAUP = ITAUQ + N
580               NWORK = ITAUP + N
581*
582*              Bidiagonalize R in A
583*              (Workspace: need 4*N, prefer 3*N+2*N*NB)
584*
585               CALL SGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
586     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
587     $                      IERR )
588               NWORK = IE + N
589*
590*              Perform bidiagonal SVD, computing singular values only
591*              (Workspace: need N+BDSPAC)
592*
593               CALL SBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
594     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
595*
596            ELSE IF( WNTQO ) THEN
597*
598*              Path 2 (M much larger than N, JOBZ = 'O')
599*              N left singular vectors to be overwritten on A and
600*              N right singular vectors to be computed in VT
601*
602               IR = 1
603*
604*              WORK(IR) is LDWRKR by N
605*
606               IF( LWORK.GE.LDA*N+N*N+3*N+BDSPAC ) THEN
607                  LDWRKR = LDA
608               ELSE
609                  LDWRKR = ( LWORK-N*N-3*N-BDSPAC ) / N
610               END IF
611               ITAU = IR + LDWRKR*N
612               NWORK = ITAU + N
613*
614*              Compute A=Q*R
615*              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
616*
617               CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
618     $                      LWORK-NWORK+1, IERR )
619*
620*              Copy R to WORK(IR), zeroing out below it
621*
622               CALL SLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
623               CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
624     $                      LDWRKR )
625*
626*              Generate Q in A
627*              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
628*
629               CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
630     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
631               IE = ITAU
632               ITAUQ = IE + N
633               ITAUP = ITAUQ + N
634               NWORK = ITAUP + N
635*
636*              Bidiagonalize R in VT, copying result to WORK(IR)
637*              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
638*
639               CALL SGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
640     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
641     $                      LWORK-NWORK+1, IERR )
642*
643*              WORK(IU) is N by N
644*
645               IU = NWORK
646               NWORK = IU + N*N
647*
648*              Perform bidiagonal SVD, computing left singular vectors
649*              of bidiagonal matrix in WORK(IU) and computing right
650*              singular vectors of bidiagonal matrix in VT
651*              (Workspace: need N+N*N+BDSPAC)
652*
653               CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
654     $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
655     $                      INFO )
656*
657*              Overwrite WORK(IU) by left singular vectors of R
658*              and VT by right singular vectors of R
659*              (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
660*
661               CALL SORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
662     $                      WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
663     $                      LWORK-NWORK+1, IERR )
664               CALL SORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
665     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
666     $                      LWORK-NWORK+1, IERR )
667*
668*              Multiply Q in A by left singular vectors of R in
669*              WORK(IU), storing result in WORK(IR) and copying to A
670*              (Workspace: need 2*N*N, prefer N*N+M*N)
671*
672               DO 10 I = 1, M, LDWRKR
673                  CHUNK = MIN( M-I+1, LDWRKR )
674                  CALL SGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
675     $                        LDA, WORK( IU ), N, ZERO, WORK( IR ),
676     $                        LDWRKR )
677                  CALL SLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
678     $                         A( I, 1 ), LDA )
679   10          CONTINUE
680*
681            ELSE IF( WNTQS ) THEN
682*
683*              Path 3 (M much larger than N, JOBZ='S')
684*              N left singular vectors to be computed in U and
685*              N right singular vectors to be computed in VT
686*
687               IR = 1
688*
689*              WORK(IR) is N by N
690*
691               LDWRKR = N
692               ITAU = IR + LDWRKR*N
693               NWORK = ITAU + N
694*
695*              Compute A=Q*R
696*              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
697*
698               CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
699     $                      LWORK-NWORK+1, IERR )
700*
701*              Copy R to WORK(IR), zeroing out below it
702*
703               CALL SLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
704               CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
705     $                      LDWRKR )
706*
707*              Generate Q in A
708*              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
709*
710               CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
711     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
712               IE = ITAU
713               ITAUQ = IE + N
714               ITAUP = ITAUQ + N
715               NWORK = ITAUP + N
716*
717*              Bidiagonalize R in WORK(IR)
718*              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
719*
720               CALL SGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
721     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
722     $                      LWORK-NWORK+1, IERR )
723*
724*              Perform bidiagonal SVD, computing left singular vectors
725*              of bidiagoal matrix in U and computing right singular
726*              vectors of bidiagonal matrix in VT
727*              (Workspace: need N+BDSPAC)
728*
729               CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
730     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
731     $                      INFO )
732*
733*              Overwrite U by left singular vectors of R and VT
734*              by right singular vectors of R
735*              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
736*
737               CALL SORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
738     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
739     $                      LWORK-NWORK+1, IERR )
740*
741               CALL SORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
742     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
743     $                      LWORK-NWORK+1, IERR )
744*
745*              Multiply Q in A by left singular vectors of R in
746*              WORK(IR), storing result in U
747*              (Workspace: need N*N)
748*
749               CALL SLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
750               CALL SGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),
751     $                     LDWRKR, ZERO, U, LDU )
752*
753            ELSE IF( WNTQA ) THEN
754*
755*              Path 4 (M much larger than N, JOBZ='A')
756*              M left singular vectors to be computed in U and
757*              N right singular vectors to be computed in VT
758*
759               IU = 1
760*
761*              WORK(IU) is N by N
762*
763               LDWRKU = N
764               ITAU = IU + LDWRKU*N
765               NWORK = ITAU + N
766*
767*              Compute A=Q*R, copying result to U
768*              (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
769*
770               CALL SGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
771     $                      LWORK-NWORK+1, IERR )
772               CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
773*
774*              Generate Q in U
775*              (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
776               CALL SORGQR( M, M, N, U, LDU, WORK( ITAU ),
777     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
778*
779*              Produce R in A, zeroing out other entries
780*
781               CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
782               IE = ITAU
783               ITAUQ = IE + N
784               ITAUP = ITAUQ + N
785               NWORK = ITAUP + N
786*
787*              Bidiagonalize R in A
788*              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
789*
790               CALL SGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
791     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
792     $                      IERR )
793*
794*              Perform bidiagonal SVD, computing left singular vectors
795*              of bidiagonal matrix in WORK(IU) and computing right
796*              singular vectors of bidiagonal matrix in VT
797*              (Workspace: need N+N*N+BDSPAC)
798*
799               CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
800     $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
801     $                      INFO )
802*
803*              Overwrite WORK(IU) by left singular vectors of R and VT
804*              by right singular vectors of R
805*              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
806*
807               CALL SORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
808     $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
809     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
810               CALL SORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
811     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
812     $                      LWORK-NWORK+1, IERR )
813*
814*              Multiply Q in U by left singular vectors of R in
815*              WORK(IU), storing result in A
816*              (Workspace: need N*N)
817*
818               CALL SGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),
819     $                     LDWRKU, ZERO, A, LDA )
820*
821*              Copy left singular vectors of A from A to U
822*
823               CALL SLACPY( 'F', M, N, A, LDA, U, LDU )
824*
825            END IF
826*
827         ELSE
828*
829*           M .LT. MNTHR
830*
831*           Path 5 (M at least N, but not much larger)
832*           Reduce to bidiagonal form without QR decomposition
833*
834            IE = 1
835            ITAUQ = IE + N
836            ITAUP = ITAUQ + N
837            NWORK = ITAUP + N
838*
839*           Bidiagonalize A
840*           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
841*
842            CALL SGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
843     $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
844     $                   IERR )
845            IF( WNTQN ) THEN
846*
847*              Perform bidiagonal SVD, only computing singular values
848*              (Workspace: need N+BDSPAC)
849*
850               CALL SBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
851     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
852            ELSE IF( WNTQO ) THEN
853               IU = NWORK
854               IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
855*
856*                 WORK( IU ) is M by N
857*
858                  LDWRKU = M
859                  NWORK = IU + LDWRKU*N
860                  CALL SLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),
861     $                         LDWRKU )
862               ELSE
863*
864*                 WORK( IU ) is N by N
865*
866                  LDWRKU = N
867                  NWORK = IU + LDWRKU*N
868*
869*                 WORK(IR) is LDWRKR by N
870*
871                  IR = NWORK
872                  LDWRKR = ( LWORK-N*N-3*N ) / N
873               END IF
874               NWORK = IU + LDWRKU*N
875*
876*              Perform bidiagonal SVD, computing left singular vectors
877*              of bidiagonal matrix in WORK(IU) and computing right
878*              singular vectors of bidiagonal matrix in VT
879*              (Workspace: need N+N*N+BDSPAC)
880*
881               CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
882     $                      LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
883     $                      IWORK, INFO )
884*
885*              Overwrite VT by right singular vectors of A
886*              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
887*
888               CALL SORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
889     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
890     $                      LWORK-NWORK+1, IERR )
891*
892               IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
893*
894*                 Overwrite WORK(IU) by left singular vectors of A
895*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
896*
897                  CALL SORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
898     $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,
899     $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
900*
901*                 Copy left singular vectors of A from WORK(IU) to A
902*
903                  CALL SLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
904               ELSE
905*
906*                 Generate Q in A
907*                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
908*
909                  CALL SORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
910     $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
911*
912*                 Multiply Q in A by left singular vectors of
913*                 bidiagonal matrix in WORK(IU), storing result in
914*                 WORK(IR) and copying to A
915*                 (Workspace: need 2*N*N, prefer N*N+M*N)
916*
917                  DO 20 I = 1, M, LDWRKR
918                     CHUNK = MIN( M-I+1, LDWRKR )
919                     CALL SGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
920     $                           LDA, WORK( IU ), LDWRKU, ZERO,
921     $                           WORK( IR ), LDWRKR )
922                     CALL SLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
923     $                            A( I, 1 ), LDA )
924   20             CONTINUE
925               END IF
926*
927            ELSE IF( WNTQS ) THEN
928*
929*              Perform bidiagonal SVD, computing left singular vectors
930*              of bidiagonal matrix in U and computing right singular
931*              vectors of bidiagonal matrix in VT
932*              (Workspace: need N+BDSPAC)
933*
934               CALL SLASET( 'F', M, N, ZERO, ZERO, U, LDU )
935               CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
936     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
937     $                      INFO )
938*
939*              Overwrite U by left singular vectors of A and VT
940*              by right singular vectors of A
941*              (Workspace: need 3*N, prefer 2*N+N*NB)
942*
943               CALL SORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
944     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
945     $                      LWORK-NWORK+1, IERR )
946               CALL SORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
947     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
948     $                      LWORK-NWORK+1, IERR )
949            ELSE IF( WNTQA ) THEN
950*
951*              Perform bidiagonal SVD, computing left singular vectors
952*              of bidiagonal matrix in U and computing right singular
953*              vectors of bidiagonal matrix in VT
954*              (Workspace: need N+BDSPAC)
955*
956               CALL SLASET( 'F', M, M, ZERO, ZERO, U, LDU )
957               CALL SBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
958     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
959     $                      INFO )
960*
961*              Set the right corner of U to identity matrix
962*
963               IF( M.GT.N ) THEN
964                  CALL SLASET( 'F', M-N, M-N, ZERO, ONE, U( N+1, N+1 ),
965     $                         LDU )
966               END IF
967*
968*              Overwrite U by left singular vectors of A and VT
969*              by right singular vectors of A
970*              (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB)
971*
972               CALL SORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
973     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
974     $                      LWORK-NWORK+1, IERR )
975               CALL SORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
976     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
977     $                      LWORK-NWORK+1, IERR )
978            END IF
979*
980         END IF
981*
982      ELSE
983*
984*        A has more columns than rows. If A has sufficiently more
985*        columns than rows, first reduce using the LQ decomposition (if
986*        sufficient workspace available)
987*
988         IF( N.GE.MNTHR ) THEN
989*
990            IF( WNTQN ) THEN
991*
992*              Path 1t (N much larger than M, JOBZ='N')
993*              No singular vectors to be computed
994*
995               ITAU = 1
996               NWORK = ITAU + M
997*
998*              Compute A=L*Q
999*              (Workspace: need 2*M, prefer M+M*NB)
1000*
1001               CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1002     $                      LWORK-NWORK+1, IERR )
1003*
1004*              Zero out above L
1005*
1006               CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1007               IE = 1
1008               ITAUQ = IE + M
1009               ITAUP = ITAUQ + M
1010               NWORK = ITAUP + M
1011*
1012*              Bidiagonalize L in A
1013*              (Workspace: need 4*M, prefer 3*M+2*M*NB)
1014*
1015               CALL SGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1016     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1017     $                      IERR )
1018               NWORK = IE + M
1019*
1020*              Perform bidiagonal SVD, computing singular values only
1021*              (Workspace: need M+BDSPAC)
1022*
1023               CALL SBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1024     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1025*
1026            ELSE IF( WNTQO ) THEN
1027*
1028*              Path 2t (N much larger than M, JOBZ='O')
1029*              M right singular vectors to be overwritten on A and
1030*              M left singular vectors to be computed in U
1031*
1032               IVT = 1
1033*
1034*              IVT is M by M
1035*
1036               IL = IVT + M*M
1037               IF( LWORK.GE.M*N+M*M+3*M+BDSPAC ) THEN
1038*
1039*                 WORK(IL) is M by N
1040*
1041                  LDWRKL = M
1042                  CHUNK = N
1043               ELSE
1044                  LDWRKL = M
1045                  CHUNK = ( LWORK-M*M ) / M
1046               END IF
1047               ITAU = IL + LDWRKL*M
1048               NWORK = ITAU + M
1049*
1050*              Compute A=L*Q
1051*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1052*
1053               CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1054     $                      LWORK-NWORK+1, IERR )
1055*
1056*              Copy L to WORK(IL), zeroing about above it
1057*
1058               CALL SLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1059               CALL SLASET( 'U', M-1, M-1, ZERO, ZERO,
1060     $                      WORK( IL+LDWRKL ), LDWRKL )
1061*
1062*              Generate Q in A
1063*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1064*
1065               CALL SORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1066     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1067               IE = ITAU
1068               ITAUQ = IE + M
1069               ITAUP = ITAUQ + M
1070               NWORK = ITAUP + M
1071*
1072*              Bidiagonalize L in WORK(IL)
1073*              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1074*
1075               CALL SGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1076     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1077     $                      LWORK-NWORK+1, IERR )
1078*
1079*              Perform bidiagonal SVD, computing left singular vectors
1080*              of bidiagonal matrix in U, and computing right singular
1081*              vectors of bidiagonal matrix in WORK(IVT)
1082*              (Workspace: need M+M*M+BDSPAC)
1083*
1084               CALL SBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1085     $                      WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
1086     $                      IWORK, INFO )
1087*
1088*              Overwrite U by left singular vectors of L and WORK(IVT)
1089*              by right singular vectors of L
1090*              (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
1091*
1092               CALL SORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1093     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1094     $                      LWORK-NWORK+1, IERR )
1095               CALL SORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1096     $                      WORK( ITAUP ), WORK( IVT ), M,
1097     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1098*
1099*              Multiply right singular vectors of L in WORK(IVT) by Q
1100*              in A, storing result in WORK(IL) and copying to A
1101*              (Workspace: need 2*M*M, prefer M*M+M*N)
1102*
1103               DO 30 I = 1, N, CHUNK
1104                  BLK = MIN( N-I+1, CHUNK )
1105                  CALL SGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,
1106     $                        A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
1107                  CALL SLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
1108     $                         A( 1, I ), LDA )
1109   30          CONTINUE
1110*
1111            ELSE IF( WNTQS ) THEN
1112*
1113*              Path 3t (N much larger than M, JOBZ='S')
1114*              M right singular vectors to be computed in VT and
1115*              M left singular vectors to be computed in U
1116*
1117               IL = 1
1118*
1119*              WORK(IL) is M by M
1120*
1121               LDWRKL = M
1122               ITAU = IL + LDWRKL*M
1123               NWORK = ITAU + M
1124*
1125*              Compute A=L*Q
1126*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1127*
1128               CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1129     $                      LWORK-NWORK+1, IERR )
1130*
1131*              Copy L to WORK(IL), zeroing out above it
1132*
1133               CALL SLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1134               CALL SLASET( 'U', M-1, M-1, ZERO, ZERO,
1135     $                      WORK( IL+LDWRKL ), LDWRKL )
1136*
1137*              Generate Q in A
1138*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1139*
1140               CALL SORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1141     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1142               IE = ITAU
1143               ITAUQ = IE + M
1144               ITAUP = ITAUQ + M
1145               NWORK = ITAUP + M
1146*
1147*              Bidiagonalize L in WORK(IU), copying result to U
1148*              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1149*
1150               CALL SGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1151     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1152     $                      LWORK-NWORK+1, IERR )
1153*
1154*              Perform bidiagonal SVD, computing left singular vectors
1155*              of bidiagonal matrix in U and computing right singular
1156*              vectors of bidiagonal matrix in VT
1157*              (Workspace: need M+BDSPAC)
1158*
1159               CALL SBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
1160     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1161     $                      INFO )
1162*
1163*              Overwrite U by left singular vectors of L and VT
1164*              by right singular vectors of L
1165*              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1166*
1167               CALL SORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1168     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1169     $                      LWORK-NWORK+1, IERR )
1170               CALL SORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1171     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1172     $                      LWORK-NWORK+1, IERR )
1173*
1174*              Multiply right singular vectors of L in WORK(IL) by
1175*              Q in A, storing result in VT
1176*              (Workspace: need M*M)
1177*
1178               CALL SLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1179               CALL SGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,
1180     $                     A, LDA, ZERO, VT, LDVT )
1181*
1182            ELSE IF( WNTQA ) THEN
1183*
1184*              Path 4t (N much larger than M, JOBZ='A')
1185*              N right singular vectors to be computed in VT and
1186*              M left singular vectors to be computed in U
1187*
1188               IVT = 1
1189*
1190*              WORK(IVT) is M by M
1191*
1192               LDWKVT = M
1193               ITAU = IVT + LDWKVT*M
1194               NWORK = ITAU + M
1195*
1196*              Compute A=L*Q, copying result to VT
1197*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1198*
1199               CALL SGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1200     $                      LWORK-NWORK+1, IERR )
1201               CALL SLACPY( 'U', M, N, A, LDA, VT, LDVT )
1202*
1203*              Generate Q in VT
1204*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1205*
1206               CALL SORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1207     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1208*
1209*              Produce L in A, zeroing out other entries
1210*
1211               CALL SLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1212               IE = ITAU
1213               ITAUQ = IE + M
1214               ITAUP = ITAUQ + M
1215               NWORK = ITAUP + M
1216*
1217*              Bidiagonalize L in A
1218*              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1219*
1220               CALL SGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1221     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1222     $                      IERR )
1223*
1224*              Perform bidiagonal SVD, computing left singular vectors
1225*              of bidiagonal matrix in U and computing right singular
1226*              vectors of bidiagonal matrix in WORK(IVT)
1227*              (Workspace: need M+M*M+BDSPAC)
1228*
1229               CALL SBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1230     $                      WORK( IVT ), LDWKVT, DUM, IDUM,
1231     $                      WORK( NWORK ), IWORK, INFO )
1232*
1233*              Overwrite U by left singular vectors of L and WORK(IVT)
1234*              by right singular vectors of L
1235*              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1236*
1237               CALL SORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1238     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1239     $                      LWORK-NWORK+1, IERR )
1240               CALL SORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
1241     $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
1242     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1243*
1244*              Multiply right singular vectors of L in WORK(IVT) by
1245*              Q in VT, storing result in A
1246*              (Workspace: need M*M)
1247*
1248               CALL SGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,
1249     $                     VT, LDVT, ZERO, A, LDA )
1250*
1251*              Copy right singular vectors of A from A to VT
1252*
1253               CALL SLACPY( 'F', M, N, A, LDA, VT, LDVT )
1254*
1255            END IF
1256*
1257         ELSE
1258*
1259*           N .LT. MNTHR
1260*
1261*           Path 5t (N greater than M, but not much larger)
1262*           Reduce to bidiagonal form without LQ decomposition
1263*
1264            IE = 1
1265            ITAUQ = IE + M
1266            ITAUP = ITAUQ + M
1267            NWORK = ITAUP + M
1268*
1269*           Bidiagonalize A
1270*           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
1271*
1272            CALL SGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1273     $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1274     $                   IERR )
1275            IF( WNTQN ) THEN
1276*
1277*              Perform bidiagonal SVD, only computing singular values
1278*              (Workspace: need M+BDSPAC)
1279*
1280               CALL SBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1281     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1282            ELSE IF( WNTQO ) THEN
1283               LDWKVT = M
1284               IVT = NWORK
1285               IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
1286*
1287*                 WORK( IVT ) is M by N
1288*
1289                  CALL SLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),
1290     $                         LDWKVT )
1291                  NWORK = IVT + LDWKVT*N
1292               ELSE
1293*
1294*                 WORK( IVT ) is M by M
1295*
1296                  NWORK = IVT + LDWKVT*M
1297                  IL = NWORK
1298*
1299*                 WORK(IL) is M by CHUNK
1300*
1301                  CHUNK = ( LWORK-M*M-3*M ) / M
1302               END IF
1303*
1304*              Perform bidiagonal SVD, computing left singular vectors
1305*              of bidiagonal matrix in U and computing right singular
1306*              vectors of bidiagonal matrix in WORK(IVT)
1307*              (Workspace: need M*M+BDSPAC)
1308*
1309               CALL SBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,
1310     $                      WORK( IVT ), LDWKVT, DUM, IDUM,
1311     $                      WORK( NWORK ), IWORK, INFO )
1312*
1313*              Overwrite U by left singular vectors of A
1314*              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1315*
1316               CALL SORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1317     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1318     $                      LWORK-NWORK+1, IERR )
1319*
1320               IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
1321*
1322*                 Overwrite WORK(IVT) by left singular vectors of A
1323*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1324*
1325                  CALL SORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1326     $                         WORK( ITAUP ), WORK( IVT ), LDWKVT,
1327     $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1328*
1329*                 Copy right singular vectors of A from WORK(IVT) to A
1330*
1331                  CALL SLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
1332               ELSE
1333*
1334*                 Generate P**T in A
1335*                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1336*
1337                  CALL SORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1338     $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1339*
1340*                 Multiply Q in A by right singular vectors of
1341*                 bidiagonal matrix in WORK(IVT), storing result in
1342*                 WORK(IL) and copying to A
1343*                 (Workspace: need 2*M*M, prefer M*M+M*N)
1344*
1345                  DO 40 I = 1, N, CHUNK
1346                     BLK = MIN( N-I+1, CHUNK )
1347                     CALL SGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),
1348     $                           LDWKVT, A( 1, I ), LDA, ZERO,
1349     $                           WORK( IL ), M )
1350                     CALL SLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ),
1351     $                            LDA )
1352   40             CONTINUE
1353               END IF
1354            ELSE IF( WNTQS ) THEN
1355*
1356*              Perform bidiagonal SVD, computing left singular vectors
1357*              of bidiagonal matrix in U and computing right singular
1358*              vectors of bidiagonal matrix in VT
1359*              (Workspace: need M+BDSPAC)
1360*
1361               CALL SLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
1362               CALL SBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1363     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1364     $                      INFO )
1365*
1366*              Overwrite U by left singular vectors of A and VT
1367*              by right singular vectors of A
1368*              (Workspace: need 3*M, prefer 2*M+M*NB)
1369*
1370               CALL SORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1371     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1372     $                      LWORK-NWORK+1, IERR )
1373               CALL SORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1374     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1375     $                      LWORK-NWORK+1, IERR )
1376            ELSE IF( WNTQA ) THEN
1377*
1378*              Perform bidiagonal SVD, computing left singular vectors
1379*              of bidiagonal matrix in U and computing right singular
1380*              vectors of bidiagonal matrix in VT
1381*              (Workspace: need M+BDSPAC)
1382*
1383               CALL SLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
1384               CALL SBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1385     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1386     $                      INFO )
1387*
1388*              Set the right corner of VT to identity matrix
1389*
1390               IF( N.GT.M ) THEN
1391                  CALL SLASET( 'F', N-M, N-M, ZERO, ONE, VT( M+1, M+1 ),
1392     $                         LDVT )
1393               END IF
1394*
1395*              Overwrite U by left singular vectors of A and VT
1396*              by right singular vectors of A
1397*              (Workspace: need 2*M+N, prefer 2*M+N*NB)
1398*
1399               CALL SORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1400     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1401     $                      LWORK-NWORK+1, IERR )
1402               CALL SORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1403     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1404     $                      LWORK-NWORK+1, IERR )
1405            END IF
1406*
1407         END IF
1408*
1409      END IF
1410*
1411*     Undo scaling if necessary
1412*
1413      IF( ISCL.EQ.1 ) THEN
1414         IF( ANRM.GT.BIGNUM )
1415     $      CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
1416     $                   IERR )
1417         IF( ANRM.LT.SMLNUM )
1418     $      CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
1419     $                   IERR )
1420      END IF
1421*
1422*     Return optimal workspace in WORK(1)
1423*
1424      WORK( 1 ) = MAXWRK
1425*
1426      RETURN
1427*
1428*     End of SGESDD
1429*
1430      END
1431