KTH framework for Nek5000 toolboxes; testing version  0.0.1
dgesvd.f
Go to the documentation of this file.
1  SUBROUTINE dgesvd( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
2  $ WORK, LWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.0) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 * Courant Institute, Argonne National Lab, and Rice University
7 * October 31, 1999
8 *
9 * .. Scalar Arguments ..
10  CHARACTER JOBU, JOBVT
11  INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
12 * ..
13 * .. Array Arguments ..
14  DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
15  $ vt( ldvt, * ), work( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DGESVD computes the singular value decomposition (SVD) of a real
22 * M-by-N matrix A, optionally computing the left and/or right singular
23 * vectors. The SVD is written
24 *
25 * A = U * SIGMA * transpose(V)
26 *
27 * where SIGMA is an M-by-N matrix which is zero except for its
28 * min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
29 * V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
30 * are the singular values of A; they are real and non-negative, and
31 * are returned in descending order. The first min(m,n) columns of
32 * U and V are the left and right singular vectors of A.
33 *
34 * Note that the routine returns V**T, not V.
35 *
36 * Arguments
37 * =========
38 *
39 * JOBU (input) CHARACTER*1
40 * Specifies options for computing all or part of the matrix U:
41 * = 'A': all M columns of U are returned in array U:
42 * = 'S': the first min(m,n) columns of U (the left singular
43 * vectors) are returned in the array U;
44 * = 'O': the first min(m,n) columns of U (the left singular
45 * vectors) are overwritten on the array A;
46 * = 'N': no columns of U (no left singular vectors) are
47 * computed.
48 *
49 * JOBVT (input) CHARACTER*1
50 * Specifies options for computing all or part of the matrix
51 * V**T:
52 * = 'A': all N rows of V**T are returned in the array VT;
53 * = 'S': the first min(m,n) rows of V**T (the right singular
54 * vectors) are returned in the array VT;
55 * = 'O': the first min(m,n) rows of V**T (the right singular
56 * vectors) are overwritten on the array A;
57 * = 'N': no rows of V**T (no right singular vectors) are
58 * computed.
59 *
60 * JOBVT and JOBU cannot both be 'O'.
61 *
62 * M (input) INTEGER
63 * The number of rows of the input matrix A. M >= 0.
64 *
65 * N (input) INTEGER
66 * The number of columns of the input matrix A. N >= 0.
67 *
68 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
69 * On entry, the M-by-N matrix A.
70 * On exit,
71 * if JOBU = 'O', A is overwritten with the first min(m,n)
72 * columns of U (the left singular vectors,
73 * stored columnwise);
74 * if JOBVT = 'O', A is overwritten with the first min(m,n)
75 * rows of V**T (the right singular vectors,
76 * stored rowwise);
77 * if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
78 * are destroyed.
79 *
80 * LDA (input) INTEGER
81 * The leading dimension of the array A. LDA >= max(1,M).
82 *
83 * S (output) DOUBLE PRECISION array, dimension (min(M,N))
84 * The singular values of A, sorted so that S(i) >= S(i+1).
85 *
86 * U (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
87 * (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
88 * If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
89 * if JOBU = 'S', U contains the first min(m,n) columns of U
90 * (the left singular vectors, stored columnwise);
91 * if JOBU = 'N' or 'O', U is not referenced.
92 *
93 * LDU (input) INTEGER
94 * The leading dimension of the array U. LDU >= 1; if
95 * JOBU = 'S' or 'A', LDU >= M.
96 *
97 * VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
98 * If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
99 * V**T;
100 * if JOBVT = 'S', VT contains the first min(m,n) rows of
101 * V**T (the right singular vectors, stored rowwise);
102 * if JOBVT = 'N' or 'O', VT is not referenced.
103 *
104 * LDVT (input) INTEGER
105 * The leading dimension of the array VT. LDVT >= 1; if
106 * JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
107 *
108 * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
109 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
110 * if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
111 * superdiagonal elements of an upper bidiagonal matrix B
112 * whose diagonal is in S (not necessarily sorted). B
113 * satisfies A = U * B * VT, so it has the same singular values
114 * as A, and singular vectors related by U and VT.
115 *
116 * LWORK (input) INTEGER
117 * The dimension of the array WORK. LWORK >= 1.
118 * LWORK >= MAX(3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
119 * For good performance, LWORK should generally be larger.
120 *
121 * If LWORK = -1, then a workspace query is assumed; the routine
122 * only calculates the optimal size of the WORK array, returns
123 * this value as the first entry of the WORK array, and no error
124 * message related to LWORK is issued by XERBLA.
125 *
126 * INFO (output) INTEGER
127 * = 0: successful exit.
128 * < 0: if INFO = -i, the i-th argument had an illegal value.
129 * > 0: if DBDSQR did not converge, INFO specifies how many
130 * superdiagonals of an intermediate bidiagonal form B
131 * did not converge to zero. See the description of WORK
132 * above for details.
133 *
134 * =====================================================================
135 *
136 * .. Parameters ..
137  DOUBLE PRECISION ZERO, ONE
138  parameter( zero = 0.0d0, one = 1.0d0 )
139 * ..
140 * .. Local Scalars ..
141  LOGICAL LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS,
142  $ wntva, wntvas, wntvn, wntvo, wntvs
143  INTEGER BDSPAC, BLK, CHUNK, I, IE, IERR, IR, ISCL,
144  $ itau, itaup, itauq, iu, iwork, ldwrkr, ldwrku,
145  $ maxwrk, minmn, minwrk, mnthr, ncu, ncvt, nru,
146  $ nrvt, wrkbl
147  DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
148 * ..
149 * .. Local Arrays ..
150  DOUBLE PRECISION DUM( 1 )
151 * ..
152 * .. External Subroutines ..
153  EXTERNAL dbdsqr, dgebrd, dgelqf, dgemm, dgeqrf, dlacpy,
154  $ dlascl, dlaset, dorgbr, dorglq, dorgqr, dormbr,
155  $ xerbla
156 * ..
157 * .. External Functions ..
158  LOGICAL LSAME
159  INTEGER ILAENV
160  DOUBLE PRECISION DLAMCH, DLANGE
161  EXTERNAL lsame, ilaenv, dlamch, dlange
162 * ..
163 * .. Intrinsic Functions ..
164  INTRINSIC max, min, sqrt
165 * ..
166 * .. Executable Statements ..
167 *
168 * Test the input arguments
169 *
170  info = 0
171  minmn = min( m, n )
172  mnthr = ilaenv( 6, 'DGESVD', jobu // jobvt, m, n, 0, 0 )
173  wntua = lsame( jobu, 'A' )
174  wntus = lsame( jobu, 'S' )
175  wntuas = wntua .OR. wntus
176  wntuo = lsame( jobu, 'O' )
177  wntun = lsame( jobu, 'N' )
178  wntva = lsame( jobvt, 'A' )
179  wntvs = lsame( jobvt, 'S' )
180  wntvas = wntva .OR. wntvs
181  wntvo = lsame( jobvt, 'O' )
182  wntvn = lsame( jobvt, 'N' )
183  minwrk = 1
184  lquery = ( lwork.EQ.-1 )
185 *
186  IF( .NOT.( wntua .OR. wntus .OR. wntuo .OR. wntun ) ) THEN
187  info = -1
188  ELSE IF( .NOT.( wntva .OR. wntvs .OR. wntvo .OR. wntvn ) .OR.
189  $ ( wntvo .AND. wntuo ) ) THEN
190  info = -2
191  ELSE IF( m.LT.0 ) THEN
192  info = -3
193  ELSE IF( n.LT.0 ) THEN
194  info = -4
195  ELSE IF( lda.LT.max( 1, m ) ) THEN
196  info = -6
197  ELSE IF( ldu.LT.1 .OR. ( wntuas .AND. ldu.LT.m ) ) THEN
198  info = -9
199  ELSE IF( ldvt.LT.1 .OR. ( wntva .AND. ldvt.LT.n ) .OR.
200  $ ( wntvs .AND. ldvt.LT.minmn ) ) THEN
201  info = -11
202  END IF
203 *
204 * Compute workspace
205 * (Note: Comments in the code beginning "Workspace:" describe the
206 * minimal amount of workspace needed at that point in the code,
207 * as well as the preferred amount for good performance.
208 * NB refers to the optimal block size for the immediately
209 * following subroutine, as returned by ILAENV.)
210 *
211  IF( info.EQ.0 .AND. ( lwork.GE.1 .OR. lquery ) .AND. m.GT.0 .AND.
212  $ n.GT.0 ) THEN
213  IF( m.GE.n ) THEN
214 *
215 * Compute space needed for DBDSQR
216 *
217  bdspac = 5*n
218  IF( m.GE.mnthr ) THEN
219  IF( wntun ) THEN
220 *
221 * Path 1 (M much larger than N, JOBU='N')
222 *
223  maxwrk = n + n*ilaenv( 1, 'DGEQRF', ' ', m, n, -1,
224  $ -1 )
225  maxwrk = max( maxwrk, 3*n+2*n*
226  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
227  IF( wntvo .OR. wntvas )
228  $ maxwrk = max( maxwrk, 3*n+( n-1 )*
229  $ ilaenv( 1, 'DORGBR', 'P', n, n, n, -1 ) )
230  maxwrk = max( maxwrk, bdspac )
231  minwrk = max( 4*n, bdspac )
232  maxwrk = max( maxwrk, minwrk )
233  ELSE IF( wntuo .AND. wntvn ) THEN
234 *
235 * Path 2 (M much larger than N, JOBU='O', JOBVT='N')
236 *
237  wrkbl = n + n*ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
238  wrkbl = max( wrkbl, n+n*ilaenv( 1, 'DORGQR', ' ', m,
239  $ n, n, -1 ) )
240  wrkbl = max( wrkbl, 3*n+2*n*
241  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
242  wrkbl = max( wrkbl, 3*n+n*
243  $ ilaenv( 1, 'DORGBR', 'Q', n, n, n, -1 ) )
244  wrkbl = max( wrkbl, bdspac )
245  maxwrk = max( n*n+wrkbl, n*n+m*n+n )
246  minwrk = max( 3*n+m, bdspac )
247  maxwrk = max( maxwrk, minwrk )
248  ELSE IF( wntuo .AND. wntvas ) THEN
249 *
250 * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
251 * 'A')
252 *
253  wrkbl = n + n*ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
254  wrkbl = max( wrkbl, n+n*ilaenv( 1, 'DORGQR', ' ', m,
255  $ n, n, -1 ) )
256  wrkbl = max( wrkbl, 3*n+2*n*
257  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
258  wrkbl = max( wrkbl, 3*n+n*
259  $ ilaenv( 1, 'DORGBR', 'Q', n, n, n, -1 ) )
260  wrkbl = max( wrkbl, 3*n+( n-1 )*
261  $ ilaenv( 1, 'DORGBR', 'P', n, n, n, -1 ) )
262  wrkbl = max( wrkbl, bdspac )
263  maxwrk = max( n*n+wrkbl, n*n+m*n+n )
264  minwrk = max( 3*n+m, bdspac )
265  maxwrk = max( maxwrk, minwrk )
266  ELSE IF( wntus .AND. wntvn ) THEN
267 *
268 * Path 4 (M much larger than N, JOBU='S', JOBVT='N')
269 *
270  wrkbl = n + n*ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
271  wrkbl = max( wrkbl, n+n*ilaenv( 1, 'DORGQR', ' ', m,
272  $ n, n, -1 ) )
273  wrkbl = max( wrkbl, 3*n+2*n*
274  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
275  wrkbl = max( wrkbl, 3*n+n*
276  $ ilaenv( 1, 'DORGBR', 'Q', n, n, n, -1 ) )
277  wrkbl = max( wrkbl, bdspac )
278  maxwrk = n*n + wrkbl
279  minwrk = max( 3*n+m, bdspac )
280  maxwrk = max( maxwrk, minwrk )
281  ELSE IF( wntus .AND. wntvo ) THEN
282 *
283 * Path 5 (M much larger than N, JOBU='S', JOBVT='O')
284 *
285  wrkbl = n + n*ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
286  wrkbl = max( wrkbl, n+n*ilaenv( 1, 'DORGQR', ' ', m,
287  $ n, n, -1 ) )
288  wrkbl = max( wrkbl, 3*n+2*n*
289  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
290  wrkbl = max( wrkbl, 3*n+n*
291  $ ilaenv( 1, 'DORGBR', 'Q', n, n, n, -1 ) )
292  wrkbl = max( wrkbl, 3*n+( n-1 )*
293  $ ilaenv( 1, 'DORGBR', 'P', n, n, n, -1 ) )
294  wrkbl = max( wrkbl, bdspac )
295  maxwrk = 2*n*n + wrkbl
296  minwrk = max( 3*n+m, bdspac )
297  maxwrk = max( maxwrk, minwrk )
298  ELSE IF( wntus .AND. wntvas ) THEN
299 *
300 * Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
301 * 'A')
302 *
303  wrkbl = n + n*ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
304  wrkbl = max( wrkbl, n+n*ilaenv( 1, 'DORGQR', ' ', m,
305  $ n, n, -1 ) )
306  wrkbl = max( wrkbl, 3*n+2*n*
307  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
308  wrkbl = max( wrkbl, 3*n+n*
309  $ ilaenv( 1, 'DORGBR', 'Q', n, n, n, -1 ) )
310  wrkbl = max( wrkbl, 3*n+( n-1 )*
311  $ ilaenv( 1, 'DORGBR', 'P', n, n, n, -1 ) )
312  wrkbl = max( wrkbl, bdspac )
313  maxwrk = n*n + wrkbl
314  minwrk = max( 3*n+m, bdspac )
315  maxwrk = max( maxwrk, minwrk )
316  ELSE IF( wntua .AND. wntvn ) THEN
317 *
318 * Path 7 (M much larger than N, JOBU='A', JOBVT='N')
319 *
320  wrkbl = n + n*ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
321  wrkbl = max( wrkbl, n+m*ilaenv( 1, 'DORGQR', ' ', m,
322  $ m, n, -1 ) )
323  wrkbl = max( wrkbl, 3*n+2*n*
324  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
325  wrkbl = max( wrkbl, 3*n+n*
326  $ ilaenv( 1, 'DORGBR', 'Q', n, n, n, -1 ) )
327  wrkbl = max( wrkbl, bdspac )
328  maxwrk = n*n + wrkbl
329  minwrk = max( 3*n+m, bdspac )
330  maxwrk = max( maxwrk, minwrk )
331  ELSE IF( wntua .AND. wntvo ) THEN
332 *
333 * Path 8 (M much larger than N, JOBU='A', JOBVT='O')
334 *
335  wrkbl = n + n*ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
336  wrkbl = max( wrkbl, n+m*ilaenv( 1, 'DORGQR', ' ', m,
337  $ m, n, -1 ) )
338  wrkbl = max( wrkbl, 3*n+2*n*
339  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
340  wrkbl = max( wrkbl, 3*n+n*
341  $ ilaenv( 1, 'DORGBR', 'Q', n, n, n, -1 ) )
342  wrkbl = max( wrkbl, 3*n+( n-1 )*
343  $ ilaenv( 1, 'DORGBR', 'P', n, n, n, -1 ) )
344  wrkbl = max( wrkbl, bdspac )
345  maxwrk = 2*n*n + wrkbl
346  minwrk = max( 3*n+m, bdspac )
347  maxwrk = max( maxwrk, minwrk )
348  ELSE IF( wntua .AND. wntvas ) THEN
349 *
350 * Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
351 * 'A')
352 *
353  wrkbl = n + n*ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
354  wrkbl = max( wrkbl, n+m*ilaenv( 1, 'DORGQR', ' ', m,
355  $ m, n, -1 ) )
356  wrkbl = max( wrkbl, 3*n+2*n*
357  $ ilaenv( 1, 'DGEBRD', ' ', n, n, -1, -1 ) )
358  wrkbl = max( wrkbl, 3*n+n*
359  $ ilaenv( 1, 'DORGBR', 'Q', n, n, n, -1 ) )
360  wrkbl = max( wrkbl, 3*n+( n-1 )*
361  $ ilaenv( 1, 'DORGBR', 'P', n, n, n, -1 ) )
362  wrkbl = max( wrkbl, bdspac )
363  maxwrk = n*n + wrkbl
364  minwrk = max( 3*n+m, bdspac )
365  maxwrk = max( maxwrk, minwrk )
366  END IF
367  ELSE
368 *
369 * Path 10 (M at least N, but not much larger)
370 *
371  maxwrk = 3*n + ( m+n )*ilaenv( 1, 'DGEBRD', ' ', m, n,
372  $ -1, -1 )
373  IF( wntus .OR. wntuo )
374  $ maxwrk = max( maxwrk, 3*n+n*
375  $ ilaenv( 1, 'DORGBR', 'Q', m, n, n, -1 ) )
376  IF( wntua )
377  $ maxwrk = max( maxwrk, 3*n+m*
378  $ ilaenv( 1, 'DORGBR', 'Q', m, m, n, -1 ) )
379  IF( .NOT.wntvn )
380  $ maxwrk = max( maxwrk, 3*n+( n-1 )*
381  $ ilaenv( 1, 'DORGBR', 'P', n, n, n, -1 ) )
382  maxwrk = max( maxwrk, bdspac )
383  minwrk = max( 3*n+m, bdspac )
384  maxwrk = max( maxwrk, minwrk )
385  END IF
386  ELSE
387 *
388 * Compute space needed for DBDSQR
389 *
390  bdspac = 5*m
391  IF( n.GE.mnthr ) THEN
392  IF( wntvn ) THEN
393 *
394 * Path 1t(N much larger than M, JOBVT='N')
395 *
396  maxwrk = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1,
397  $ -1 )
398  maxwrk = max( maxwrk, 3*m+2*m*
399  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
400  IF( wntuo .OR. wntuas )
401  $ maxwrk = max( maxwrk, 3*m+m*
402  $ ilaenv( 1, 'DORGBR', 'Q', m, m, m, -1 ) )
403  maxwrk = max( maxwrk, bdspac )
404  minwrk = max( 4*m, bdspac )
405  maxwrk = max( maxwrk, minwrk )
406  ELSE IF( wntvo .AND. wntun ) THEN
407 *
408 * Path 2t(N much larger than M, JOBU='N', JOBVT='O')
409 *
410  wrkbl = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
411  wrkbl = max( wrkbl, m+m*ilaenv( 1, 'DORGLQ', ' ', m,
412  $ n, m, -1 ) )
413  wrkbl = max( wrkbl, 3*m+2*m*
414  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
415  wrkbl = max( wrkbl, 3*m+( m-1 )*
416  $ ilaenv( 1, 'DORGBR', 'P', m, m, m, -1 ) )
417  wrkbl = max( wrkbl, bdspac )
418  maxwrk = max( m*m+wrkbl, m*m+m*n+m )
419  minwrk = max( 3*m+n, bdspac )
420  maxwrk = max( maxwrk, minwrk )
421  ELSE IF( wntvo .AND. wntuas ) THEN
422 *
423 * Path 3t(N much larger than M, JOBU='S' or 'A',
424 * JOBVT='O')
425 *
426  wrkbl = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
427  wrkbl = max( wrkbl, m+m*ilaenv( 1, 'DORGLQ', ' ', m,
428  $ n, m, -1 ) )
429  wrkbl = max( wrkbl, 3*m+2*m*
430  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
431  wrkbl = max( wrkbl, 3*m+( m-1 )*
432  $ ilaenv( 1, 'DORGBR', 'P', m, m, m, -1 ) )
433  wrkbl = max( wrkbl, 3*m+m*
434  $ ilaenv( 1, 'DORGBR', 'Q', m, m, m, -1 ) )
435  wrkbl = max( wrkbl, bdspac )
436  maxwrk = max( m*m+wrkbl, m*m+m*n+m )
437  minwrk = max( 3*m+n, bdspac )
438  maxwrk = max( maxwrk, minwrk )
439  ELSE IF( wntvs .AND. wntun ) THEN
440 *
441 * Path 4t(N much larger than M, JOBU='N', JOBVT='S')
442 *
443  wrkbl = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
444  wrkbl = max( wrkbl, m+m*ilaenv( 1, 'DORGLQ', ' ', m,
445  $ n, m, -1 ) )
446  wrkbl = max( wrkbl, 3*m+2*m*
447  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
448  wrkbl = max( wrkbl, 3*m+( m-1 )*
449  $ ilaenv( 1, 'DORGBR', 'P', m, m, m, -1 ) )
450  wrkbl = max( wrkbl, bdspac )
451  maxwrk = m*m + wrkbl
452  minwrk = max( 3*m+n, bdspac )
453  maxwrk = max( maxwrk, minwrk )
454  ELSE IF( wntvs .AND. wntuo ) THEN
455 *
456 * Path 5t(N much larger than M, JOBU='O', JOBVT='S')
457 *
458  wrkbl = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
459  wrkbl = max( wrkbl, m+m*ilaenv( 1, 'DORGLQ', ' ', m,
460  $ n, m, -1 ) )
461  wrkbl = max( wrkbl, 3*m+2*m*
462  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
463  wrkbl = max( wrkbl, 3*m+( m-1 )*
464  $ ilaenv( 1, 'DORGBR', 'P', m, m, m, -1 ) )
465  wrkbl = max( wrkbl, 3*m+m*
466  $ ilaenv( 1, 'DORGBR', 'Q', m, m, m, -1 ) )
467  wrkbl = max( wrkbl, bdspac )
468  maxwrk = 2*m*m + wrkbl
469  minwrk = max( 3*m+n, bdspac )
470  maxwrk = max( maxwrk, minwrk )
471  ELSE IF( wntvs .AND. wntuas ) THEN
472 *
473 * Path 6t(N much larger than M, JOBU='S' or 'A',
474 * JOBVT='S')
475 *
476  wrkbl = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
477  wrkbl = max( wrkbl, m+m*ilaenv( 1, 'DORGLQ', ' ', m,
478  $ n, m, -1 ) )
479  wrkbl = max( wrkbl, 3*m+2*m*
480  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
481  wrkbl = max( wrkbl, 3*m+( m-1 )*
482  $ ilaenv( 1, 'DORGBR', 'P', m, m, m, -1 ) )
483  wrkbl = max( wrkbl, 3*m+m*
484  $ ilaenv( 1, 'DORGBR', 'Q', m, m, m, -1 ) )
485  wrkbl = max( wrkbl, bdspac )
486  maxwrk = m*m + wrkbl
487  minwrk = max( 3*m+n, bdspac )
488  maxwrk = max( maxwrk, minwrk )
489  ELSE IF( wntva .AND. wntun ) THEN
490 *
491 * Path 7t(N much larger than M, JOBU='N', JOBVT='A')
492 *
493  wrkbl = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
494  wrkbl = max( wrkbl, m+n*ilaenv( 1, 'DORGLQ', ' ', n,
495  $ n, m, -1 ) )
496  wrkbl = max( wrkbl, 3*m+2*m*
497  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
498  wrkbl = max( wrkbl, 3*m+( m-1 )*
499  $ ilaenv( 1, 'DORGBR', 'P', m, m, m, -1 ) )
500  wrkbl = max( wrkbl, bdspac )
501  maxwrk = m*m + wrkbl
502  minwrk = max( 3*m+n, bdspac )
503  maxwrk = max( maxwrk, minwrk )
504  ELSE IF( wntva .AND. wntuo ) THEN
505 *
506 * Path 8t(N much larger than M, JOBU='O', JOBVT='A')
507 *
508  wrkbl = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
509  wrkbl = max( wrkbl, m+n*ilaenv( 1, 'DORGLQ', ' ', n,
510  $ n, m, -1 ) )
511  wrkbl = max( wrkbl, 3*m+2*m*
512  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
513  wrkbl = max( wrkbl, 3*m+( m-1 )*
514  $ ilaenv( 1, 'DORGBR', 'P', m, m, m, -1 ) )
515  wrkbl = max( wrkbl, 3*m+m*
516  $ ilaenv( 1, 'DORGBR', 'Q', m, m, m, -1 ) )
517  wrkbl = max( wrkbl, bdspac )
518  maxwrk = 2*m*m + wrkbl
519  minwrk = max( 3*m+n, bdspac )
520  maxwrk = max( maxwrk, minwrk )
521  ELSE IF( wntva .AND. wntuas ) THEN
522 *
523 * Path 9t(N much larger than M, JOBU='S' or 'A',
524 * JOBVT='A')
525 *
526  wrkbl = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
527  wrkbl = max( wrkbl, m+n*ilaenv( 1, 'DORGLQ', ' ', n,
528  $ n, m, -1 ) )
529  wrkbl = max( wrkbl, 3*m+2*m*
530  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
531  wrkbl = max( wrkbl, 3*m+( m-1 )*
532  $ ilaenv( 1, 'DORGBR', 'P', m, m, m, -1 ) )
533  wrkbl = max( wrkbl, 3*m+m*
534  $ ilaenv( 1, 'DORGBR', 'Q', m, m, m, -1 ) )
535  wrkbl = max( wrkbl, bdspac )
536  maxwrk = m*m + wrkbl
537  minwrk = max( 3*m+n, bdspac )
538  maxwrk = max( maxwrk, minwrk )
539  END IF
540  ELSE
541 *
542 * Path 10t(N greater than M, but not much larger)
543 *
544  maxwrk = 3*m + ( m+n )*ilaenv( 1, 'DGEBRD', ' ', m, n,
545  $ -1, -1 )
546  IF( wntvs .OR. wntvo )
547  $ maxwrk = max( maxwrk, 3*m+m*
548  $ ilaenv( 1, 'DORGBR', 'P', m, n, m, -1 ) )
549  IF( wntva )
550  $ maxwrk = max( maxwrk, 3*m+n*
551  $ ilaenv( 1, 'DORGBR', 'P', n, n, m, -1 ) )
552  IF( .NOT.wntun )
553  $ maxwrk = max( maxwrk, 3*m+( m-1 )*
554  $ ilaenv( 1, 'DORGBR', 'Q', m, m, m, -1 ) )
555  maxwrk = max( maxwrk, bdspac )
556  minwrk = max( 3*m+n, bdspac )
557  maxwrk = max( maxwrk, minwrk )
558  END IF
559  END IF
560  work( 1 ) = maxwrk
561  END IF
562 *
563  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
564  info = -13
565  END IF
566  IF( info.NE.0 ) THEN
567  CALL xerbla( 'DGESVD', -info )
568  RETURN
569  ELSE IF( lquery ) THEN
570  RETURN
571  END IF
572 *
573 * Quick return if possible
574 *
575  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
576  IF( lwork.GE.1 )
577  $ work( 1 ) = one
578  RETURN
579  END IF
580 *
581 * Get machine constants
582 *
583  eps = dlamch( 'P' )
584  smlnum = sqrt( dlamch( 'S' ) ) / eps
585  bignum = one / smlnum
586 *
587 * Scale A if max element outside range [SMLNUM,BIGNUM]
588 *
589  anrm = dlange( 'M', m, n, a, lda, dum )
590  iscl = 0
591  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
592  iscl = 1
593  CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
594  ELSE IF( anrm.GT.bignum ) THEN
595  iscl = 1
596  CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
597  END IF
598 *
599  IF( m.GE.n ) THEN
600 *
601 * A has at least as many rows as columns. If A has sufficiently
602 * more rows than columns, first reduce using the QR
603 * decomposition (if sufficient workspace available)
604 *
605  IF( m.GE.mnthr ) THEN
606 *
607  IF( wntun ) THEN
608 *
609 * Path 1 (M much larger than N, JOBU='N')
610 * No left singular vectors to be computed
611 *
612  itau = 1
613  iwork = itau + n
614 *
615 * Compute A=Q*R
616 * (Workspace: need 2*N, prefer N+N*NB)
617 *
618  CALL dgeqrf( m, n, a, lda, work( itau ), work( iwork ),
619  $ lwork-iwork+1, ierr )
620 *
621 * Zero out below R
622 *
623  CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
624  ie = 1
625  itauq = ie + n
626  itaup = itauq + n
627  iwork = itaup + n
628 *
629 * Bidiagonalize R in A
630 * (Workspace: need 4*N, prefer 3*N+2*N*NB)
631 *
632  CALL dgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
633  $ work( itaup ), work( iwork ), lwork-iwork+1,
634  $ ierr )
635  ncvt = 0
636  IF( wntvo .OR. wntvas ) THEN
637 *
638 * If right singular vectors desired, generate P'.
639 * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
640 *
641  CALL dorgbr( 'P', n, n, n, a, lda, work( itaup ),
642  $ work( iwork ), lwork-iwork+1, ierr )
643  ncvt = n
644  END IF
645  iwork = ie + n
646 *
647 * Perform bidiagonal QR iteration, computing right
648 * singular vectors of A in A if desired
649 * (Workspace: need BDSPAC)
650 *
651  CALL dbdsqr( 'U', n, ncvt, 0, 0, s, work( ie ), a, lda,
652  $ dum, 1, dum, 1, work( iwork ), info )
653 *
654 * If right singular vectors desired in VT, copy them there
655 *
656  IF( wntvas )
657  $ CALL dlacpy( 'F', n, n, a, lda, vt, ldvt )
658 *
659  ELSE IF( wntuo .AND. wntvn ) THEN
660 *
661 * Path 2 (M much larger than N, JOBU='O', JOBVT='N')
662 * N left singular vectors to be overwritten on A and
663 * no right singular vectors to be computed
664 *
665  IF( lwork.GE.n*n+max( 4*n, bdspac ) ) THEN
666 *
667 * Sufficient workspace for a fast algorithm
668 *
669  ir = 1
670  IF( lwork.GE.max( wrkbl, lda*n+n )+lda*n ) THEN
671 *
672 * WORK(IU) is LDA by N, WORK(IR) is LDA by N
673 *
674  ldwrku = lda
675  ldwrkr = lda
676  ELSE IF( lwork.GE.max( wrkbl, lda*n+n )+n*n ) THEN
677 *
678 * WORK(IU) is LDA by N, WORK(IR) is N by N
679 *
680  ldwrku = lda
681  ldwrkr = n
682  ELSE
683 *
684 * WORK(IU) is LDWRKU by N, WORK(IR) is N by N
685 *
686  ldwrku = ( lwork-n*n-n ) / n
687  ldwrkr = n
688  END IF
689  itau = ir + ldwrkr*n
690  iwork = itau + n
691 *
692 * Compute A=Q*R
693 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
694 *
695  CALL dgeqrf( m, n, a, lda, work( itau ),
696  $ work( iwork ), lwork-iwork+1, ierr )
697 *
698 * Copy R to WORK(IR) and zero out below it
699 *
700  CALL dlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
701  CALL dlaset( 'L', n-1, n-1, zero, zero, work( ir+1 ),
702  $ ldwrkr )
703 *
704 * Generate Q in A
705 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
706 *
707  CALL dorgqr( m, n, n, a, lda, work( itau ),
708  $ work( iwork ), lwork-iwork+1, ierr )
709  ie = itau
710  itauq = ie + n
711  itaup = itauq + n
712  iwork = itaup + n
713 *
714 * Bidiagonalize R in WORK(IR)
715 * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
716 *
717  CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
718  $ work( itauq ), work( itaup ),
719  $ work( iwork ), lwork-iwork+1, ierr )
720 *
721 * Generate left vectors bidiagonalizing R
722 * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
723 *
724  CALL dorgbr( 'Q', n, n, n, work( ir ), ldwrkr,
725  $ work( itauq ), work( iwork ),
726  $ lwork-iwork+1, ierr )
727  iwork = ie + n
728 *
729 * Perform bidiagonal QR iteration, computing left
730 * singular vectors of R in WORK(IR)
731 * (Workspace: need N*N+BDSPAC)
732 *
733  CALL dbdsqr( 'U', n, 0, n, 0, s, work( ie ), dum, 1,
734  $ work( ir ), ldwrkr, dum, 1,
735  $ work( iwork ), info )
736  iu = ie + n
737 *
738 * Multiply Q in A by left singular vectors of R in
739 * WORK(IR), storing result in WORK(IU) and copying to A
740 * (Workspace: need N*N+2*N, prefer N*N+M*N+N)
741 *
742  DO 10 i = 1, m, ldwrku
743  chunk = min( m-i+1, ldwrku )
744  CALL dgemm( 'N', 'N', chunk, n, n, one, a( i, 1 ),
745  $ lda, work( ir ), ldwrkr, zero,
746  $ work( iu ), ldwrku )
747  CALL dlacpy( 'F', chunk, n, work( iu ), ldwrku,
748  $ a( i, 1 ), lda )
749  10 CONTINUE
750 *
751  ELSE
752 *
753 * Insufficient workspace for a fast algorithm
754 *
755  ie = 1
756  itauq = ie + n
757  itaup = itauq + n
758  iwork = itaup + n
759 *
760 * Bidiagonalize A
761 * (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
762 *
763  CALL dgebrd( m, n, a, lda, s, work( ie ),
764  $ work( itauq ), work( itaup ),
765  $ work( iwork ), lwork-iwork+1, ierr )
766 *
767 * Generate left vectors bidiagonalizing A
768 * (Workspace: need 4*N, prefer 3*N+N*NB)
769 *
770  CALL dorgbr( 'Q', m, n, n, a, lda, work( itauq ),
771  $ work( iwork ), lwork-iwork+1, ierr )
772  iwork = ie + n
773 *
774 * Perform bidiagonal QR iteration, computing left
775 * singular vectors of A in A
776 * (Workspace: need BDSPAC)
777 *
778  CALL dbdsqr( 'U', n, 0, m, 0, s, work( ie ), dum, 1,
779  $ a, lda, dum, 1, work( iwork ), info )
780 *
781  END IF
782 *
783  ELSE IF( wntuo .AND. wntvas ) THEN
784 *
785 * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
786 * N left singular vectors to be overwritten on A and
787 * N right singular vectors to be computed in VT
788 *
789  IF( lwork.GE.n*n+max( 4*n, bdspac ) ) THEN
790 *
791 * Sufficient workspace for a fast algorithm
792 *
793  ir = 1
794  IF( lwork.GE.max( wrkbl, lda*n+n )+lda*n ) THEN
795 *
796 * WORK(IU) is LDA by N and WORK(IR) is LDA by N
797 *
798  ldwrku = lda
799  ldwrkr = lda
800  ELSE IF( lwork.GE.max( wrkbl, lda*n+n )+n*n ) THEN
801 *
802 * WORK(IU) is LDA by N and WORK(IR) is N by N
803 *
804  ldwrku = lda
805  ldwrkr = n
806  ELSE
807 *
808 * WORK(IU) is LDWRKU by N and WORK(IR) is N by N
809 *
810  ldwrku = ( lwork-n*n-n ) / n
811  ldwrkr = n
812  END IF
813  itau = ir + ldwrkr*n
814  iwork = itau + n
815 *
816 * Compute A=Q*R
817 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
818 *
819  CALL dgeqrf( m, n, a, lda, work( itau ),
820  $ work( iwork ), lwork-iwork+1, ierr )
821 *
822 * Copy R to VT, zeroing out below it
823 *
824  CALL dlacpy( 'U', n, n, a, lda, vt, ldvt )
825  CALL dlaset( 'L', n-1, n-1, zero, zero, vt( 2, 1 ),
826  $ ldvt )
827 *
828 * Generate Q in A
829 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
830 *
831  CALL dorgqr( m, n, n, a, lda, work( itau ),
832  $ work( iwork ), lwork-iwork+1, ierr )
833  ie = itau
834  itauq = ie + n
835  itaup = itauq + n
836  iwork = itaup + n
837 *
838 * Bidiagonalize R in VT, copying result to WORK(IR)
839 * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
840 *
841  CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
842  $ work( itauq ), work( itaup ),
843  $ work( iwork ), lwork-iwork+1, ierr )
844  CALL dlacpy( 'L', n, n, vt, ldvt, work( ir ), ldwrkr )
845 *
846 * Generate left vectors bidiagonalizing R in WORK(IR)
847 * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
848 *
849  CALL dorgbr( 'Q', n, n, n, work( ir ), ldwrkr,
850  $ work( itauq ), work( iwork ),
851  $ lwork-iwork+1, ierr )
852 *
853 * Generate right vectors bidiagonalizing R in VT
854 * (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB)
855 *
856  CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
857  $ work( iwork ), lwork-iwork+1, ierr )
858  iwork = ie + n
859 *
860 * Perform bidiagonal QR iteration, computing left
861 * singular vectors of R in WORK(IR) and computing right
862 * singular vectors of R in VT
863 * (Workspace: need N*N+BDSPAC)
864 *
865  CALL dbdsqr( 'U', n, n, n, 0, s, work( ie ), vt, ldvt,
866  $ work( ir ), ldwrkr, dum, 1,
867  $ work( iwork ), info )
868  iu = ie + n
869 *
870 * Multiply Q in A by left singular vectors of R in
871 * WORK(IR), storing result in WORK(IU) and copying to A
872 * (Workspace: need N*N+2*N, prefer N*N+M*N+N)
873 *
874  DO 20 i = 1, m, ldwrku
875  chunk = min( m-i+1, ldwrku )
876  CALL dgemm( 'N', 'N', chunk, n, n, one, a( i, 1 ),
877  $ lda, work( ir ), ldwrkr, zero,
878  $ work( iu ), ldwrku )
879  CALL dlacpy( 'F', chunk, n, work( iu ), ldwrku,
880  $ a( i, 1 ), lda )
881  20 CONTINUE
882 *
883  ELSE
884 *
885 * Insufficient workspace for a fast algorithm
886 *
887  itau = 1
888  iwork = itau + n
889 *
890 * Compute A=Q*R
891 * (Workspace: need 2*N, prefer N+N*NB)
892 *
893  CALL dgeqrf( m, n, a, lda, work( itau ),
894  $ work( iwork ), lwork-iwork+1, ierr )
895 *
896 * Copy R to VT, zeroing out below it
897 *
898  CALL dlacpy( 'U', n, n, a, lda, vt, ldvt )
899  CALL dlaset( 'L', n-1, n-1, zero, zero, vt( 2, 1 ),
900  $ ldvt )
901 *
902 * Generate Q in A
903 * (Workspace: need 2*N, prefer N+N*NB)
904 *
905  CALL dorgqr( m, n, n, a, lda, work( itau ),
906  $ work( iwork ), lwork-iwork+1, ierr )
907  ie = itau
908  itauq = ie + n
909  itaup = itauq + n
910  iwork = itaup + n
911 *
912 * Bidiagonalize R in VT
913 * (Workspace: need 4*N, prefer 3*N+2*N*NB)
914 *
915  CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
916  $ work( itauq ), work( itaup ),
917  $ work( iwork ), lwork-iwork+1, ierr )
918 *
919 * Multiply Q in A by left vectors bidiagonalizing R
920 * (Workspace: need 3*N+M, prefer 3*N+M*NB)
921 *
922  CALL dormbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
923  $ work( itauq ), a, lda, work( iwork ),
924  $ lwork-iwork+1, ierr )
925 *
926 * Generate right vectors bidiagonalizing R in VT
927 * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
928 *
929  CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
930  $ work( iwork ), lwork-iwork+1, ierr )
931  iwork = ie + n
932 *
933 * Perform bidiagonal QR iteration, computing left
934 * singular vectors of A in A and computing right
935 * singular vectors of A in VT
936 * (Workspace: need BDSPAC)
937 *
938  CALL dbdsqr( 'U', n, n, m, 0, s, work( ie ), vt, ldvt,
939  $ a, lda, dum, 1, work( iwork ), info )
940 *
941  END IF
942 *
943  ELSE IF( wntus ) THEN
944 *
945  IF( wntvn ) THEN
946 *
947 * Path 4 (M much larger than N, JOBU='S', JOBVT='N')
948 * N left singular vectors to be computed in U and
949 * no right singular vectors to be computed
950 *
951  IF( lwork.GE.n*n+max( 4*n, bdspac ) ) THEN
952 *
953 * Sufficient workspace for a fast algorithm
954 *
955  ir = 1
956  IF( lwork.GE.wrkbl+lda*n ) THEN
957 *
958 * WORK(IR) is LDA by N
959 *
960  ldwrkr = lda
961  ELSE
962 *
963 * WORK(IR) is N by N
964 *
965  ldwrkr = n
966  END IF
967  itau = ir + ldwrkr*n
968  iwork = itau + n
969 *
970 * Compute A=Q*R
971 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
972 *
973  CALL dgeqrf( m, n, a, lda, work( itau ),
974  $ work( iwork ), lwork-iwork+1, ierr )
975 *
976 * Copy R to WORK(IR), zeroing out below it
977 *
978  CALL dlacpy( 'U', n, n, a, lda, work( ir ),
979  $ ldwrkr )
980  CALL dlaset( 'L', n-1, n-1, zero, zero,
981  $ work( ir+1 ), ldwrkr )
982 *
983 * Generate Q in A
984 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
985 *
986  CALL dorgqr( m, n, n, a, lda, work( itau ),
987  $ work( iwork ), lwork-iwork+1, ierr )
988  ie = itau
989  itauq = ie + n
990  itaup = itauq + n
991  iwork = itaup + n
992 *
993 * Bidiagonalize R in WORK(IR)
994 * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
995 *
996  CALL dgebrd( n, n, work( ir ), ldwrkr, s,
997  $ work( ie ), work( itauq ),
998  $ work( itaup ), work( iwork ),
999  $ lwork-iwork+1, ierr )
1000 *
1001 * Generate left vectors bidiagonalizing R in WORK(IR)
1002 * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1003 *
1004  CALL dorgbr( 'Q', n, n, n, work( ir ), ldwrkr,
1005  $ work( itauq ), work( iwork ),
1006  $ lwork-iwork+1, ierr )
1007  iwork = ie + n
1008 *
1009 * Perform bidiagonal QR iteration, computing left
1010 * singular vectors of R in WORK(IR)
1011 * (Workspace: need N*N+BDSPAC)
1012 *
1013  CALL dbdsqr( 'U', n, 0, n, 0, s, work( ie ), dum,
1014  $ 1, work( ir ), ldwrkr, dum, 1,
1015  $ work( iwork ), info )
1016 *
1017 * Multiply Q in A by left singular vectors of R in
1018 * WORK(IR), storing result in U
1019 * (Workspace: need N*N)
1020 *
1021  CALL dgemm( 'N', 'N', m, n, n, one, a, lda,
1022  $ work( ir ), ldwrkr, zero, u, ldu )
1023 *
1024  ELSE
1025 *
1026 * Insufficient workspace for a fast algorithm
1027 *
1028  itau = 1
1029  iwork = itau + n
1030 *
1031 * Compute A=Q*R, copying result to U
1032 * (Workspace: need 2*N, prefer N+N*NB)
1033 *
1034  CALL dgeqrf( m, n, a, lda, work( itau ),
1035  $ work( iwork ), lwork-iwork+1, ierr )
1036  CALL dlacpy( 'L', m, n, a, lda, u, ldu )
1037 *
1038 * Generate Q in U
1039 * (Workspace: need 2*N, prefer N+N*NB)
1040 *
1041  CALL dorgqr( m, n, n, u, ldu, work( itau ),
1042  $ work( iwork ), lwork-iwork+1, ierr )
1043  ie = itau
1044  itauq = ie + n
1045  itaup = itauq + n
1046  iwork = itaup + n
1047 *
1048 * Zero out below R in A
1049 *
1050  CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ),
1051  $ lda )
1052 *
1053 * Bidiagonalize R in A
1054 * (Workspace: need 4*N, prefer 3*N+2*N*NB)
1055 *
1056  CALL dgebrd( n, n, a, lda, s, work( ie ),
1057  $ work( itauq ), work( itaup ),
1058  $ work( iwork ), lwork-iwork+1, ierr )
1059 *
1060 * Multiply Q in U by left vectors bidiagonalizing R
1061 * (Workspace: need 3*N+M, prefer 3*N+M*NB)
1062 *
1063  CALL dormbr( 'Q', 'R', 'N', m, n, n, a, lda,
1064  $ work( itauq ), u, ldu, work( iwork ),
1065  $ lwork-iwork+1, ierr )
1066  iwork = ie + n
1067 *
1068 * Perform bidiagonal QR iteration, computing left
1069 * singular vectors of A in U
1070 * (Workspace: need BDSPAC)
1071 *
1072  CALL dbdsqr( 'U', n, 0, m, 0, s, work( ie ), dum,
1073  $ 1, u, ldu, dum, 1, work( iwork ),
1074  $ info )
1075 *
1076  END IF
1077 *
1078  ELSE IF( wntvo ) THEN
1079 *
1080 * Path 5 (M much larger than N, JOBU='S', JOBVT='O')
1081 * N left singular vectors to be computed in U and
1082 * N right singular vectors to be overwritten on A
1083 *
1084  IF( lwork.GE.2*n*n+max( 4*n, bdspac ) ) THEN
1085 *
1086 * Sufficient workspace for a fast algorithm
1087 *
1088  iu = 1
1089  IF( lwork.GE.wrkbl+2*lda*n ) THEN
1090 *
1091 * WORK(IU) is LDA by N and WORK(IR) is LDA by N
1092 *
1093  ldwrku = lda
1094  ir = iu + ldwrku*n
1095  ldwrkr = lda
1096  ELSE IF( lwork.GE.wrkbl+( lda+n )*n ) THEN
1097 *
1098 * WORK(IU) is LDA by N and WORK(IR) is N by N
1099 *
1100  ldwrku = lda
1101  ir = iu + ldwrku*n
1102  ldwrkr = n
1103  ELSE
1104 *
1105 * WORK(IU) is N by N and WORK(IR) is N by N
1106 *
1107  ldwrku = n
1108  ir = iu + ldwrku*n
1109  ldwrkr = n
1110  END IF
1111  itau = ir + ldwrkr*n
1112  iwork = itau + n
1113 *
1114 * Compute A=Q*R
1115 * (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1116 *
1117  CALL dgeqrf( m, n, a, lda, work( itau ),
1118  $ work( iwork ), lwork-iwork+1, ierr )
1119 *
1120 * Copy R to WORK(IU), zeroing out below it
1121 *
1122  CALL dlacpy( 'U', n, n, a, lda, work( iu ),
1123  $ ldwrku )
1124  CALL dlaset( 'L', n-1, n-1, zero, zero,
1125  $ work( iu+1 ), ldwrku )
1126 *
1127 * Generate Q in A
1128 * (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1129 *
1130  CALL dorgqr( m, n, n, a, lda, work( itau ),
1131  $ work( iwork ), lwork-iwork+1, ierr )
1132  ie = itau
1133  itauq = ie + n
1134  itaup = itauq + n
1135  iwork = itaup + n
1136 *
1137 * Bidiagonalize R in WORK(IU), copying result to
1138 * WORK(IR)
1139 * (Workspace: need 2*N*N+4*N,
1140 * prefer 2*N*N+3*N+2*N*NB)
1141 *
1142  CALL dgebrd( n, n, work( iu ), ldwrku, s,
1143  $ work( ie ), work( itauq ),
1144  $ work( itaup ), work( iwork ),
1145  $ lwork-iwork+1, ierr )
1146  CALL dlacpy( 'U', n, n, work( iu ), ldwrku,
1147  $ work( ir ), ldwrkr )
1148 *
1149 * Generate left bidiagonalizing vectors in WORK(IU)
1150 * (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
1151 *
1152  CALL dorgbr( 'Q', n, n, n, work( iu ), ldwrku,
1153  $ work( itauq ), work( iwork ),
1154  $ lwork-iwork+1, ierr )
1155 *
1156 * Generate right bidiagonalizing vectors in WORK(IR)
1157 * (Workspace: need 2*N*N+4*N-1,
1158 * prefer 2*N*N+3*N+(N-1)*NB)
1159 *
1160  CALL dorgbr( 'P', n, n, n, work( ir ), ldwrkr,
1161  $ work( itaup ), work( iwork ),
1162  $ lwork-iwork+1, ierr )
1163  iwork = ie + n
1164 *
1165 * Perform bidiagonal QR iteration, computing left
1166 * singular vectors of R in WORK(IU) and computing
1167 * right singular vectors of R in WORK(IR)
1168 * (Workspace: need 2*N*N+BDSPAC)
1169 *
1170  CALL dbdsqr( 'U', n, n, n, 0, s, work( ie ),
1171  $ work( ir ), ldwrkr, work( iu ),
1172  $ ldwrku, dum, 1, work( iwork ), info )
1173 *
1174 * Multiply Q in A by left singular vectors of R in
1175 * WORK(IU), storing result in U
1176 * (Workspace: need N*N)
1177 *
1178  CALL dgemm( 'N', 'N', m, n, n, one, a, lda,
1179  $ work( iu ), ldwrku, zero, u, ldu )
1180 *
1181 * Copy right singular vectors of R to A
1182 * (Workspace: need N*N)
1183 *
1184  CALL dlacpy( 'F', n, n, work( ir ), ldwrkr, a,
1185  $ lda )
1186 *
1187  ELSE
1188 *
1189 * Insufficient workspace for a fast algorithm
1190 *
1191  itau = 1
1192  iwork = itau + n
1193 *
1194 * Compute A=Q*R, copying result to U
1195 * (Workspace: need 2*N, prefer N+N*NB)
1196 *
1197  CALL dgeqrf( m, n, a, lda, work( itau ),
1198  $ work( iwork ), lwork-iwork+1, ierr )
1199  CALL dlacpy( 'L', m, n, a, lda, u, ldu )
1200 *
1201 * Generate Q in U
1202 * (Workspace: need 2*N, prefer N+N*NB)
1203 *
1204  CALL dorgqr( m, n, n, u, ldu, work( itau ),
1205  $ work( iwork ), lwork-iwork+1, ierr )
1206  ie = itau
1207  itauq = ie + n
1208  itaup = itauq + n
1209  iwork = itaup + n
1210 *
1211 * Zero out below R in A
1212 *
1213  CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ),
1214  $ lda )
1215 *
1216 * Bidiagonalize R in A
1217 * (Workspace: need 4*N, prefer 3*N+2*N*NB)
1218 *
1219  CALL dgebrd( n, n, a, lda, s, work( ie ),
1220  $ work( itauq ), work( itaup ),
1221  $ work( iwork ), lwork-iwork+1, ierr )
1222 *
1223 * Multiply Q in U by left vectors bidiagonalizing R
1224 * (Workspace: need 3*N+M, prefer 3*N+M*NB)
1225 *
1226  CALL dormbr( 'Q', 'R', 'N', m, n, n, a, lda,
1227  $ work( itauq ), u, ldu, work( iwork ),
1228  $ lwork-iwork+1, ierr )
1229 *
1230 * Generate right vectors bidiagonalizing R in A
1231 * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1232 *
1233  CALL dorgbr( 'P', n, n, n, a, lda, work( itaup ),
1234  $ work( iwork ), lwork-iwork+1, ierr )
1235  iwork = ie + n
1236 *
1237 * Perform bidiagonal QR iteration, computing left
1238 * singular vectors of A in U and computing right
1239 * singular vectors of A in A
1240 * (Workspace: need BDSPAC)
1241 *
1242  CALL dbdsqr( 'U', n, n, m, 0, s, work( ie ), a,
1243  $ lda, u, ldu, dum, 1, work( iwork ),
1244  $ info )
1245 *
1246  END IF
1247 *
1248  ELSE IF( wntvas ) THEN
1249 *
1250 * Path 6 (M much larger than N, JOBU='S', JOBVT='S'
1251 * or 'A')
1252 * N left singular vectors to be computed in U and
1253 * N right singular vectors to be computed in VT
1254 *
1255  IF( lwork.GE.n*n+max( 4*n, bdspac ) ) THEN
1256 *
1257 * Sufficient workspace for a fast algorithm
1258 *
1259  iu = 1
1260  IF( lwork.GE.wrkbl+lda*n ) THEN
1261 *
1262 * WORK(IU) is LDA by N
1263 *
1264  ldwrku = lda
1265  ELSE
1266 *
1267 * WORK(IU) is N by N
1268 *
1269  ldwrku = n
1270  END IF
1271  itau = iu + ldwrku*n
1272  iwork = itau + n
1273 *
1274 * Compute A=Q*R
1275 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1276 *
1277  CALL dgeqrf( m, n, a, lda, work( itau ),
1278  $ work( iwork ), lwork-iwork+1, ierr )
1279 *
1280 * Copy R to WORK(IU), zeroing out below it
1281 *
1282  CALL dlacpy( 'U', n, n, a, lda, work( iu ),
1283  $ ldwrku )
1284  CALL dlaset( 'L', n-1, n-1, zero, zero,
1285  $ work( iu+1 ), ldwrku )
1286 *
1287 * Generate Q in A
1288 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1289 *
1290  CALL dorgqr( m, n, n, a, lda, work( itau ),
1291  $ work( iwork ), lwork-iwork+1, ierr )
1292  ie = itau
1293  itauq = ie + n
1294  itaup = itauq + n
1295  iwork = itaup + n
1296 *
1297 * Bidiagonalize R in WORK(IU), copying result to VT
1298 * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1299 *
1300  CALL dgebrd( n, n, work( iu ), ldwrku, s,
1301  $ work( ie ), work( itauq ),
1302  $ work( itaup ), work( iwork ),
1303  $ lwork-iwork+1, ierr )
1304  CALL dlacpy( 'U', n, n, work( iu ), ldwrku, vt,
1305  $ ldvt )
1306 *
1307 * Generate left bidiagonalizing vectors in WORK(IU)
1308 * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1309 *
1310  CALL dorgbr( 'Q', n, n, n, work( iu ), ldwrku,
1311  $ work( itauq ), work( iwork ),
1312  $ lwork-iwork+1, ierr )
1313 *
1314 * Generate right bidiagonalizing vectors in VT
1315 * (Workspace: need N*N+4*N-1,
1316 * prefer N*N+3*N+(N-1)*NB)
1317 *
1318  CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1319  $ work( iwork ), lwork-iwork+1, ierr )
1320  iwork = ie + n
1321 *
1322 * Perform bidiagonal QR iteration, computing left
1323 * singular vectors of R in WORK(IU) and computing
1324 * right singular vectors of R in VT
1325 * (Workspace: need N*N+BDSPAC)
1326 *
1327  CALL dbdsqr( 'U', n, n, n, 0, s, work( ie ), vt,
1328  $ ldvt, work( iu ), ldwrku, dum, 1,
1329  $ work( iwork ), info )
1330 *
1331 * Multiply Q in A by left singular vectors of R in
1332 * WORK(IU), storing result in U
1333 * (Workspace: need N*N)
1334 *
1335  CALL dgemm( 'N', 'N', m, n, n, one, a, lda,
1336  $ work( iu ), ldwrku, zero, u, ldu )
1337 *
1338  ELSE
1339 *
1340 * Insufficient workspace for a fast algorithm
1341 *
1342  itau = 1
1343  iwork = itau + n
1344 *
1345 * Compute A=Q*R, copying result to U
1346 * (Workspace: need 2*N, prefer N+N*NB)
1347 *
1348  CALL dgeqrf( m, n, a, lda, work( itau ),
1349  $ work( iwork ), lwork-iwork+1, ierr )
1350  CALL dlacpy( 'L', m, n, a, lda, u, ldu )
1351 *
1352 * Generate Q in U
1353 * (Workspace: need 2*N, prefer N+N*NB)
1354 *
1355  CALL dorgqr( m, n, n, u, ldu, work( itau ),
1356  $ work( iwork ), lwork-iwork+1, ierr )
1357 *
1358 * Copy R to VT, zeroing out below it
1359 *
1360  CALL dlacpy( 'U', n, n, a, lda, vt, ldvt )
1361  CALL dlaset( 'L', n-1, n-1, zero, zero, vt( 2, 1 ),
1362  $ ldvt )
1363  ie = itau
1364  itauq = ie + n
1365  itaup = itauq + n
1366  iwork = itaup + n
1367 *
1368 * Bidiagonalize R in VT
1369 * (Workspace: need 4*N, prefer 3*N+2*N*NB)
1370 *
1371  CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
1372  $ work( itauq ), work( itaup ),
1373  $ work( iwork ), lwork-iwork+1, ierr )
1374 *
1375 * Multiply Q in U by left bidiagonalizing vectors
1376 * in VT
1377 * (Workspace: need 3*N+M, prefer 3*N+M*NB)
1378 *
1379  CALL dormbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
1380  $ work( itauq ), u, ldu, work( iwork ),
1381  $ lwork-iwork+1, ierr )
1382 *
1383 * Generate right bidiagonalizing vectors in VT
1384 * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1385 *
1386  CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1387  $ work( iwork ), lwork-iwork+1, ierr )
1388  iwork = ie + n
1389 *
1390 * Perform bidiagonal QR iteration, computing left
1391 * singular vectors of A in U and computing right
1392 * singular vectors of A in VT
1393 * (Workspace: need BDSPAC)
1394 *
1395  CALL dbdsqr( 'U', n, n, m, 0, s, work( ie ), vt,
1396  $ ldvt, u, ldu, dum, 1, work( iwork ),
1397  $ info )
1398 *
1399  END IF
1400 *
1401  END IF
1402 *
1403  ELSE IF( wntua ) THEN
1404 *
1405  IF( wntvn ) THEN
1406 *
1407 * Path 7 (M much larger than N, JOBU='A', JOBVT='N')
1408 * M left singular vectors to be computed in U and
1409 * no right singular vectors to be computed
1410 *
1411  IF( lwork.GE.n*n+max( n+m, 4*n, bdspac ) ) THEN
1412 *
1413 * Sufficient workspace for a fast algorithm
1414 *
1415  ir = 1
1416  IF( lwork.GE.wrkbl+lda*n ) THEN
1417 *
1418 * WORK(IR) is LDA by N
1419 *
1420  ldwrkr = lda
1421  ELSE
1422 *
1423 * WORK(IR) is N by N
1424 *
1425  ldwrkr = n
1426  END IF
1427  itau = ir + ldwrkr*n
1428  iwork = itau + n
1429 *
1430 * Compute A=Q*R, copying result to U
1431 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1432 *
1433  CALL dgeqrf( m, n, a, lda, work( itau ),
1434  $ work( iwork ), lwork-iwork+1, ierr )
1435  CALL dlacpy( 'L', m, n, a, lda, u, ldu )
1436 *
1437 * Copy R to WORK(IR), zeroing out below it
1438 *
1439  CALL dlacpy( 'U', n, n, a, lda, work( ir ),
1440  $ ldwrkr )
1441  CALL dlaset( 'L', n-1, n-1, zero, zero,
1442  $ work( ir+1 ), ldwrkr )
1443 *
1444 * Generate Q in U
1445 * (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
1446 *
1447  CALL dorgqr( m, m, n, u, ldu, work( itau ),
1448  $ work( iwork ), lwork-iwork+1, ierr )
1449  ie = itau
1450  itauq = ie + n
1451  itaup = itauq + n
1452  iwork = itaup + n
1453 *
1454 * Bidiagonalize R in WORK(IR)
1455 * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1456 *
1457  CALL dgebrd( n, n, work( ir ), ldwrkr, s,
1458  $ work( ie ), work( itauq ),
1459  $ work( itaup ), work( iwork ),
1460  $ lwork-iwork+1, ierr )
1461 *
1462 * Generate left bidiagonalizing vectors in WORK(IR)
1463 * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1464 *
1465  CALL dorgbr( 'Q', n, n, n, work( ir ), ldwrkr,
1466  $ work( itauq ), work( iwork ),
1467  $ lwork-iwork+1, ierr )
1468  iwork = ie + n
1469 *
1470 * Perform bidiagonal QR iteration, computing left
1471 * singular vectors of R in WORK(IR)
1472 * (Workspace: need N*N+BDSPAC)
1473 *
1474  CALL dbdsqr( 'U', n, 0, n, 0, s, work( ie ), dum,
1475  $ 1, work( ir ), ldwrkr, dum, 1,
1476  $ work( iwork ), info )
1477 *
1478 * Multiply Q in U by left singular vectors of R in
1479 * WORK(IR), storing result in A
1480 * (Workspace: need N*N)
1481 *
1482  CALL dgemm( 'N', 'N', m, n, n, one, u, ldu,
1483  $ work( ir ), ldwrkr, zero, a, lda )
1484 *
1485 * Copy left singular vectors of A from A to U
1486 *
1487  CALL dlacpy( 'F', m, n, a, lda, u, ldu )
1488 *
1489  ELSE
1490 *
1491 * Insufficient workspace for a fast algorithm
1492 *
1493  itau = 1
1494  iwork = itau + n
1495 *
1496 * Compute A=Q*R, copying result to U
1497 * (Workspace: need 2*N, prefer N+N*NB)
1498 *
1499  CALL dgeqrf( m, n, a, lda, work( itau ),
1500  $ work( iwork ), lwork-iwork+1, ierr )
1501  CALL dlacpy( 'L', m, n, a, lda, u, ldu )
1502 *
1503 * Generate Q in U
1504 * (Workspace: need N+M, prefer N+M*NB)
1505 *
1506  CALL dorgqr( m, m, n, u, ldu, work( itau ),
1507  $ work( iwork ), lwork-iwork+1, ierr )
1508  ie = itau
1509  itauq = ie + n
1510  itaup = itauq + n
1511  iwork = itaup + n
1512 *
1513 * Zero out below R in A
1514 *
1515  CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ),
1516  $ lda )
1517 *
1518 * Bidiagonalize R in A
1519 * (Workspace: need 4*N, prefer 3*N+2*N*NB)
1520 *
1521  CALL dgebrd( n, n, a, lda, s, work( ie ),
1522  $ work( itauq ), work( itaup ),
1523  $ work( iwork ), lwork-iwork+1, ierr )
1524 *
1525 * Multiply Q in U by left bidiagonalizing vectors
1526 * in A
1527 * (Workspace: need 3*N+M, prefer 3*N+M*NB)
1528 *
1529  CALL dormbr( 'Q', 'R', 'N', m, n, n, a, lda,
1530  $ work( itauq ), u, ldu, work( iwork ),
1531  $ lwork-iwork+1, ierr )
1532  iwork = ie + n
1533 *
1534 * Perform bidiagonal QR iteration, computing left
1535 * singular vectors of A in U
1536 * (Workspace: need BDSPAC)
1537 *
1538  CALL dbdsqr( 'U', n, 0, m, 0, s, work( ie ), dum,
1539  $ 1, u, ldu, dum, 1, work( iwork ),
1540  $ info )
1541 *
1542  END IF
1543 *
1544  ELSE IF( wntvo ) THEN
1545 *
1546 * Path 8 (M much larger than N, JOBU='A', JOBVT='O')
1547 * M left singular vectors to be computed in U and
1548 * N right singular vectors to be overwritten on A
1549 *
1550  IF( lwork.GE.2*n*n+max( n+m, 4*n, bdspac ) ) THEN
1551 *
1552 * Sufficient workspace for a fast algorithm
1553 *
1554  iu = 1
1555  IF( lwork.GE.wrkbl+2*lda*n ) THEN
1556 *
1557 * WORK(IU) is LDA by N and WORK(IR) is LDA by N
1558 *
1559  ldwrku = lda
1560  ir = iu + ldwrku*n
1561  ldwrkr = lda
1562  ELSE IF( lwork.GE.wrkbl+( lda+n )*n ) THEN
1563 *
1564 * WORK(IU) is LDA by N and WORK(IR) is N by N
1565 *
1566  ldwrku = lda
1567  ir = iu + ldwrku*n
1568  ldwrkr = n
1569  ELSE
1570 *
1571 * WORK(IU) is N by N and WORK(IR) is N by N
1572 *
1573  ldwrku = n
1574  ir = iu + ldwrku*n
1575  ldwrkr = n
1576  END IF
1577  itau = ir + ldwrkr*n
1578  iwork = itau + n
1579 *
1580 * Compute A=Q*R, copying result to U
1581 * (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1582 *
1583  CALL dgeqrf( m, n, a, lda, work( itau ),
1584  $ work( iwork ), lwork-iwork+1, ierr )
1585  CALL dlacpy( 'L', m, n, a, lda, u, ldu )
1586 *
1587 * Generate Q in U
1588 * (Workspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
1589 *
1590  CALL dorgqr( m, m, n, u, ldu, work( itau ),
1591  $ work( iwork ), lwork-iwork+1, ierr )
1592 *
1593 * Copy R to WORK(IU), zeroing out below it
1594 *
1595  CALL dlacpy( 'U', n, n, a, lda, work( iu ),
1596  $ ldwrku )
1597  CALL dlaset( 'L', n-1, n-1, zero, zero,
1598  $ work( iu+1 ), ldwrku )
1599  ie = itau
1600  itauq = ie + n
1601  itaup = itauq + n
1602  iwork = itaup + n
1603 *
1604 * Bidiagonalize R in WORK(IU), copying result to
1605 * WORK(IR)
1606 * (Workspace: need 2*N*N+4*N,
1607 * prefer 2*N*N+3*N+2*N*NB)
1608 *
1609  CALL dgebrd( n, n, work( iu ), ldwrku, s,
1610  $ work( ie ), work( itauq ),
1611  $ work( itaup ), work( iwork ),
1612  $ lwork-iwork+1, ierr )
1613  CALL dlacpy( 'U', n, n, work( iu ), ldwrku,
1614  $ work( ir ), ldwrkr )
1615 *
1616 * Generate left bidiagonalizing vectors in WORK(IU)
1617 * (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
1618 *
1619  CALL dorgbr( 'Q', n, n, n, work( iu ), ldwrku,
1620  $ work( itauq ), work( iwork ),
1621  $ lwork-iwork+1, ierr )
1622 *
1623 * Generate right bidiagonalizing vectors in WORK(IR)
1624 * (Workspace: need 2*N*N+4*N-1,
1625 * prefer 2*N*N+3*N+(N-1)*NB)
1626 *
1627  CALL dorgbr( 'P', n, n, n, work( ir ), ldwrkr,
1628  $ work( itaup ), work( iwork ),
1629  $ lwork-iwork+1, ierr )
1630  iwork = ie + n
1631 *
1632 * Perform bidiagonal QR iteration, computing left
1633 * singular vectors of R in WORK(IU) and computing
1634 * right singular vectors of R in WORK(IR)
1635 * (Workspace: need 2*N*N+BDSPAC)
1636 *
1637  CALL dbdsqr( 'U', n, n, n, 0, s, work( ie ),
1638  $ work( ir ), ldwrkr, work( iu ),
1639  $ ldwrku, dum, 1, work( iwork ), info )
1640 *
1641 * Multiply Q in U by left singular vectors of R in
1642 * WORK(IU), storing result in A
1643 * (Workspace: need N*N)
1644 *
1645  CALL dgemm( 'N', 'N', m, n, n, one, u, ldu,
1646  $ work( iu ), ldwrku, zero, a, lda )
1647 *
1648 * Copy left singular vectors of A from A to U
1649 *
1650  CALL dlacpy( 'F', m, n, a, lda, u, ldu )
1651 *
1652 * Copy right singular vectors of R from WORK(IR) to A
1653 *
1654  CALL dlacpy( 'F', n, n, work( ir ), ldwrkr, a,
1655  $ lda )
1656 *
1657  ELSE
1658 *
1659 * Insufficient workspace for a fast algorithm
1660 *
1661  itau = 1
1662  iwork = itau + n
1663 *
1664 * Compute A=Q*R, copying result to U
1665 * (Workspace: need 2*N, prefer N+N*NB)
1666 *
1667  CALL dgeqrf( m, n, a, lda, work( itau ),
1668  $ work( iwork ), lwork-iwork+1, ierr )
1669  CALL dlacpy( 'L', m, n, a, lda, u, ldu )
1670 *
1671 * Generate Q in U
1672 * (Workspace: need N+M, prefer N+M*NB)
1673 *
1674  CALL dorgqr( m, m, n, u, ldu, work( itau ),
1675  $ work( iwork ), lwork-iwork+1, ierr )
1676  ie = itau
1677  itauq = ie + n
1678  itaup = itauq + n
1679  iwork = itaup + n
1680 *
1681 * Zero out below R in A
1682 *
1683  CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ),
1684  $ lda )
1685 *
1686 * Bidiagonalize R in A
1687 * (Workspace: need 4*N, prefer 3*N+2*N*NB)
1688 *
1689  CALL dgebrd( n, n, a, lda, s, work( ie ),
1690  $ work( itauq ), work( itaup ),
1691  $ work( iwork ), lwork-iwork+1, ierr )
1692 *
1693 * Multiply Q in U by left bidiagonalizing vectors
1694 * in A
1695 * (Workspace: need 3*N+M, prefer 3*N+M*NB)
1696 *
1697  CALL dormbr( 'Q', 'R', 'N', m, n, n, a, lda,
1698  $ work( itauq ), u, ldu, work( iwork ),
1699  $ lwork-iwork+1, ierr )
1700 *
1701 * Generate right bidiagonalizing vectors in A
1702 * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1703 *
1704  CALL dorgbr( 'P', n, n, n, a, lda, work( itaup ),
1705  $ work( iwork ), lwork-iwork+1, ierr )
1706  iwork = ie + n
1707 *
1708 * Perform bidiagonal QR iteration, computing left
1709 * singular vectors of A in U and computing right
1710 * singular vectors of A in A
1711 * (Workspace: need BDSPAC)
1712 *
1713  CALL dbdsqr( 'U', n, n, m, 0, s, work( ie ), a,
1714  $ lda, u, ldu, dum, 1, work( iwork ),
1715  $ info )
1716 *
1717  END IF
1718 *
1719  ELSE IF( wntvas ) THEN
1720 *
1721 * Path 9 (M much larger than N, JOBU='A', JOBVT='S'
1722 * or 'A')
1723 * M left singular vectors to be computed in U and
1724 * N right singular vectors to be computed in VT
1725 *
1726  IF( lwork.GE.n*n+max( n+m, 4*n, bdspac ) ) THEN
1727 *
1728 * Sufficient workspace for a fast algorithm
1729 *
1730  iu = 1
1731  IF( lwork.GE.wrkbl+lda*n ) THEN
1732 *
1733 * WORK(IU) is LDA by N
1734 *
1735  ldwrku = lda
1736  ELSE
1737 *
1738 * WORK(IU) is N by N
1739 *
1740  ldwrku = n
1741  END IF
1742  itau = iu + ldwrku*n
1743  iwork = itau + n
1744 *
1745 * Compute A=Q*R, copying result to U
1746 * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1747 *
1748  CALL dgeqrf( m, n, a, lda, work( itau ),
1749  $ work( iwork ), lwork-iwork+1, ierr )
1750  CALL dlacpy( 'L', m, n, a, lda, u, ldu )
1751 *
1752 * Generate Q in U
1753 * (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
1754 *
1755  CALL dorgqr( m, m, n, u, ldu, work( itau ),
1756  $ work( iwork ), lwork-iwork+1, ierr )
1757 *
1758 * Copy R to WORK(IU), zeroing out below it
1759 *
1760  CALL dlacpy( 'U', n, n, a, lda, work( iu ),
1761  $ ldwrku )
1762  CALL dlaset( 'L', n-1, n-1, zero, zero,
1763  $ work( iu+1 ), ldwrku )
1764  ie = itau
1765  itauq = ie + n
1766  itaup = itauq + n
1767  iwork = itaup + n
1768 *
1769 * Bidiagonalize R in WORK(IU), copying result to VT
1770 * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1771 *
1772  CALL dgebrd( n, n, work( iu ), ldwrku, s,
1773  $ work( ie ), work( itauq ),
1774  $ work( itaup ), work( iwork ),
1775  $ lwork-iwork+1, ierr )
1776  CALL dlacpy( 'U', n, n, work( iu ), ldwrku, vt,
1777  $ ldvt )
1778 *
1779 * Generate left bidiagonalizing vectors in WORK(IU)
1780 * (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1781 *
1782  CALL dorgbr( 'Q', n, n, n, work( iu ), ldwrku,
1783  $ work( itauq ), work( iwork ),
1784  $ lwork-iwork+1, ierr )
1785 *
1786 * Generate right bidiagonalizing vectors in VT
1787 * (Workspace: need N*N+4*N-1,
1788 * prefer N*N+3*N+(N-1)*NB)
1789 *
1790  CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1791  $ work( iwork ), lwork-iwork+1, ierr )
1792  iwork = ie + n
1793 *
1794 * Perform bidiagonal QR iteration, computing left
1795 * singular vectors of R in WORK(IU) and computing
1796 * right singular vectors of R in VT
1797 * (Workspace: need N*N+BDSPAC)
1798 *
1799  CALL dbdsqr( 'U', n, n, n, 0, s, work( ie ), vt,
1800  $ ldvt, work( iu ), ldwrku, dum, 1,
1801  $ work( iwork ), info )
1802 *
1803 * Multiply Q in U by left singular vectors of R in
1804 * WORK(IU), storing result in A
1805 * (Workspace: need N*N)
1806 *
1807  CALL dgemm( 'N', 'N', m, n, n, one, u, ldu,
1808  $ work( iu ), ldwrku, zero, a, lda )
1809 *
1810 * Copy left singular vectors of A from A to U
1811 *
1812  CALL dlacpy( 'F', m, n, a, lda, u, ldu )
1813 *
1814  ELSE
1815 *
1816 * Insufficient workspace for a fast algorithm
1817 *
1818  itau = 1
1819  iwork = itau + n
1820 *
1821 * Compute A=Q*R, copying result to U
1822 * (Workspace: need 2*N, prefer N+N*NB)
1823 *
1824  CALL dgeqrf( m, n, a, lda, work( itau ),
1825  $ work( iwork ), lwork-iwork+1, ierr )
1826  CALL dlacpy( 'L', m, n, a, lda, u, ldu )
1827 *
1828 * Generate Q in U
1829 * (Workspace: need N+M, prefer N+M*NB)
1830 *
1831  CALL dorgqr( m, m, n, u, ldu, work( itau ),
1832  $ work( iwork ), lwork-iwork+1, ierr )
1833 *
1834 * Copy R from A to VT, zeroing out below it
1835 *
1836  CALL dlacpy( 'U', n, n, a, lda, vt, ldvt )
1837  CALL dlaset( 'L', n-1, n-1, zero, zero, vt( 2, 1 ),
1838  $ ldvt )
1839  ie = itau
1840  itauq = ie + n
1841  itaup = itauq + n
1842  iwork = itaup + n
1843 *
1844 * Bidiagonalize R in VT
1845 * (Workspace: need 4*N, prefer 3*N+2*N*NB)
1846 *
1847  CALL dgebrd( n, n, vt, ldvt, s, work( ie ),
1848  $ work( itauq ), work( itaup ),
1849  $ work( iwork ), lwork-iwork+1, ierr )
1850 *
1851 * Multiply Q in U by left bidiagonalizing vectors
1852 * in VT
1853 * (Workspace: need 3*N+M, prefer 3*N+M*NB)
1854 *
1855  CALL dormbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
1856  $ work( itauq ), u, ldu, work( iwork ),
1857  $ lwork-iwork+1, ierr )
1858 *
1859 * Generate right bidiagonalizing vectors in VT
1860 * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1861 *
1862  CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1863  $ work( iwork ), lwork-iwork+1, ierr )
1864  iwork = ie + n
1865 *
1866 * Perform bidiagonal QR iteration, computing left
1867 * singular vectors of A in U and computing right
1868 * singular vectors of A in VT
1869 * (Workspace: need BDSPAC)
1870 *
1871  CALL dbdsqr( 'U', n, n, m, 0, s, work( ie ), vt,
1872  $ ldvt, u, ldu, dum, 1, work( iwork ),
1873  $ info )
1874 *
1875  END IF
1876 *
1877  END IF
1878 *
1879  END IF
1880 *
1881  ELSE
1882 *
1883 * M .LT. MNTHR
1884 *
1885 * Path 10 (M at least N, but not much larger)
1886 * Reduce to bidiagonal form without QR decomposition
1887 *
1888  ie = 1
1889  itauq = ie + n
1890  itaup = itauq + n
1891  iwork = itaup + n
1892 *
1893 * Bidiagonalize A
1894 * (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
1895 *
1896  CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
1897  $ work( itaup ), work( iwork ), lwork-iwork+1,
1898  $ ierr )
1899  IF( wntuas ) THEN
1900 *
1901 * If left singular vectors desired in U, copy result to U
1902 * and generate left bidiagonalizing vectors in U
1903 * (Workspace: need 3*N+NCU, prefer 3*N+NCU*NB)
1904 *
1905  CALL dlacpy( 'L', m, n, a, lda, u, ldu )
1906  IF( wntus )
1907  $ ncu = n
1908  IF( wntua )
1909  $ ncu = m
1910  CALL dorgbr( 'Q', m, ncu, n, u, ldu, work( itauq ),
1911  $ work( iwork ), lwork-iwork+1, ierr )
1912  END IF
1913  IF( wntvas ) THEN
1914 *
1915 * If right singular vectors desired in VT, copy result to
1916 * VT and generate right bidiagonalizing vectors in VT
1917 * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1918 *
1919  CALL dlacpy( 'U', n, n, a, lda, vt, ldvt )
1920  CALL dorgbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1921  $ work( iwork ), lwork-iwork+1, ierr )
1922  END IF
1923  IF( wntuo ) THEN
1924 *
1925 * If left singular vectors desired in A, generate left
1926 * bidiagonalizing vectors in A
1927 * (Workspace: need 4*N, prefer 3*N+N*NB)
1928 *
1929  CALL dorgbr( 'Q', m, n, n, a, lda, work( itauq ),
1930  $ work( iwork ), lwork-iwork+1, ierr )
1931  END IF
1932  IF( wntvo ) THEN
1933 *
1934 * If right singular vectors desired in A, generate right
1935 * bidiagonalizing vectors in A
1936 * (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1937 *
1938  CALL dorgbr( 'P', n, n, n, a, lda, work( itaup ),
1939  $ work( iwork ), lwork-iwork+1, ierr )
1940  END IF
1941  iwork = ie + n
1942  IF( wntuas .OR. wntuo )
1943  $ nru = m
1944  IF( wntun )
1945  $ nru = 0
1946  IF( wntvas .OR. wntvo )
1947  $ ncvt = n
1948  IF( wntvn )
1949  $ ncvt = 0
1950  IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) ) THEN
1951 *
1952 * Perform bidiagonal QR iteration, if desired, computing
1953 * left singular vectors in U and computing right singular
1954 * vectors in VT
1955 * (Workspace: need BDSPAC)
1956 *
1957  CALL dbdsqr( 'U', n, ncvt, nru, 0, s, work( ie ), vt,
1958  $ ldvt, u, ldu, dum, 1, work( iwork ), info )
1959  ELSE IF( ( .NOT.wntuo ) .AND. wntvo ) THEN
1960 *
1961 * Perform bidiagonal QR iteration, if desired, computing
1962 * left singular vectors in U and computing right singular
1963 * vectors in A
1964 * (Workspace: need BDSPAC)
1965 *
1966  CALL dbdsqr( 'U', n, ncvt, nru, 0, s, work( ie ), a, lda,
1967  $ u, ldu, dum, 1, work( iwork ), info )
1968  ELSE
1969 *
1970 * Perform bidiagonal QR iteration, if desired, computing
1971 * left singular vectors in A and computing right singular
1972 * vectors in VT
1973 * (Workspace: need BDSPAC)
1974 *
1975  CALL dbdsqr( 'U', n, ncvt, nru, 0, s, work( ie ), vt,
1976  $ ldvt, a, lda, dum, 1, work( iwork ), info )
1977  END IF
1978 *
1979  END IF
1980 *
1981  ELSE
1982 *
1983 * A has more columns than rows. If A has sufficiently more
1984 * columns than rows, first reduce using the LQ decomposition (if
1985 * sufficient workspace available)
1986 *
1987  IF( n.GE.mnthr ) THEN
1988 *
1989  IF( wntvn ) THEN
1990 *
1991 * Path 1t(N much larger than M, JOBVT='N')
1992 * No right singular vectors to be computed
1993 *
1994  itau = 1
1995  iwork = itau + m
1996 *
1997 * Compute A=L*Q
1998 * (Workspace: need 2*M, prefer M+M*NB)
1999 *
2000  CALL dgelqf( m, n, a, lda, work( itau ), work( iwork ),
2001  $ lwork-iwork+1, ierr )
2002 *
2003 * Zero out above L
2004 *
2005  CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
2006  ie = 1
2007  itauq = ie + m
2008  itaup = itauq + m
2009  iwork = itaup + m
2010 *
2011 * Bidiagonalize L in A
2012 * (Workspace: need 4*M, prefer 3*M+2*M*NB)
2013 *
2014  CALL dgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
2015  $ work( itaup ), work( iwork ), lwork-iwork+1,
2016  $ ierr )
2017  IF( wntuo .OR. wntuas ) THEN
2018 *
2019 * If left singular vectors desired, generate Q
2020 * (Workspace: need 4*M, prefer 3*M+M*NB)
2021 *
2022  CALL dorgbr( 'Q', m, m, m, a, lda, work( itauq ),
2023  $ work( iwork ), lwork-iwork+1, ierr )
2024  END IF
2025  iwork = ie + m
2026  nru = 0
2027  IF( wntuo .OR. wntuas )
2028  $ nru = m
2029 *
2030 * Perform bidiagonal QR iteration, computing left singular
2031 * vectors of A in A if desired
2032 * (Workspace: need BDSPAC)
2033 *
2034  CALL dbdsqr( 'U', m, 0, nru, 0, s, work( ie ), dum, 1, a,
2035  $ lda, dum, 1, work( iwork ), info )
2036 *
2037 * If left singular vectors desired in U, copy them there
2038 *
2039  IF( wntuas )
2040  $ CALL dlacpy( 'F', m, m, a, lda, u, ldu )
2041 *
2042  ELSE IF( wntvo .AND. wntun ) THEN
2043 *
2044 * Path 2t(N much larger than M, JOBU='N', JOBVT='O')
2045 * M right singular vectors to be overwritten on A and
2046 * no left singular vectors to be computed
2047 *
2048  IF( lwork.GE.m*m+max( 4*m, bdspac ) ) THEN
2049 *
2050 * Sufficient workspace for a fast algorithm
2051 *
2052  ir = 1
2053  IF( lwork.GE.max( wrkbl, lda*n+m )+lda*m ) THEN
2054 *
2055 * WORK(IU) is LDA by N and WORK(IR) is LDA by M
2056 *
2057  ldwrku = lda
2058  chunk = n
2059  ldwrkr = lda
2060  ELSE IF( lwork.GE.max( wrkbl, lda*n+m )+m*m ) THEN
2061 *
2062 * WORK(IU) is LDA by N and WORK(IR) is M by M
2063 *
2064  ldwrku = lda
2065  chunk = n
2066  ldwrkr = m
2067  ELSE
2068 *
2069 * WORK(IU) is M by CHUNK and WORK(IR) is M by M
2070 *
2071  ldwrku = m
2072  chunk = ( lwork-m*m-m ) / m
2073  ldwrkr = m
2074  END IF
2075  itau = ir + ldwrkr*m
2076  iwork = itau + m
2077 *
2078 * Compute A=L*Q
2079 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2080 *
2081  CALL dgelqf( m, n, a, lda, work( itau ),
2082  $ work( iwork ), lwork-iwork+1, ierr )
2083 *
2084 * Copy L to WORK(IR) and zero out above it
2085 *
2086  CALL dlacpy( 'L', m, m, a, lda, work( ir ), ldwrkr )
2087  CALL dlaset( 'U', m-1, m-1, zero, zero,
2088  $ work( ir+ldwrkr ), ldwrkr )
2089 *
2090 * Generate Q in A
2091 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2092 *
2093  CALL dorglq( m, n, m, a, lda, work( itau ),
2094  $ work( iwork ), lwork-iwork+1, ierr )
2095  ie = itau
2096  itauq = ie + m
2097  itaup = itauq + m
2098  iwork = itaup + m
2099 *
2100 * Bidiagonalize L in WORK(IR)
2101 * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
2102 *
2103  CALL dgebrd( m, m, work( ir ), ldwrkr, s, work( ie ),
2104  $ work( itauq ), work( itaup ),
2105  $ work( iwork ), lwork-iwork+1, ierr )
2106 *
2107 * Generate right vectors bidiagonalizing L
2108 * (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
2109 *
2110  CALL dorgbr( 'P', m, m, m, work( ir ), ldwrkr,
2111  $ work( itaup ), work( iwork ),
2112  $ lwork-iwork+1, ierr )
2113  iwork = ie + m
2114 *
2115 * Perform bidiagonal QR iteration, computing right
2116 * singular vectors of L in WORK(IR)
2117 * (Workspace: need M*M+BDSPAC)
2118 *
2119  CALL dbdsqr( 'U', m, m, 0, 0, s, work( ie ),
2120  $ work( ir ), ldwrkr, dum, 1, dum, 1,
2121  $ work( iwork ), info )
2122  iu = ie + m
2123 *
2124 * Multiply right singular vectors of L in WORK(IR) by Q
2125 * in A, storing result in WORK(IU) and copying to A
2126 * (Workspace: need M*M+2*M, prefer M*M+M*N+M)
2127 *
2128  DO 30 i = 1, n, chunk
2129  blk = min( n-i+1, chunk )
2130  CALL dgemm( 'N', 'N', m, blk, m, one, work( ir ),
2131  $ ldwrkr, a( 1, i ), lda, zero,
2132  $ work( iu ), ldwrku )
2133  CALL dlacpy( 'F', m, blk, work( iu ), ldwrku,
2134  $ a( 1, i ), lda )
2135  30 CONTINUE
2136 *
2137  ELSE
2138 *
2139 * Insufficient workspace for a fast algorithm
2140 *
2141  ie = 1
2142  itauq = ie + m
2143  itaup = itauq + m
2144  iwork = itaup + m
2145 *
2146 * Bidiagonalize A
2147 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
2148 *
2149  CALL dgebrd( m, n, a, lda, s, work( ie ),
2150  $ work( itauq ), work( itaup ),
2151  $ work( iwork ), lwork-iwork+1, ierr )
2152 *
2153 * Generate right vectors bidiagonalizing A
2154 * (Workspace: need 4*M, prefer 3*M+M*NB)
2155 *
2156  CALL dorgbr( 'P', m, n, m, a, lda, work( itaup ),
2157  $ work( iwork ), lwork-iwork+1, ierr )
2158  iwork = ie + m
2159 *
2160 * Perform bidiagonal QR iteration, computing right
2161 * singular vectors of A in A
2162 * (Workspace: need BDSPAC)
2163 *
2164  CALL dbdsqr( 'L', m, n, 0, 0, s, work( ie ), a, lda,
2165  $ dum, 1, dum, 1, work( iwork ), info )
2166 *
2167  END IF
2168 *
2169  ELSE IF( wntvo .AND. wntuas ) THEN
2170 *
2171 * Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
2172 * M right singular vectors to be overwritten on A and
2173 * M left singular vectors to be computed in U
2174 *
2175  IF( lwork.GE.m*m+max( 4*m, bdspac ) ) THEN
2176 *
2177 * Sufficient workspace for a fast algorithm
2178 *
2179  ir = 1
2180  IF( lwork.GE.max( wrkbl, lda*n+m )+lda*m ) THEN
2181 *
2182 * WORK(IU) is LDA by N and WORK(IR) is LDA by M
2183 *
2184  ldwrku = lda
2185  chunk = n
2186  ldwrkr = lda
2187  ELSE IF( lwork.GE.max( wrkbl, lda*n+m )+m*m ) THEN
2188 *
2189 * WORK(IU) is LDA by N and WORK(IR) is M by M
2190 *
2191  ldwrku = lda
2192  chunk = n
2193  ldwrkr = m
2194  ELSE
2195 *
2196 * WORK(IU) is M by CHUNK and WORK(IR) is M by M
2197 *
2198  ldwrku = m
2199  chunk = ( lwork-m*m-m ) / m
2200  ldwrkr = m
2201  END IF
2202  itau = ir + ldwrkr*m
2203  iwork = itau + m
2204 *
2205 * Compute A=L*Q
2206 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2207 *
2208  CALL dgelqf( m, n, a, lda, work( itau ),
2209  $ work( iwork ), lwork-iwork+1, ierr )
2210 *
2211 * Copy L to U, zeroing about above it
2212 *
2213  CALL dlacpy( 'L', m, m, a, lda, u, ldu )
2214  CALL dlaset( 'U', m-1, m-1, zero, zero, u( 1, 2 ),
2215  $ ldu )
2216 *
2217 * Generate Q in A
2218 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2219 *
2220  CALL dorglq( m, n, m, a, lda, work( itau ),
2221  $ work( iwork ), lwork-iwork+1, ierr )
2222  ie = itau
2223  itauq = ie + m
2224  itaup = itauq + m
2225  iwork = itaup + m
2226 *
2227 * Bidiagonalize L in U, copying result to WORK(IR)
2228 * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
2229 *
2230  CALL dgebrd( m, m, u, ldu, s, work( ie ),
2231  $ work( itauq ), work( itaup ),
2232  $ work( iwork ), lwork-iwork+1, ierr )
2233  CALL dlacpy( 'U', m, m, u, ldu, work( ir ), ldwrkr )
2234 *
2235 * Generate right vectors bidiagonalizing L in WORK(IR)
2236 * (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
2237 *
2238  CALL dorgbr( 'P', m, m, m, work( ir ), ldwrkr,
2239  $ work( itaup ), work( iwork ),
2240  $ lwork-iwork+1, ierr )
2241 *
2242 * Generate left vectors bidiagonalizing L in U
2243 * (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
2244 *
2245  CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
2246  $ work( iwork ), lwork-iwork+1, ierr )
2247  iwork = ie + m
2248 *
2249 * Perform bidiagonal QR iteration, computing left
2250 * singular vectors of L in U, and computing right
2251 * singular vectors of L in WORK(IR)
2252 * (Workspace: need M*M+BDSPAC)
2253 *
2254  CALL dbdsqr( 'U', m, m, m, 0, s, work( ie ),
2255  $ work( ir ), ldwrkr, u, ldu, dum, 1,
2256  $ work( iwork ), info )
2257  iu = ie + m
2258 *
2259 * Multiply right singular vectors of L in WORK(IR) by Q
2260 * in A, storing result in WORK(IU) and copying to A
2261 * (Workspace: need M*M+2*M, prefer M*M+M*N+M))
2262 *
2263  DO 40 i = 1, n, chunk
2264  blk = min( n-i+1, chunk )
2265  CALL dgemm( 'N', 'N', m, blk, m, one, work( ir ),
2266  $ ldwrkr, a( 1, i ), lda, zero,
2267  $ work( iu ), ldwrku )
2268  CALL dlacpy( 'F', m, blk, work( iu ), ldwrku,
2269  $ a( 1, i ), lda )
2270  40 CONTINUE
2271 *
2272  ELSE
2273 *
2274 * Insufficient workspace for a fast algorithm
2275 *
2276  itau = 1
2277  iwork = itau + m
2278 *
2279 * Compute A=L*Q
2280 * (Workspace: need 2*M, prefer M+M*NB)
2281 *
2282  CALL dgelqf( m, n, a, lda, work( itau ),
2283  $ work( iwork ), lwork-iwork+1, ierr )
2284 *
2285 * Copy L to U, zeroing out above it
2286 *
2287  CALL dlacpy( 'L', m, m, a, lda, u, ldu )
2288  CALL dlaset( 'U', m-1, m-1, zero, zero, u( 1, 2 ),
2289  $ ldu )
2290 *
2291 * Generate Q in A
2292 * (Workspace: need 2*M, prefer M+M*NB)
2293 *
2294  CALL dorglq( m, n, m, a, lda, work( itau ),
2295  $ work( iwork ), lwork-iwork+1, ierr )
2296  ie = itau
2297  itauq = ie + m
2298  itaup = itauq + m
2299  iwork = itaup + m
2300 *
2301 * Bidiagonalize L in U
2302 * (Workspace: need 4*M, prefer 3*M+2*M*NB)
2303 *
2304  CALL dgebrd( m, m, u, ldu, s, work( ie ),
2305  $ work( itauq ), work( itaup ),
2306  $ work( iwork ), lwork-iwork+1, ierr )
2307 *
2308 * Multiply right vectors bidiagonalizing L by Q in A
2309 * (Workspace: need 3*M+N, prefer 3*M+N*NB)
2310 *
2311  CALL dormbr( 'P', 'L', 'T', m, n, m, u, ldu,
2312  $ work( itaup ), a, lda, work( iwork ),
2313  $ lwork-iwork+1, ierr )
2314 *
2315 * Generate left vectors bidiagonalizing L in U
2316 * (Workspace: need 4*M, prefer 3*M+M*NB)
2317 *
2318  CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
2319  $ work( iwork ), lwork-iwork+1, ierr )
2320  iwork = ie + m
2321 *
2322 * Perform bidiagonal QR iteration, computing left
2323 * singular vectors of A in U and computing right
2324 * singular vectors of A in A
2325 * (Workspace: need BDSPAC)
2326 *
2327  CALL dbdsqr( 'U', m, n, m, 0, s, work( ie ), a, lda,
2328  $ u, ldu, dum, 1, work( iwork ), info )
2329 *
2330  END IF
2331 *
2332  ELSE IF( wntvs ) THEN
2333 *
2334  IF( wntun ) THEN
2335 *
2336 * Path 4t(N much larger than M, JOBU='N', JOBVT='S')
2337 * M right singular vectors to be computed in VT and
2338 * no left singular vectors to be computed
2339 *
2340  IF( lwork.GE.m*m+max( 4*m, bdspac ) ) THEN
2341 *
2342 * Sufficient workspace for a fast algorithm
2343 *
2344  ir = 1
2345  IF( lwork.GE.wrkbl+lda*m ) THEN
2346 *
2347 * WORK(IR) is LDA by M
2348 *
2349  ldwrkr = lda
2350  ELSE
2351 *
2352 * WORK(IR) is M by M
2353 *
2354  ldwrkr = m
2355  END IF
2356  itau = ir + ldwrkr*m
2357  iwork = itau + m
2358 *
2359 * Compute A=L*Q
2360 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2361 *
2362  CALL dgelqf( m, n, a, lda, work( itau ),
2363  $ work( iwork ), lwork-iwork+1, ierr )
2364 *
2365 * Copy L to WORK(IR), zeroing out above it
2366 *
2367  CALL dlacpy( 'L', m, m, a, lda, work( ir ),
2368  $ ldwrkr )
2369  CALL dlaset( 'U', m-1, m-1, zero, zero,
2370  $ work( ir+ldwrkr ), ldwrkr )
2371 *
2372 * Generate Q in A
2373 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2374 *
2375  CALL dorglq( m, n, m, a, lda, work( itau ),
2376  $ work( iwork ), lwork-iwork+1, ierr )
2377  ie = itau
2378  itauq = ie + m
2379  itaup = itauq + m
2380  iwork = itaup + m
2381 *
2382 * Bidiagonalize L in WORK(IR)
2383 * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
2384 *
2385  CALL dgebrd( m, m, work( ir ), ldwrkr, s,
2386  $ work( ie ), work( itauq ),
2387  $ work( itaup ), work( iwork ),
2388  $ lwork-iwork+1, ierr )
2389 *
2390 * Generate right vectors bidiagonalizing L in
2391 * WORK(IR)
2392 * (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
2393 *
2394  CALL dorgbr( 'P', m, m, m, work( ir ), ldwrkr,
2395  $ work( itaup ), work( iwork ),
2396  $ lwork-iwork+1, ierr )
2397  iwork = ie + m
2398 *
2399 * Perform bidiagonal QR iteration, computing right
2400 * singular vectors of L in WORK(IR)
2401 * (Workspace: need M*M+BDSPAC)
2402 *
2403  CALL dbdsqr( 'U', m, m, 0, 0, s, work( ie ),
2404  $ work( ir ), ldwrkr, dum, 1, dum, 1,
2405  $ work( iwork ), info )
2406 *
2407 * Multiply right singular vectors of L in WORK(IR) by
2408 * Q in A, storing result in VT
2409 * (Workspace: need M*M)
2410 *
2411  CALL dgemm( 'N', 'N', m, n, m, one, work( ir ),
2412  $ ldwrkr, a, lda, zero, vt, ldvt )
2413 *
2414  ELSE
2415 *
2416 * Insufficient workspace for a fast algorithm
2417 *
2418  itau = 1
2419  iwork = itau + m
2420 *
2421 * Compute A=L*Q
2422 * (Workspace: need 2*M, prefer M+M*NB)
2423 *
2424  CALL dgelqf( m, n, a, lda, work( itau ),
2425  $ work( iwork ), lwork-iwork+1, ierr )
2426 *
2427 * Copy result to VT
2428 *
2429  CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
2430 *
2431 * Generate Q in VT
2432 * (Workspace: need 2*M, prefer M+M*NB)
2433 *
2434  CALL dorglq( m, n, m, vt, ldvt, work( itau ),
2435  $ work( iwork ), lwork-iwork+1, ierr )
2436  ie = itau
2437  itauq = ie + m
2438  itaup = itauq + m
2439  iwork = itaup + m
2440 *
2441 * Zero out above L in A
2442 *
2443  CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ),
2444  $ lda )
2445 *
2446 * Bidiagonalize L in A
2447 * (Workspace: need 4*M, prefer 3*M+2*M*NB)
2448 *
2449  CALL dgebrd( m, m, a, lda, s, work( ie ),
2450  $ work( itauq ), work( itaup ),
2451  $ work( iwork ), lwork-iwork+1, ierr )
2452 *
2453 * Multiply right vectors bidiagonalizing L by Q in VT
2454 * (Workspace: need 3*M+N, prefer 3*M+N*NB)
2455 *
2456  CALL dormbr( 'P', 'L', 'T', m, n, m, a, lda,
2457  $ work( itaup ), vt, ldvt,
2458  $ work( iwork ), lwork-iwork+1, ierr )
2459  iwork = ie + m
2460 *
2461 * Perform bidiagonal QR iteration, computing right
2462 * singular vectors of A in VT
2463 * (Workspace: need BDSPAC)
2464 *
2465  CALL dbdsqr( 'U', m, n, 0, 0, s, work( ie ), vt,
2466  $ ldvt, dum, 1, dum, 1, work( iwork ),
2467  $ info )
2468 *
2469  END IF
2470 *
2471  ELSE IF( wntuo ) THEN
2472 *
2473 * Path 5t(N much larger than M, JOBU='O', JOBVT='S')
2474 * M right singular vectors to be computed in VT and
2475 * M left singular vectors to be overwritten on A
2476 *
2477  IF( lwork.GE.2*m*m+max( 4*m, bdspac ) ) THEN
2478 *
2479 * Sufficient workspace for a fast algorithm
2480 *
2481  iu = 1
2482  IF( lwork.GE.wrkbl+2*lda*m ) THEN
2483 *
2484 * WORK(IU) is LDA by M and WORK(IR) is LDA by M
2485 *
2486  ldwrku = lda
2487  ir = iu + ldwrku*m
2488  ldwrkr = lda
2489  ELSE IF( lwork.GE.wrkbl+( lda+m )*m ) THEN
2490 *
2491 * WORK(IU) is LDA by M and WORK(IR) is M by M
2492 *
2493  ldwrku = lda
2494  ir = iu + ldwrku*m
2495  ldwrkr = m
2496  ELSE
2497 *
2498 * WORK(IU) is M by M and WORK(IR) is M by M
2499 *
2500  ldwrku = m
2501  ir = iu + ldwrku*m
2502  ldwrkr = m
2503  END IF
2504  itau = ir + ldwrkr*m
2505  iwork = itau + m
2506 *
2507 * Compute A=L*Q
2508 * (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
2509 *
2510  CALL dgelqf( m, n, a, lda, work( itau ),
2511  $ work( iwork ), lwork-iwork+1, ierr )
2512 *
2513 * Copy L to WORK(IU), zeroing out below it
2514 *
2515  CALL dlacpy( 'L', m, m, a, lda, work( iu ),
2516  $ ldwrku )
2517  CALL dlaset( 'U', m-1, m-1, zero, zero,
2518  $ work( iu+ldwrku ), ldwrku )
2519 *
2520 * Generate Q in A
2521 * (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
2522 *
2523  CALL dorglq( m, n, m, a, lda, work( itau ),
2524  $ work( iwork ), lwork-iwork+1, ierr )
2525  ie = itau
2526  itauq = ie + m
2527  itaup = itauq + m
2528  iwork = itaup + m
2529 *
2530 * Bidiagonalize L in WORK(IU), copying result to
2531 * WORK(IR)
2532 * (Workspace: need 2*M*M+4*M,
2533 * prefer 2*M*M+3*M+2*M*NB)
2534 *
2535  CALL dgebrd( m, m, work( iu ), ldwrku, s,
2536  $ work( ie ), work( itauq ),
2537  $ work( itaup ), work( iwork ),
2538  $ lwork-iwork+1, ierr )
2539  CALL dlacpy( 'L', m, m, work( iu ), ldwrku,
2540  $ work( ir ), ldwrkr )
2541 *
2542 * Generate right bidiagonalizing vectors in WORK(IU)
2543 * (Workspace: need 2*M*M+4*M-1,
2544 * prefer 2*M*M+3*M+(M-1)*NB)
2545 *
2546  CALL dorgbr( 'P', m, m, m, work( iu ), ldwrku,
2547  $ work( itaup ), work( iwork ),
2548  $ lwork-iwork+1, ierr )
2549 *
2550 * Generate left bidiagonalizing vectors in WORK(IR)
2551 * (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
2552 *
2553  CALL dorgbr( 'Q', m, m, m, work( ir ), ldwrkr,
2554  $ work( itauq ), work( iwork ),
2555  $ lwork-iwork+1, ierr )
2556  iwork = ie + m
2557 *
2558 * Perform bidiagonal QR iteration, computing left
2559 * singular vectors of L in WORK(IR) and computing
2560 * right singular vectors of L in WORK(IU)
2561 * (Workspace: need 2*M*M+BDSPAC)
2562 *
2563  CALL dbdsqr( 'U', m, m, m, 0, s, work( ie ),
2564  $ work( iu ), ldwrku, work( ir ),
2565  $ ldwrkr, dum, 1, work( iwork ), info )
2566 *
2567 * Multiply right singular vectors of L in WORK(IU) by
2568 * Q in A, storing result in VT
2569 * (Workspace: need M*M)
2570 *
2571  CALL dgemm( 'N', 'N', m, n, m, one, work( iu ),
2572  $ ldwrku, a, lda, zero, vt, ldvt )
2573 *
2574 * Copy left singular vectors of L to A
2575 * (Workspace: need M*M)
2576 *
2577  CALL dlacpy( 'F', m, m, work( ir ), ldwrkr, a,
2578  $ lda )
2579 *
2580  ELSE
2581 *
2582 * Insufficient workspace for a fast algorithm
2583 *
2584  itau = 1
2585  iwork = itau + m
2586 *
2587 * Compute A=L*Q, copying result to VT
2588 * (Workspace: need 2*M, prefer M+M*NB)
2589 *
2590  CALL dgelqf( m, n, a, lda, work( itau ),
2591  $ work( iwork ), lwork-iwork+1, ierr )
2592  CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
2593 *
2594 * Generate Q in VT
2595 * (Workspace: need 2*M, prefer M+M*NB)
2596 *
2597  CALL dorglq( m, n, m, vt, ldvt, work( itau ),
2598  $ work( iwork ), lwork-iwork+1, ierr )
2599  ie = itau
2600  itauq = ie + m
2601  itaup = itauq + m
2602  iwork = itaup + m
2603 *
2604 * Zero out above L in A
2605 *
2606  CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ),
2607  $ lda )
2608 *
2609 * Bidiagonalize L in A
2610 * (Workspace: need 4*M, prefer 3*M+2*M*NB)
2611 *
2612  CALL dgebrd( m, m, a, lda, s, work( ie ),
2613  $ work( itauq ), work( itaup ),
2614  $ work( iwork ), lwork-iwork+1, ierr )
2615 *
2616 * Multiply right vectors bidiagonalizing L by Q in VT
2617 * (Workspace: need 3*M+N, prefer 3*M+N*NB)
2618 *
2619  CALL dormbr( 'P', 'L', 'T', m, n, m, a, lda,
2620  $ work( itaup ), vt, ldvt,
2621  $ work( iwork ), lwork-iwork+1, ierr )
2622 *
2623 * Generate left bidiagonalizing vectors of L in A
2624 * (Workspace: need 4*M, prefer 3*M+M*NB)
2625 *
2626  CALL dorgbr( 'Q', m, m, m, a, lda, work( itauq ),
2627  $ work( iwork ), lwork-iwork+1, ierr )
2628  iwork = ie + m
2629 *
2630 * Perform bidiagonal QR iteration, compute left
2631 * singular vectors of A in A and compute right
2632 * singular vectors of A in VT
2633 * (Workspace: need BDSPAC)
2634 *
2635  CALL dbdsqr( 'U', m, n, m, 0, s, work( ie ), vt,
2636  $ ldvt, a, lda, dum, 1, work( iwork ),
2637  $ info )
2638 *
2639  END IF
2640 *
2641  ELSE IF( wntuas ) THEN
2642 *
2643 * Path 6t(N much larger than M, JOBU='S' or 'A',
2644 * JOBVT='S')
2645 * M right singular vectors to be computed in VT and
2646 * M left singular vectors to be computed in U
2647 *
2648  IF( lwork.GE.m*m+max( 4*m, bdspac ) ) THEN
2649 *
2650 * Sufficient workspace for a fast algorithm
2651 *
2652  iu = 1
2653  IF( lwork.GE.wrkbl+lda*m ) THEN
2654 *
2655 * WORK(IU) is LDA by N
2656 *
2657  ldwrku = lda
2658  ELSE
2659 *
2660 * WORK(IU) is LDA by M
2661 *
2662  ldwrku = m
2663  END IF
2664  itau = iu + ldwrku*m
2665  iwork = itau + m
2666 *
2667 * Compute A=L*Q
2668 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2669 *
2670  CALL dgelqf( m, n, a, lda, work( itau ),
2671  $ work( iwork ), lwork-iwork+1, ierr )
2672 *
2673 * Copy L to WORK(IU), zeroing out above it
2674 *
2675  CALL dlacpy( 'L', m, m, a, lda, work( iu ),
2676  $ ldwrku )
2677  CALL dlaset( 'U', m-1, m-1, zero, zero,
2678  $ work( iu+ldwrku ), ldwrku )
2679 *
2680 * Generate Q in A
2681 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2682 *
2683  CALL dorglq( m, n, m, a, lda, work( itau ),
2684  $ work( iwork ), lwork-iwork+1, ierr )
2685  ie = itau
2686  itauq = ie + m
2687  itaup = itauq + m
2688  iwork = itaup + m
2689 *
2690 * Bidiagonalize L in WORK(IU), copying result to U
2691 * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
2692 *
2693  CALL dgebrd( m, m, work( iu ), ldwrku, s,
2694  $ work( ie ), work( itauq ),
2695  $ work( itaup ), work( iwork ),
2696  $ lwork-iwork+1, ierr )
2697  CALL dlacpy( 'L', m, m, work( iu ), ldwrku, u,
2698  $ ldu )
2699 *
2700 * Generate right bidiagonalizing vectors in WORK(IU)
2701 * (Workspace: need M*M+4*M-1,
2702 * prefer M*M+3*M+(M-1)*NB)
2703 *
2704  CALL dorgbr( 'P', m, m, m, work( iu ), ldwrku,
2705  $ work( itaup ), work( iwork ),
2706  $ lwork-iwork+1, ierr )
2707 *
2708 * Generate left bidiagonalizing vectors in U
2709 * (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
2710 *
2711  CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
2712  $ work( iwork ), lwork-iwork+1, ierr )
2713  iwork = ie + m
2714 *
2715 * Perform bidiagonal QR iteration, computing left
2716 * singular vectors of L in U and computing right
2717 * singular vectors of L in WORK(IU)
2718 * (Workspace: need M*M+BDSPAC)
2719 *
2720  CALL dbdsqr( 'U', m, m, m, 0, s, work( ie ),
2721  $ work( iu ), ldwrku, u, ldu, dum, 1,
2722  $ work( iwork ), info )
2723 *
2724 * Multiply right singular vectors of L in WORK(IU) by
2725 * Q in A, storing result in VT
2726 * (Workspace: need M*M)
2727 *
2728  CALL dgemm( 'N', 'N', m, n, m, one, work( iu ),
2729  $ ldwrku, a, lda, zero, vt, ldvt )
2730 *
2731  ELSE
2732 *
2733 * Insufficient workspace for a fast algorithm
2734 *
2735  itau = 1
2736  iwork = itau + m
2737 *
2738 * Compute A=L*Q, copying result to VT
2739 * (Workspace: need 2*M, prefer M+M*NB)
2740 *
2741  CALL dgelqf( m, n, a, lda, work( itau ),
2742  $ work( iwork ), lwork-iwork+1, ierr )
2743  CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
2744 *
2745 * Generate Q in VT
2746 * (Workspace: need 2*M, prefer M+M*NB)
2747 *
2748  CALL dorglq( m, n, m, vt, ldvt, work( itau ),
2749  $ work( iwork ), lwork-iwork+1, ierr )
2750 *
2751 * Copy L to U, zeroing out above it
2752 *
2753  CALL dlacpy( 'L', m, m, a, lda, u, ldu )
2754  CALL dlaset( 'U', m-1, m-1, zero, zero, u( 1, 2 ),
2755  $ ldu )
2756  ie = itau
2757  itauq = ie + m
2758  itaup = itauq + m
2759  iwork = itaup + m
2760 *
2761 * Bidiagonalize L in U
2762 * (Workspace: need 4*M, prefer 3*M+2*M*NB)
2763 *
2764  CALL dgebrd( m, m, u, ldu, s, work( ie ),
2765  $ work( itauq ), work( itaup ),
2766  $ work( iwork ), lwork-iwork+1, ierr )
2767 *
2768 * Multiply right bidiagonalizing vectors in U by Q
2769 * in VT
2770 * (Workspace: need 3*M+N, prefer 3*M+N*NB)
2771 *
2772  CALL dormbr( 'P', 'L', 'T', m, n, m, u, ldu,
2773  $ work( itaup ), vt, ldvt,
2774  $ work( iwork ), lwork-iwork+1, ierr )
2775 *
2776 * Generate left bidiagonalizing vectors in U
2777 * (Workspace: need 4*M, prefer 3*M+M*NB)
2778 *
2779  CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
2780  $ work( iwork ), lwork-iwork+1, ierr )
2781  iwork = ie + m
2782 *
2783 * Perform bidiagonal QR iteration, computing left
2784 * singular vectors of A in U and computing right
2785 * singular vectors of A in VT
2786 * (Workspace: need BDSPAC)
2787 *
2788  CALL dbdsqr( 'U', m, n, m, 0, s, work( ie ), vt,
2789  $ ldvt, u, ldu, dum, 1, work( iwork ),
2790  $ info )
2791 *
2792  END IF
2793 *
2794  END IF
2795 *
2796  ELSE IF( wntva ) THEN
2797 *
2798  IF( wntun ) THEN
2799 *
2800 * Path 7t(N much larger than M, JOBU='N', JOBVT='A')
2801 * N right singular vectors to be computed in VT and
2802 * no left singular vectors to be computed
2803 *
2804  IF( lwork.GE.m*m+max( n+m, 4*m, bdspac ) ) THEN
2805 *
2806 * Sufficient workspace for a fast algorithm
2807 *
2808  ir = 1
2809  IF( lwork.GE.wrkbl+lda*m ) THEN
2810 *
2811 * WORK(IR) is LDA by M
2812 *
2813  ldwrkr = lda
2814  ELSE
2815 *
2816 * WORK(IR) is M by M
2817 *
2818  ldwrkr = m
2819  END IF
2820  itau = ir + ldwrkr*m
2821  iwork = itau + m
2822 *
2823 * Compute A=L*Q, copying result to VT
2824 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2825 *
2826  CALL dgelqf( m, n, a, lda, work( itau ),
2827  $ work( iwork ), lwork-iwork+1, ierr )
2828  CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
2829 *
2830 * Copy L to WORK(IR), zeroing out above it
2831 *
2832  CALL dlacpy( 'L', m, m, a, lda, work( ir ),
2833  $ ldwrkr )
2834  CALL dlaset( 'U', m-1, m-1, zero, zero,
2835  $ work( ir+ldwrkr ), ldwrkr )
2836 *
2837 * Generate Q in VT
2838 * (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
2839 *
2840  CALL dorglq( n, n, m, vt, ldvt, work( itau ),
2841  $ work( iwork ), lwork-iwork+1, ierr )
2842  ie = itau
2843  itauq = ie + m
2844  itaup = itauq + m
2845  iwork = itaup + m
2846 *
2847 * Bidiagonalize L in WORK(IR)
2848 * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
2849 *
2850  CALL dgebrd( m, m, work( ir ), ldwrkr, s,
2851  $ work( ie ), work( itauq ),
2852  $ work( itaup ), work( iwork ),
2853  $ lwork-iwork+1, ierr )
2854 *
2855 * Generate right bidiagonalizing vectors in WORK(IR)
2856 * (Workspace: need M*M+4*M-1,
2857 * prefer M*M+3*M+(M-1)*NB)
2858 *
2859  CALL dorgbr( 'P', m, m, m, work( ir ), ldwrkr,
2860  $ work( itaup ), work( iwork ),
2861  $ lwork-iwork+1, ierr )
2862  iwork = ie + m
2863 *
2864 * Perform bidiagonal QR iteration, computing right
2865 * singular vectors of L in WORK(IR)
2866 * (Workspace: need M*M+BDSPAC)
2867 *
2868  CALL dbdsqr( 'U', m, m, 0, 0, s, work( ie ),
2869  $ work( ir ), ldwrkr, dum, 1, dum, 1,
2870  $ work( iwork ), info )
2871 *
2872 * Multiply right singular vectors of L in WORK(IR) by
2873 * Q in VT, storing result in A
2874 * (Workspace: need M*M)
2875 *
2876  CALL dgemm( 'N', 'N', m, n, m, one, work( ir ),
2877  $ ldwrkr, vt, ldvt, zero, a, lda )
2878 *
2879 * Copy right singular vectors of A from A to VT
2880 *
2881  CALL dlacpy( 'F', m, n, a, lda, vt, ldvt )
2882 *
2883  ELSE
2884 *
2885 * Insufficient workspace for a fast algorithm
2886 *
2887  itau = 1
2888  iwork = itau + m
2889 *
2890 * Compute A=L*Q, copying result to VT
2891 * (Workspace: need 2*M, prefer M+M*NB)
2892 *
2893  CALL dgelqf( m, n, a, lda, work( itau ),
2894  $ work( iwork ), lwork-iwork+1, ierr )
2895  CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
2896 *
2897 * Generate Q in VT
2898 * (Workspace: need M+N, prefer M+N*NB)
2899 *
2900  CALL dorglq( n, n, m, vt, ldvt, work( itau ),
2901  $ work( iwork ), lwork-iwork+1, ierr )
2902  ie = itau
2903  itauq = ie + m
2904  itaup = itauq + m
2905  iwork = itaup + m
2906 *
2907 * Zero out above L in A
2908 *
2909  CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ),
2910  $ lda )
2911 *
2912 * Bidiagonalize L in A
2913 * (Workspace: need 4*M, prefer 3*M+2*M*NB)
2914 *
2915  CALL dgebrd( m, m, a, lda, s, work( ie ),
2916  $ work( itauq ), work( itaup ),
2917  $ work( iwork ), lwork-iwork+1, ierr )
2918 *
2919 * Multiply right bidiagonalizing vectors in A by Q
2920 * in VT
2921 * (Workspace: need 3*M+N, prefer 3*M+N*NB)
2922 *
2923  CALL dormbr( 'P', 'L', 'T', m, n, m, a, lda,
2924  $ work( itaup ), vt, ldvt,
2925  $ work( iwork ), lwork-iwork+1, ierr )
2926  iwork = ie + m
2927 *
2928 * Perform bidiagonal QR iteration, computing right
2929 * singular vectors of A in VT
2930 * (Workspace: need BDSPAC)
2931 *
2932  CALL dbdsqr( 'U', m, n, 0, 0, s, work( ie ), vt,
2933  $ ldvt, dum, 1, dum, 1, work( iwork ),
2934  $ info )
2935 *
2936  END IF
2937 *
2938  ELSE IF( wntuo ) THEN
2939 *
2940 * Path 8t(N much larger than M, JOBU='O', JOBVT='A')
2941 * N right singular vectors to be computed in VT and
2942 * M left singular vectors to be overwritten on A
2943 *
2944  IF( lwork.GE.2*m*m+max( n+m, 4*m, bdspac ) ) THEN
2945 *
2946 * Sufficient workspace for a fast algorithm
2947 *
2948  iu = 1
2949  IF( lwork.GE.wrkbl+2*lda*m ) THEN
2950 *
2951 * WORK(IU) is LDA by M and WORK(IR) is LDA by M
2952 *
2953  ldwrku = lda
2954  ir = iu + ldwrku*m
2955  ldwrkr = lda
2956  ELSE IF( lwork.GE.wrkbl+( lda+m )*m ) THEN
2957 *
2958 * WORK(IU) is LDA by M and WORK(IR) is M by M
2959 *
2960  ldwrku = lda
2961  ir = iu + ldwrku*m
2962  ldwrkr = m
2963  ELSE
2964 *
2965 * WORK(IU) is M by M and WORK(IR) is M by M
2966 *
2967  ldwrku = m
2968  ir = iu + ldwrku*m
2969  ldwrkr = m
2970  END IF
2971  itau = ir + ldwrkr*m
2972  iwork = itau + m
2973 *
2974 * Compute A=L*Q, copying result to VT
2975 * (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
2976 *
2977  CALL dgelqf( m, n, a, lda, work( itau ),
2978  $ work( iwork ), lwork-iwork+1, ierr )
2979  CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
2980 *
2981 * Generate Q in VT
2982 * (Workspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
2983 *
2984  CALL dorglq( n, n, m, vt, ldvt, work( itau ),
2985  $ work( iwork ), lwork-iwork+1, ierr )
2986 *
2987 * Copy L to WORK(IU), zeroing out above it
2988 *
2989  CALL dlacpy( 'L', m, m, a, lda, work( iu ),
2990  $ ldwrku )
2991  CALL dlaset( 'U', m-1, m-1, zero, zero,
2992  $ work( iu+ldwrku ), ldwrku )
2993  ie = itau
2994  itauq = ie + m
2995  itaup = itauq + m
2996  iwork = itaup + m
2997 *
2998 * Bidiagonalize L in WORK(IU), copying result to
2999 * WORK(IR)
3000 * (Workspace: need 2*M*M+4*M,
3001 * prefer 2*M*M+3*M+2*M*NB)
3002 *
3003  CALL dgebrd( m, m, work( iu ), ldwrku, s,
3004  $ work( ie ), work( itauq ),
3005  $ work( itaup ), work( iwork ),
3006  $ lwork-iwork+1, ierr )
3007  CALL dlacpy( 'L', m, m, work( iu ), ldwrku,
3008  $ work( ir ), ldwrkr )
3009 *
3010 * Generate right bidiagonalizing vectors in WORK(IU)
3011 * (Workspace: need 2*M*M+4*M-1,
3012 * prefer 2*M*M+3*M+(M-1)*NB)
3013 *
3014  CALL dorgbr( 'P', m, m, m, work( iu ), ldwrku,
3015  $ work( itaup ), work( iwork ),
3016  $ lwork-iwork+1, ierr )
3017 *
3018 * Generate left bidiagonalizing vectors in WORK(IR)
3019 * (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
3020 *
3021  CALL dorgbr( 'Q', m, m, m, work( ir ), ldwrkr,
3022  $ work( itauq ), work( iwork ),
3023  $ lwork-iwork+1, ierr )
3024  iwork = ie + m
3025 *
3026 * Perform bidiagonal QR iteration, computing left
3027 * singular vectors of L in WORK(IR) and computing
3028 * right singular vectors of L in WORK(IU)
3029 * (Workspace: need 2*M*M+BDSPAC)
3030 *
3031  CALL dbdsqr( 'U', m, m, m, 0, s, work( ie ),
3032  $ work( iu ), ldwrku, work( ir ),
3033  $ ldwrkr, dum, 1, work( iwork ), info )
3034 *
3035 * Multiply right singular vectors of L in WORK(IU) by
3036 * Q in VT, storing result in A
3037 * (Workspace: need M*M)
3038 *
3039  CALL dgemm( 'N', 'N', m, n, m, one, work( iu ),
3040  $ ldwrku, vt, ldvt, zero, a, lda )
3041 *
3042 * Copy right singular vectors of A from A to VT
3043 *
3044  CALL dlacpy( 'F', m, n, a, lda, vt, ldvt )
3045 *
3046 * Copy left singular vectors of A from WORK(IR) to A
3047 *
3048  CALL dlacpy( 'F', m, m, work( ir ), ldwrkr, a,
3049  $ lda )
3050 *
3051  ELSE
3052 *
3053 * Insufficient workspace for a fast algorithm
3054 *
3055  itau = 1
3056  iwork = itau + m
3057 *
3058 * Compute A=L*Q, copying result to VT
3059 * (Workspace: need 2*M, prefer M+M*NB)
3060 *
3061  CALL dgelqf( m, n, a, lda, work( itau ),
3062  $ work( iwork ), lwork-iwork+1, ierr )
3063  CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
3064 *
3065 * Generate Q in VT
3066 * (Workspace: need M+N, prefer M+N*NB)
3067 *
3068  CALL dorglq( n, n, m, vt, ldvt, work( itau ),
3069  $ work( iwork ), lwork-iwork+1, ierr )
3070  ie = itau
3071  itauq = ie + m
3072  itaup = itauq + m
3073  iwork = itaup + m
3074 *
3075 * Zero out above L in A
3076 *
3077  CALL dlaset( 'U', m-1, m-1, zero, zero, a( 1, 2 ),
3078  $ lda )
3079 *
3080 * Bidiagonalize L in A
3081 * (Workspace: need 4*M, prefer 3*M+2*M*NB)
3082 *
3083  CALL dgebrd( m, m, a, lda, s, work( ie ),
3084  $ work( itauq ), work( itaup ),
3085  $ work( iwork ), lwork-iwork+1, ierr )
3086 *
3087 * Multiply right bidiagonalizing vectors in A by Q
3088 * in VT
3089 * (Workspace: need 3*M+N, prefer 3*M+N*NB)
3090 *
3091  CALL dormbr( 'P', 'L', 'T', m, n, m, a, lda,
3092  $ work( itaup ), vt, ldvt,
3093  $ work( iwork ), lwork-iwork+1, ierr )
3094 *
3095 * Generate left bidiagonalizing vectors in A
3096 * (Workspace: need 4*M, prefer 3*M+M*NB)
3097 *
3098  CALL dorgbr( 'Q', m, m, m, a, lda, work( itauq ),
3099  $ work( iwork ), lwork-iwork+1, ierr )
3100  iwork = ie + m
3101 *
3102 * Perform bidiagonal QR iteration, computing left
3103 * singular vectors of A in A and computing right
3104 * singular vectors of A in VT
3105 * (Workspace: need BDSPAC)
3106 *
3107  CALL dbdsqr( 'U', m, n, m, 0, s, work( ie ), vt,
3108  $ ldvt, a, lda, dum, 1, work( iwork ),
3109  $ info )
3110 *
3111  END IF
3112 *
3113  ELSE IF( wntuas ) THEN
3114 *
3115 * Path 9t(N much larger than M, JOBU='S' or 'A',
3116 * JOBVT='A')
3117 * N right singular vectors to be computed in VT and
3118 * M left singular vectors to be computed in U
3119 *
3120  IF( lwork.GE.m*m+max( n+m, 4*m, bdspac ) ) THEN
3121 *
3122 * Sufficient workspace for a fast algorithm
3123 *
3124  iu = 1
3125  IF( lwork.GE.wrkbl+lda*m ) THEN
3126 *
3127 * WORK(IU) is LDA by M
3128 *
3129  ldwrku = lda
3130  ELSE
3131 *
3132 * WORK(IU) is M by M
3133 *
3134  ldwrku = m
3135  END IF
3136  itau = iu + ldwrku*m
3137  iwork = itau + m
3138 *
3139 * Compute A=L*Q, copying result to VT
3140 * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
3141 *
3142  CALL dgelqf( m, n, a, lda, work( itau ),
3143  $ work( iwork ), lwork-iwork+1, ierr )
3144  CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
3145 *
3146 * Generate Q in VT
3147 * (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
3148 *
3149  CALL dorglq( n, n, m, vt, ldvt, work( itau ),
3150  $ work( iwork ), lwork-iwork+1, ierr )
3151 *
3152 * Copy L to WORK(IU), zeroing out above it
3153 *
3154  CALL dlacpy( 'L', m, m, a, lda, work( iu ),
3155  $ ldwrku )
3156  CALL dlaset( 'U', m-1, m-1, zero, zero,
3157  $ work( iu+ldwrku ), ldwrku )
3158  ie = itau
3159  itauq = ie + m
3160  itaup = itauq + m
3161  iwork = itaup + m
3162 *
3163 * Bidiagonalize L in WORK(IU), copying result to U
3164 * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
3165 *
3166  CALL dgebrd( m, m, work( iu ), ldwrku, s,
3167  $ work( ie ), work( itauq ),
3168  $ work( itaup ), work( iwork ),
3169  $ lwork-iwork+1, ierr )
3170  CALL dlacpy( 'L', m, m, work( iu ), ldwrku, u,
3171  $ ldu )
3172 *
3173 * Generate right bidiagonalizing vectors in WORK(IU)
3174 * (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
3175 *
3176  CALL dorgbr( 'P', m, m, m, work( iu ), ldwrku,
3177  $ work( itaup ), work( iwork ),
3178  $ lwork-iwork+1, ierr )
3179 *
3180 * Generate left bidiagonalizing vectors in U
3181 * (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
3182 *
3183  CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
3184  $ work( iwork ), lwork-iwork+1, ierr )
3185  iwork = ie + m
3186 *
3187 * Perform bidiagonal QR iteration, computing left
3188 * singular vectors of L in U and computing right
3189 * singular vectors of L in WORK(IU)
3190 * (Workspace: need M*M+BDSPAC)
3191 *
3192  CALL dbdsqr( 'U', m, m, m, 0, s, work( ie ),
3193  $ work( iu ), ldwrku, u, ldu, dum, 1,
3194  $ work( iwork ), info )
3195 *
3196 * Multiply right singular vectors of L in WORK(IU) by
3197 * Q in VT, storing result in A
3198 * (Workspace: need M*M)
3199 *
3200  CALL dgemm( 'N', 'N', m, n, m, one, work( iu ),
3201  $ ldwrku, vt, ldvt, zero, a, lda )
3202 *
3203 * Copy right singular vectors of A from A to VT
3204 *
3205  CALL dlacpy( 'F', m, n, a, lda, vt, ldvt )
3206 *
3207  ELSE
3208 *
3209 * Insufficient workspace for a fast algorithm
3210 *
3211  itau = 1
3212  iwork = itau + m
3213 *
3214 * Compute A=L*Q, copying result to VT
3215 * (Workspace: need 2*M, prefer M+M*NB)
3216 *
3217  CALL dgelqf( m, n, a, lda, work( itau ),
3218  $ work( iwork ), lwork-iwork+1, ierr )
3219  CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
3220 *
3221 * Generate Q in VT
3222 * (Workspace: need M+N, prefer M+N*NB)
3223 *
3224  CALL dorglq( n, n, m, vt, ldvt, work( itau ),
3225  $ work( iwork ), lwork-iwork+1, ierr )
3226 *
3227 * Copy L to U, zeroing out above it
3228 *
3229  CALL dlacpy( 'L', m, m, a, lda, u, ldu )
3230  CALL dlaset( 'U', m-1, m-1, zero, zero, u( 1, 2 ),
3231  $ ldu )
3232  ie = itau
3233  itauq = ie + m
3234  itaup = itauq + m
3235  iwork = itaup + m
3236 *
3237 * Bidiagonalize L in U
3238 * (Workspace: need 4*M, prefer 3*M+2*M*NB)
3239 *
3240  CALL dgebrd( m, m, u, ldu, s, work( ie ),
3241  $ work( itauq ), work( itaup ),
3242  $ work( iwork ), lwork-iwork+1, ierr )
3243 *
3244 * Multiply right bidiagonalizing vectors in U by Q
3245 * in VT
3246 * (Workspace: need 3*M+N, prefer 3*M+N*NB)
3247 *
3248  CALL dormbr( 'P', 'L', 'T', m, n, m, u, ldu,
3249  $ work( itaup ), vt, ldvt,
3250  $ work( iwork ), lwork-iwork+1, ierr )
3251 *
3252 * Generate left bidiagonalizing vectors in U
3253 * (Workspace: need 4*M, prefer 3*M+M*NB)
3254 *
3255  CALL dorgbr( 'Q', m, m, m, u, ldu, work( itauq ),
3256  $ work( iwork ), lwork-iwork+1, ierr )
3257  iwork = ie + m
3258 *
3259 * Perform bidiagonal QR iteration, computing left
3260 * singular vectors of A in U and computing right
3261 * singular vectors of A in VT
3262 * (Workspace: need BDSPAC)
3263 *
3264  CALL dbdsqr( 'U', m, n, m, 0, s, work( ie ), vt,
3265  $ ldvt, u, ldu, dum, 1, work( iwork ),
3266  $ info )
3267 *
3268  END IF
3269 *
3270  END IF
3271 *
3272  END IF
3273 *
3274  ELSE
3275 *
3276 * N .LT. MNTHR
3277 *
3278 * Path 10t(N greater than M, but not much larger)
3279 * Reduce to bidiagonal form without LQ decomposition
3280 *
3281  ie = 1
3282  itauq = ie + m
3283  itaup = itauq + m
3284  iwork = itaup + m
3285 *
3286 * Bidiagonalize A
3287 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
3288 *
3289  CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
3290  $ work( itaup ), work( iwork ), lwork-iwork+1,
3291  $ ierr )
3292  IF( wntuas ) THEN
3293 *
3294 * If left singular vectors desired in U, copy result to U
3295 * and generate left bidiagonalizing vectors in U
3296 * (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
3297 *
3298  CALL dlacpy( 'L', m, m, a, lda, u, ldu )
3299  CALL dorgbr( 'Q', m, m, n, u, ldu, work( itauq ),
3300  $ work( iwork ), lwork-iwork+1, ierr )
3301  END IF
3302  IF( wntvas ) THEN
3303 *
3304 * If right singular vectors desired in VT, copy result to
3305 * VT and generate right bidiagonalizing vectors in VT
3306 * (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB)
3307 *
3308  CALL dlacpy( 'U', m, n, a, lda, vt, ldvt )
3309  IF( wntva )
3310  $ nrvt = n
3311  IF( wntvs )
3312  $ nrvt = m
3313  CALL dorgbr( 'P', nrvt, n, m, vt, ldvt, work( itaup ),
3314  $ work( iwork ), lwork-iwork+1, ierr )
3315  END IF
3316  IF( wntuo ) THEN
3317 *
3318 * If left singular vectors desired in A, generate left
3319 * bidiagonalizing vectors in A
3320 * (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
3321 *
3322  CALL dorgbr( 'Q', m, m, n, a, lda, work( itauq ),
3323  $ work( iwork ), lwork-iwork+1, ierr )
3324  END IF
3325  IF( wntvo ) THEN
3326 *
3327 * If right singular vectors desired in A, generate right
3328 * bidiagonalizing vectors in A
3329 * (Workspace: need 4*M, prefer 3*M+M*NB)
3330 *
3331  CALL dorgbr( 'P', m, n, m, a, lda, work( itaup ),
3332  $ work( iwork ), lwork-iwork+1, ierr )
3333  END IF
3334  iwork = ie + m
3335  IF( wntuas .OR. wntuo )
3336  $ nru = m
3337  IF( wntun )
3338  $ nru = 0
3339  IF( wntvas .OR. wntvo )
3340  $ ncvt = n
3341  IF( wntvn )
3342  $ ncvt = 0
3343  IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) ) THEN
3344 *
3345 * Perform bidiagonal QR iteration, if desired, computing
3346 * left singular vectors in U and computing right singular
3347 * vectors in VT
3348 * (Workspace: need BDSPAC)
3349 *
3350  CALL dbdsqr( 'L', m, ncvt, nru, 0, s, work( ie ), vt,
3351  $ ldvt, u, ldu, dum, 1, work( iwork ), info )
3352  ELSE IF( ( .NOT.wntuo ) .AND. wntvo ) THEN
3353 *
3354 * Perform bidiagonal QR iteration, if desired, computing
3355 * left singular vectors in U and computing right singular
3356 * vectors in A
3357 * (Workspace: need BDSPAC)
3358 *
3359  CALL dbdsqr( 'L', m, ncvt, nru, 0, s, work( ie ), a, lda,
3360  $ u, ldu, dum, 1, work( iwork ), info )
3361  ELSE
3362 *
3363 * Perform bidiagonal QR iteration, if desired, computing
3364 * left singular vectors in A and computing right singular
3365 * vectors in VT
3366 * (Workspace: need BDSPAC)
3367 *
3368  CALL dbdsqr( 'L', m, ncvt, nru, 0, s, work( ie ), vt,
3369  $ ldvt, a, lda, dum, 1, work( iwork ), info )
3370  END IF
3371 *
3372  END IF
3373 *
3374  END IF
3375 *
3376 * If DBDSQR failed to converge, copy unconverged superdiagonals
3377 * to WORK( 2:MINMN )
3378 *
3379  IF( info.NE.0 ) THEN
3380  IF( ie.GT.2 ) THEN
3381  DO 50 i = 1, minmn - 1
3382  work( i+1 ) = work( i+ie-1 )
3383  50 CONTINUE
3384  END IF
3385  IF( ie.LT.2 ) THEN
3386  DO 60 i = minmn - 1, 1, -1
3387  work( i+1 ) = work( i+ie-1 )
3388  60 CONTINUE
3389  END IF
3390  END IF
3391 *
3392 * Undo scaling if necessary
3393 *
3394  IF( iscl.EQ.1 ) THEN
3395  IF( anrm.GT.bignum )
3396  $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3397  $ ierr )
3398  IF( info.NE.0 .AND. anrm.GT.bignum )
3399  $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn-1, 1, work( 2 ),
3400  $ minmn, ierr )
3401  IF( anrm.LT.smlnum )
3402  $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3403  $ ierr )
3404  IF( info.NE.0 .AND. anrm.LT.smlnum )
3405  $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn-1, 1, work( 2 ),
3406  $ minmn, ierr )
3407  END IF
3408 *
3409 * Return optimal workspace in WORK(1)
3410 *
3411  work( 1 ) = maxwrk
3412 *
3413  RETURN
3414 *
3415 * End of DGESVD
3416 *
3417  END
dbdsqr
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
Definition: dbdsqr.f:3
dgebrd
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
Definition: dgebrd.f:3
dgelqf
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
Definition: dgelqf.f:2
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
Definition: dgemm.f:3
dgeqrf
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
Definition: dgeqrf.f:2
dgesvd
subroutine dgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
Definition: dgesvd.f:3
dlacpy
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
Definition: dlacpy.f:2
dlascl
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
Definition: dlascl.f:2
dlaset
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
Definition: dlaset.f:2
dorgbr
subroutine dorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
Definition: dorgbr.f:2
dorglq
subroutine dorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
Definition: dorglq.f:2
dorgqr
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
Definition: dorgqr.f:2
dormbr
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
Definition: dormbr.f:3
xerbla
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2