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