LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dgesdd.f
Go to the documentation of this file.
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, 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 * 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; 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 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 JOBZ = 'N',
173 *> LWORK >= 3*min(M,N) + max(max(M,N),7*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)*(6+4*min(M,N))+max(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: DBDSDC 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 2013
206 *
207 *> \ingroup doubleGEsing
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 dgesdd( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
217  $ lwork, iwork, info )
218 *
219 * -- LAPACK driver routine (version 3.5.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 2013
223 *
224 * .. Scalar Arguments ..
225  CHARACTER jobz
226  INTEGER info, lda, ldu, ldvt, lwork, m, n
227 * ..
228 * .. Array Arguments ..
229  INTEGER iwork( * )
230  DOUBLE PRECISION a( lda, * ), s( * ), u( ldu, * ),
231  $ vt( ldvt, * ), work( * )
232 * ..
233 *
234 * =====================================================================
235 *
236 * .. Parameters ..
237  DOUBLE PRECISION zero, one
238  parameter( zero = 0.0d0, one = 1.0d0 )
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  DOUBLE PRECISION anrm, bignum, eps, smlnum
247 * ..
248 * .. Local Arrays ..
249  INTEGER idum( 1 )
250  DOUBLE PRECISION dum( 1 )
251 * ..
252 * .. External Subroutines ..
253  EXTERNAL dbdsdc, dgebrd, dgelqf, dgemm, dgeqrf, dlacpy,
255  $ xerbla
256 * ..
257 * .. External Functions ..
258  LOGICAL lsame
259  INTEGER ilaenv
260  DOUBLE PRECISION dlamch, dlange
261  EXTERNAL dlamch, dlange, ilaenv, lsame
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 DBDSDC
309 *
310  mnthr = int( minmn*11.0d0 / 6.0d0 )
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, 'DGEQRF', ' ', m, n, -1,
322  $ -1 )
323  wrkbl = max( wrkbl, 3*n+2*n*
324  $ ilaenv( 1, 'DGEBRD', ' ', 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, 'DGEQRF', ' ', m, n, -1, -1 )
332  wrkbl = max( wrkbl, n+n*ilaenv( 1, 'DORGQR', ' ', m,
333  $ n, n, -1 ) )
334  wrkbl = max( wrkbl, 3*n+2*n*
335  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
336  wrkbl = max( wrkbl, 3*n+n*
337  $ ilaenv( 1, 'DORMBR', 'QLN', n, n, n, -1 ) )
338  wrkbl = max( wrkbl, 3*n+n*
339  $ ilaenv( 1, 'DORMBR', '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, 'DGEQRF', ' ', m, n, -1, -1 )
348  wrkbl = max( wrkbl, n+n*ilaenv( 1, 'DORGQR', ' ', m,
349  $ n, n, -1 ) )
350  wrkbl = max( wrkbl, 3*n+2*n*
351  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
352  wrkbl = max( wrkbl, 3*n+n*
353  $ ilaenv( 1, 'DORMBR', 'QLN', n, n, n, -1 ) )
354  wrkbl = max( wrkbl, 3*n+n*
355  $ ilaenv( 1, 'DORMBR', '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, 'DGEQRF', ' ', m, n, -1, -1 )
364  wrkbl = max( wrkbl, n+m*ilaenv( 1, 'DORGQR', ' ', m,
365  $ m, n, -1 ) )
366  wrkbl = max( wrkbl, 3*n+2*n*
367  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
368  wrkbl = max( wrkbl, 3*n+n*
369  $ ilaenv( 1, 'DORMBR', 'QLN', n, n, n, -1 ) )
370  wrkbl = max( wrkbl, 3*n+n*
371  $ ilaenv( 1, 'DORMBR', '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, 'DGEBRD', ' ', 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, 'DORMBR', 'QLN', m, n, n, -1 ) )
388  wrkbl = max( wrkbl, 3*n+n*
389  $ ilaenv( 1, 'DORMBR', '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, 'DORMBR', 'QLN', m, n, n, -1 ) )
396  wrkbl = max( wrkbl, 3*n+n*
397  $ ilaenv( 1, 'DORMBR', '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, 'DORMBR', 'QLN', m, m, n, -1 ) )
403  wrkbl = max( wrkbl, 3*n+n*
404  $ ilaenv( 1, 'DORMBR', '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 DBDSDC
412 *
413  mnthr = int( minmn*11.0d0 / 6.0d0 )
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, 'DGELQF', ' ', m, n, -1,
425  $ -1 )
426  wrkbl = max( wrkbl, 3*m+2*m*
427  $ ilaenv( 1, 'DGEBRD', ' ', 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, 'DGELQF', ' ', m, n, -1, -1 )
435  wrkbl = max( wrkbl, m+m*ilaenv( 1, 'DORGLQ', ' ', m,
436  $ n, m, -1 ) )
437  wrkbl = max( wrkbl, 3*m+2*m*
438  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
439  wrkbl = max( wrkbl, 3*m+m*
440  $ ilaenv( 1, 'DORMBR', 'QLN', m, m, m, -1 ) )
441  wrkbl = max( wrkbl, 3*m+m*
442  $ ilaenv( 1, 'DORMBR', '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, 'DGELQF', ' ', m, n, -1, -1 )
451  wrkbl = max( wrkbl, m+m*ilaenv( 1, 'DORGLQ', ' ', m,
452  $ n, m, -1 ) )
453  wrkbl = max( wrkbl, 3*m+2*m*
454  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
455  wrkbl = max( wrkbl, 3*m+m*
456  $ ilaenv( 1, 'DORMBR', 'QLN', m, m, m, -1 ) )
457  wrkbl = max( wrkbl, 3*m+m*
458  $ ilaenv( 1, 'DORMBR', '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, 'DGELQF', ' ', m, n, -1, -1 )
467  wrkbl = max( wrkbl, m+n*ilaenv( 1, 'DORGLQ', ' ', n,
468  $ n, m, -1 ) )
469  wrkbl = max( wrkbl, 3*m+2*m*
470  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
471  wrkbl = max( wrkbl, 3*m+m*
472  $ ilaenv( 1, 'DORMBR', 'QLN', m, m, m, -1 ) )
473  wrkbl = max( wrkbl, 3*m+m*
474  $ ilaenv( 1, 'DORMBR', '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, 'DGEBRD', ' ', 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, 'DORMBR', 'QLN', m, m, n, -1 ) )
491  wrkbl = max( wrkbl, 3*m+m*
492  $ ilaenv( 1, 'DORMBR', '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, 'DORMBR', 'QLN', m, m, n, -1 ) )
499  wrkbl = max( wrkbl, 3*m+m*
500  $ ilaenv( 1, 'DORMBR', '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, 'DORMBR', 'QLN', m, m, n, -1 ) )
506  wrkbl = max( wrkbl, 3*m+m*
507  $ ilaenv( 1, 'DORMBR', '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( 'DGESDD', -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 = dlamch( 'P' )
537  smlnum = sqrt( dlamch( 'S' ) ) / eps
538  bignum = one / smlnum
539 *
540 * Scale A if max element outside range [SMLNUM,BIGNUM]
541 *
542  anrm = dlange( 'M', m, n, a, lda, dum )
543  iscl = 0
544  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
545  iscl = 1
546  CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
547  ELSE IF( anrm.GT.bignum ) THEN
548  iscl = 1
549  CALL dlascl( '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 dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
572  $ lwork-nwork+1, ierr )
573 *
574 * Zero out below R
575 *
576  CALL dlaset( '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 dgebrd( 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 dbdsdc( '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 dgeqrf( 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 dlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
623  CALL dlaset( '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 dorgqr( 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 dgebrd( 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 dbdsdc( '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 dormbr( 'Q', 'L', 'N', n, n, n, work( ir ), ldwrkr,
662  $ work( itauq ), work( iu ), n, work( nwork ),
663  $ lwork-nwork+1, ierr )
664  CALL dormbr( '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 dgemm( 'N', 'N', chunk, n, n, one, a( i, 1 ),
675  $ lda, work( iu ), n, zero, work( ir ),
676  $ ldwrkr )
677  CALL dlacpy( '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 dgeqrf( 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 dlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
704  CALL dlaset( '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 dorgqr( 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 dgebrd( 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 dbdsdc( '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 dormbr( 'Q', 'L', 'N', n, n, n, work( ir ), ldwrkr,
738  $ work( itauq ), u, ldu, work( nwork ),
739  $ lwork-nwork+1, ierr )
740 *
741  CALL dormbr( '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 dlacpy( 'F', n, n, u, ldu, work( ir ), ldwrkr )
750  CALL dgemm( '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 dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
771  $ lwork-nwork+1, ierr )
772  CALL dlacpy( '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 dorgqr( 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 dlaset( '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 dgebrd( 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 dbdsdc( '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 dormbr( 'Q', 'L', 'N', n, n, n, a, lda,
808  $ work( itauq ), work( iu ), ldwrku,
809  $ work( nwork ), lwork-nwork+1, ierr )
810  CALL dormbr( '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 dgemm( '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 dlacpy( '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 dgebrd( 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 dbdsdc( '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 dlaset( '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 dbdsdc( '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 dormbr( '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 dormbr( '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 dlacpy( '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 dorgbr( '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 dgemm( 'N', 'N', chunk, n, n, one, a( i, 1 ),
920  $ lda, work( iu ), ldwrku, zero,
921  $ work( ir ), ldwrkr )
922  CALL dlacpy( '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 dlaset( 'F', m, n, zero, zero, u, ldu )
935  CALL dbdsdc( '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 dormbr( 'Q', 'L', 'N', m, n, n, a, lda,
944  $ work( itauq ), u, ldu, work( nwork ),
945  $ lwork-nwork+1, ierr )
946  CALL dormbr( '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 dlaset( 'F', m, m, zero, zero, u, ldu )
957  CALL dbdsdc( '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 dlaset( '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 dormbr( 'Q', 'L', 'N', m, m, n, a, lda,
973  $ work( itauq ), u, ldu, work( nwork ),
974  $ lwork-nwork+1, ierr )
975  CALL dormbr( '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 dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1002  $ lwork-nwork+1, ierr )
1003 *
1004 * Zero out above L
1005 *
1006  CALL dlaset( '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 dgebrd( 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 dbdsdc( '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 dgelqf( 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 dlacpy( 'L', m, m, a, lda, work( il ), ldwrkl )
1059  CALL dlaset( '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 dorglq( 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 dgebrd( 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 dbdsdc( '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 dormbr( 'Q', 'L', 'N', m, m, m, work( il ), ldwrkl,
1093  $ work( itauq ), u, ldu, work( nwork ),
1094  $ lwork-nwork+1, ierr )
1095  CALL dormbr( '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 dgemm( 'N', 'N', m, blk, m, one, work( ivt ), m,
1106  $ a( 1, i ), lda, zero, work( il ), ldwrkl )
1107  CALL dlacpy( '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 dgelqf( 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 dlacpy( 'L', m, m, a, lda, work( il ), ldwrkl )
1134  CALL dlaset( '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 dorglq( 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 dgebrd( 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 dbdsdc( '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 dormbr( 'Q', 'L', 'N', m, m, m, work( il ), ldwrkl,
1168  $ work( itauq ), u, ldu, work( nwork ),
1169  $ lwork-nwork+1, ierr )
1170  CALL dormbr( '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 dlacpy( 'F', m, m, vt, ldvt, work( il ), ldwrkl )
1179  CALL dgemm( '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 dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1200  $ lwork-nwork+1, ierr )
1201  CALL dlacpy( '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 dorglq( 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 dlaset( '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 dgebrd( 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 dbdsdc( '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 dormbr( 'Q', 'L', 'N', m, m, m, a, lda,
1238  $ work( itauq ), u, ldu, work( nwork ),
1239  $ lwork-nwork+1, ierr )
1240  CALL dormbr( '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 dgemm( '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 dlacpy( '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 dgebrd( 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 dbdsdc( '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 dlaset( '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 dbdsdc( '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 dormbr( '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 dormbr( '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 dlacpy( '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 dorgbr( '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 dgemm( 'N', 'N', m, blk, m, one, work( ivt ),
1348  $ ldwkvt, a( 1, i ), lda, zero,
1349  $ work( il ), m )
1350  CALL dlacpy( '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 dlaset( 'F', m, n, zero, zero, vt, ldvt )
1362  CALL dbdsdc( '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 dormbr( 'Q', 'L', 'N', m, m, n, a, lda,
1371  $ work( itauq ), u, ldu, work( nwork ),
1372  $ lwork-nwork+1, ierr )
1373  CALL dormbr( '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 dlaset( 'F', n, n, zero, zero, vt, ldvt )
1384  CALL dbdsdc( '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 dlaset( '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 dormbr( 'Q', 'L', 'N', m, m, n, a, lda,
1400  $ work( itauq ), u, ldu, work( nwork ),
1401  $ lwork-nwork+1, ierr )
1402  CALL dormbr( '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 dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
1416  $ ierr )
1417  IF( anrm.LT.smlnum )
1418  $ CALL dlascl( '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 DGESDD
1429 *
1430  END