LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
zgesdd.f
Go to the documentation of this file.
1 *> \brief \b ZGESDD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZGESDD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesdd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesdd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesdd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
22 * WORK, LWORK, RWORK, IWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBZ
26 * INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IWORK( * )
30 * DOUBLE PRECISION RWORK( * ), S( * )
31 * COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
32 * $ WORK( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> ZGESDD computes the singular value decomposition (SVD) of a complex
42 *> M-by-N matrix A, optionally computing the left and/or right singular
43 *> vectors, by using divide-and-conquer method. The SVD is written
44 *>
45 *> A = U * SIGMA * conjugate-transpose(V)
46 *>
47 *> where SIGMA is an M-by-N matrix which is zero except for its
48 *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
49 *> V is an N-by-N unitary matrix. The diagonal elements of SIGMA
50 *> are the singular values of A; they are real and non-negative, and
51 *> are returned in descending order. The first min(m,n) columns of
52 *> U and V are the left and right singular vectors of A.
53 *>
54 *> Note that the routine returns VT = V**H, not V.
55 *>
56 *> The divide and conquer algorithm makes very mild assumptions about
57 *> floating point arithmetic. It will work on machines with a guard
58 *> digit in add/subtract, or on those binary machines without guard
59 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
60 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
61 *> without guard digits, but we know of none.
62 *> \endverbatim
63 *
64 * Arguments:
65 * ==========
66 *
67 *> \param[in] JOBZ
68 *> \verbatim
69 *> JOBZ is CHARACTER*1
70 *> Specifies options for computing all or part of the matrix U:
71 *> = 'A': all M columns of U and all N rows of V**H are
72 *> returned in the arrays U and VT;
73 *> = 'S': the first min(M,N) columns of U and the first
74 *> min(M,N) rows of V**H are returned in the arrays U
75 *> and VT;
76 *> = 'O': If M >= N, the first N columns of U are overwritten
77 *> in the array A and all rows of V**H are returned in
78 *> the array VT;
79 *> otherwise, all columns of U are returned in the
80 *> array U and the first M rows of V**H are overwritten
81 *> in the array A;
82 *> = 'N': no columns of U or rows of V**H are computed.
83 *> \endverbatim
84 *>
85 *> \param[in] M
86 *> \verbatim
87 *> M is INTEGER
88 *> The number of rows of the input matrix A. M >= 0.
89 *> \endverbatim
90 *>
91 *> \param[in] N
92 *> \verbatim
93 *> N is INTEGER
94 *> The number of columns of the input matrix A. N >= 0.
95 *> \endverbatim
96 *>
97 *> \param[in,out] A
98 *> \verbatim
99 *> A is COMPLEX*16 array, dimension (LDA,N)
100 *> On entry, the M-by-N matrix A.
101 *> On exit,
102 *> if JOBZ = 'O', A is overwritten with the first N columns
103 *> of U (the left singular vectors, stored
104 *> columnwise) if M >= N;
105 *> A is overwritten with the first M rows
106 *> of V**H (the right singular vectors, stored
107 *> rowwise) otherwise.
108 *> if JOBZ .ne. 'O', the contents of A are destroyed.
109 *> \endverbatim
110 *>
111 *> \param[in] LDA
112 *> \verbatim
113 *> LDA is INTEGER
114 *> The leading dimension of the array A. LDA >= max(1,M).
115 *> \endverbatim
116 *>
117 *> \param[out] S
118 *> \verbatim
119 *> S is DOUBLE PRECISION array, dimension (min(M,N))
120 *> The singular values of A, sorted so that S(i) >= S(i+1).
121 *> \endverbatim
122 *>
123 *> \param[out] U
124 *> \verbatim
125 *> U is COMPLEX*16 array, dimension (LDU,UCOL)
126 *> UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
127 *> UCOL = min(M,N) if JOBZ = 'S'.
128 *> If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
129 *> unitary matrix U;
130 *> if JOBZ = 'S', U contains the first min(M,N) columns of U
131 *> (the left singular vectors, stored columnwise);
132 *> if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
133 *> \endverbatim
134 *>
135 *> \param[in] LDU
136 *> \verbatim
137 *> LDU is INTEGER
138 *> The leading dimension of the array U. LDU >= 1;
139 *> if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
140 *> \endverbatim
141 *>
142 *> \param[out] VT
143 *> \verbatim
144 *> VT is COMPLEX*16 array, dimension (LDVT,N)
145 *> If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
146 *> N-by-N unitary matrix V**H;
147 *> if JOBZ = 'S', VT contains the first min(M,N) rows of
148 *> V**H (the right singular vectors, stored rowwise);
149 *> if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
150 *> \endverbatim
151 *>
152 *> \param[in] LDVT
153 *> \verbatim
154 *> LDVT is INTEGER
155 *> The leading dimension of the array VT. LDVT >= 1;
156 *> if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
157 *> if JOBZ = 'S', LDVT >= min(M,N).
158 *> \endverbatim
159 *>
160 *> \param[out] WORK
161 *> \verbatim
162 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
163 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
164 *> \endverbatim
165 *>
166 *> \param[in] LWORK
167 *> \verbatim
168 *> LWORK is INTEGER
169 *> The dimension of the array WORK. LWORK >= 1.
170 *> If LWORK = -1, a workspace query is assumed. The optimal
171 *> size for the WORK array is calculated and stored in WORK(1),
172 *> and no other work except argument checking is performed.
173 *>
174 *> Let mx = max(M,N) and mn = min(M,N).
175 *> If JOBZ = 'N', LWORK >= 2*mn + mx.
176 *> If JOBZ = 'O', LWORK >= 2*mn*mn + 2*mn + mx.
177 *> If JOBZ = 'S', LWORK >= mn*mn + 3*mn.
178 *> If JOBZ = 'A', LWORK >= mn*mn + 2*mn + mx.
179 *> These are not tight minimums in all cases; see comments inside code.
180 *> For good performance, LWORK should generally be larger;
181 *> a query is recommended.
182 *> \endverbatim
183 *>
184 *> \param[out] RWORK
185 *> \verbatim
186 *> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
187 *> Let mx = max(M,N) and mn = min(M,N).
188 *> If JOBZ = 'N', LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn);
189 *> else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn;
190 *> else LRWORK >= max( 5*mn*mn + 5*mn,
191 *> 2*mx*mn + 2*mn*mn + mn ).
192 *> \endverbatim
193 *>
194 *> \param[out] IWORK
195 *> \verbatim
196 *> IWORK is INTEGER array, dimension (8*min(M,N))
197 *> \endverbatim
198 *>
199 *> \param[out] INFO
200 *> \verbatim
201 *> INFO is INTEGER
202 *> = 0: successful exit.
203 *> < 0: if INFO = -i, the i-th argument had an illegal value.
204 *> > 0: The updating process of DBDSDC did not converge.
205 *> \endverbatim
206 *
207 * Authors:
208 * ========
209 *
210 *> \author Univ. of Tennessee
211 *> \author Univ. of California Berkeley
212 *> \author Univ. of Colorado Denver
213 *> \author NAG Ltd.
214 *
215 *> \date June 2016
216 *
217 *> \ingroup complex16GEsing
218 *
219 *> \par Contributors:
220 * ==================
221 *>
222 *> Ming Gu and Huan Ren, Computer Science Division, University of
223 *> California at Berkeley, USA
224 *>
225 * =====================================================================
226  SUBROUTINE zgesdd( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
227  $ WORK, LWORK, RWORK, IWORK, INFO )
228  implicit none
229 *
230 * -- LAPACK driver routine (version 3.7.0) --
231 * -- LAPACK is a software package provided by Univ. of Tennessee, --
232 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
233 * June 2016
234 *
235 * .. Scalar Arguments ..
236  CHARACTER JOBZ
237  INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
238 * ..
239 * .. Array Arguments ..
240  INTEGER IWORK( * )
241  DOUBLE PRECISION RWORK( * ), S( * )
242  COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
243  $ work( * )
244 * ..
245 *
246 * =====================================================================
247 *
248 * .. Parameters ..
249  COMPLEX*16 CZERO, CONE
250  parameter( czero = ( 0.0d+0, 0.0d+0 ),
251  $ cone = ( 1.0d+0, 0.0d+0 ) )
252  DOUBLE PRECISION ZERO, ONE
253  parameter( zero = 0.0d+0, one = 1.0d+0 )
254 * ..
255 * .. Local Scalars ..
256  LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
257  INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
258  $ iscl, itau, itaup, itauq, iu, ivt, ldwkvt,
259  $ ldwrkl, ldwrkr, ldwrku, maxwrk, minmn, minwrk,
260  $ mnthr1, mnthr2, nrwork, nwork, wrkbl
261  INTEGER LWORK_ZGEBRD_MN, LWORK_ZGEBRD_MM,
262  $ lwork_zgebrd_nn, lwork_zgelqf_mn,
263  $ lwork_zgeqrf_mn,
264  $ lwork_zungbr_p_mn, lwork_zungbr_p_nn,
265  $ lwork_zungbr_q_mn, lwork_zungbr_q_mm,
266  $ lwork_zunglq_mn, lwork_zunglq_nn,
267  $ lwork_zungqr_mm, lwork_zungqr_mn,
268  $ lwork_zunmbr_prc_mm, lwork_zunmbr_qln_mm,
269  $ lwork_zunmbr_prc_mn, lwork_zunmbr_qln_mn,
270  $ lwork_zunmbr_prc_nn, lwork_zunmbr_qln_nn
271  DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
272 * ..
273 * .. Local Arrays ..
274  INTEGER IDUM( 1 )
275  DOUBLE PRECISION DUM( 1 )
276  COMPLEX*16 CDUM( 1 )
277 * ..
278 * .. External Subroutines ..
279  EXTERNAL dbdsdc, dlascl, xerbla, zgebrd, zgelqf, zgemm,
280  $ zgeqrf, zlacp2, zlacpy, zlacrm, zlarcm, zlascl,
281  $ zlaset, zungbr, zunglq, zungqr, zunmbr
282 * ..
283 * .. External Functions ..
284  LOGICAL LSAME
285  DOUBLE PRECISION DLAMCH, ZLANGE
286  EXTERNAL lsame, dlamch, zlange
287 * ..
288 * .. Intrinsic Functions ..
289  INTRINSIC int, max, min, sqrt
290 * ..
291 * .. Executable Statements ..
292 *
293 * Test the input arguments
294 *
295  info = 0
296  minmn = min( m, n )
297  mnthr1 = int( minmn*17.0d0 / 9.0d0 )
298  mnthr2 = int( minmn*5.0d0 / 3.0d0 )
299  wntqa = lsame( jobz, 'A' )
300  wntqs = lsame( jobz, 'S' )
301  wntqas = wntqa .OR. wntqs
302  wntqo = lsame( jobz, 'O' )
303  wntqn = lsame( jobz, 'N' )
304  lquery = ( lwork.EQ.-1 )
305  minwrk = 1
306  maxwrk = 1
307 *
308  IF( .NOT.( wntqa .OR. wntqs .OR. wntqo .OR. wntqn ) ) THEN
309  info = -1
310  ELSE IF( m.LT.0 ) THEN
311  info = -2
312  ELSE IF( n.LT.0 ) THEN
313  info = -3
314  ELSE IF( lda.LT.max( 1, m ) ) THEN
315  info = -5
316  ELSE IF( ldu.LT.1 .OR. ( wntqas .AND. ldu.LT.m ) .OR.
317  $ ( wntqo .AND. m.LT.n .AND. ldu.LT.m ) ) THEN
318  info = -8
319  ELSE IF( ldvt.LT.1 .OR. ( wntqa .AND. ldvt.LT.n ) .OR.
320  $ ( wntqs .AND. ldvt.LT.minmn ) .OR.
321  $ ( wntqo .AND. m.GE.n .AND. ldvt.LT.n ) ) THEN
322  info = -10
323  END IF
324 *
325 * Compute workspace
326 * Note: Comments in the code beginning "Workspace:" describe the
327 * minimal amount of workspace allocated at that point in the code,
328 * as well as the preferred amount for good performance.
329 * CWorkspace refers to complex workspace, and RWorkspace to
330 * real workspace. NB refers to the optimal block size for the
331 * immediately following subroutine, as returned by ILAENV.)
332 *
333  IF( info.EQ.0 ) THEN
334  minwrk = 1
335  maxwrk = 1
336  IF( m.GE.n .AND. minmn.GT.0 ) THEN
337 *
338 * There is no complex work space needed for bidiagonal SVD
339 * The real work space needed for bidiagonal SVD (dbdsdc) is
340 * BDSPAC = 3*N*N + 4*N for singular values and vectors;
341 * BDSPAC = 4*N for singular values only;
342 * not including e, RU, and RVT matrices.
343 *
344 * Compute space preferred for each routine
345  CALL zgebrd( m, n, cdum(1), m, dum(1), dum(1), cdum(1),
346  $ cdum(1), cdum(1), -1, ierr )
347  lwork_zgebrd_mn = int( cdum(1) )
348 *
349  CALL zgebrd( n, n, cdum(1), n, dum(1), dum(1), cdum(1),
350  $ cdum(1), cdum(1), -1, ierr )
351  lwork_zgebrd_nn = int( cdum(1) )
352 *
353  CALL zgeqrf( m, n, cdum(1), m, cdum(1), cdum(1), -1, ierr )
354  lwork_zgeqrf_mn = int( cdum(1) )
355 *
356  CALL zungbr( 'P', n, n, n, cdum(1), n, cdum(1), cdum(1),
357  $ -1, ierr )
358  lwork_zungbr_p_nn = int( cdum(1) )
359 *
360  CALL zungbr( 'Q', m, m, n, cdum(1), m, cdum(1), cdum(1),
361  $ -1, ierr )
362  lwork_zungbr_q_mm = int( cdum(1) )
363 *
364  CALL zungbr( 'Q', m, n, n, cdum(1), m, cdum(1), cdum(1),
365  $ -1, ierr )
366  lwork_zungbr_q_mn = int( cdum(1) )
367 *
368  CALL zungqr( m, m, n, cdum(1), m, cdum(1), cdum(1),
369  $ -1, ierr )
370  lwork_zungqr_mm = int( cdum(1) )
371 *
372  CALL zungqr( m, n, n, cdum(1), m, cdum(1), cdum(1),
373  $ -1, ierr )
374  lwork_zungqr_mn = int( cdum(1) )
375 *
376  CALL zunmbr( 'P', 'R', 'C', n, n, n, cdum(1), n, cdum(1),
377  $ cdum(1), n, cdum(1), -1, ierr )
378  lwork_zunmbr_prc_nn = int( cdum(1) )
379 *
380  CALL zunmbr( 'Q', 'L', 'N', m, m, n, cdum(1), m, cdum(1),
381  $ cdum(1), m, cdum(1), -1, ierr )
382  lwork_zunmbr_qln_mm = int( cdum(1) )
383 *
384  CALL zunmbr( 'Q', 'L', 'N', m, n, n, cdum(1), m, cdum(1),
385  $ cdum(1), m, cdum(1), -1, ierr )
386  lwork_zunmbr_qln_mn = int( cdum(1) )
387 *
388  CALL zunmbr( 'Q', 'L', 'N', n, n, n, cdum(1), n, cdum(1),
389  $ cdum(1), n, cdum(1), -1, ierr )
390  lwork_zunmbr_qln_nn = int( cdum(1) )
391 *
392  IF( m.GE.mnthr1 ) THEN
393  IF( wntqn ) THEN
394 *
395 * Path 1 (M >> N, JOBZ='N')
396 *
397  maxwrk = n + lwork_zgeqrf_mn
398  maxwrk = max( maxwrk, 2*n + lwork_zgebrd_nn )
399  minwrk = 3*n
400  ELSE IF( wntqo ) THEN
401 *
402 * Path 2 (M >> N, JOBZ='O')
403 *
404  wrkbl = n + lwork_zgeqrf_mn
405  wrkbl = max( wrkbl, n + lwork_zungqr_mn )
406  wrkbl = max( wrkbl, 2*n + lwork_zgebrd_nn )
407  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_qln_nn )
408  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_prc_nn )
409  maxwrk = m*n + n*n + wrkbl
410  minwrk = 2*n*n + 3*n
411  ELSE IF( wntqs ) THEN
412 *
413 * Path 3 (M >> N, JOBZ='S')
414 *
415  wrkbl = n + lwork_zgeqrf_mn
416  wrkbl = max( wrkbl, n + lwork_zungqr_mn )
417  wrkbl = max( wrkbl, 2*n + lwork_zgebrd_nn )
418  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_qln_nn )
419  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_prc_nn )
420  maxwrk = n*n + wrkbl
421  minwrk = n*n + 3*n
422  ELSE IF( wntqa ) THEN
423 *
424 * Path 4 (M >> N, JOBZ='A')
425 *
426  wrkbl = n + lwork_zgeqrf_mn
427  wrkbl = max( wrkbl, n + lwork_zungqr_mm )
428  wrkbl = max( wrkbl, 2*n + lwork_zgebrd_nn )
429  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_qln_nn )
430  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_prc_nn )
431  maxwrk = n*n + wrkbl
432  minwrk = n*n + max( 3*n, n + m )
433  END IF
434  ELSE IF( m.GE.mnthr2 ) THEN
435 *
436 * Path 5 (M >> N, but not as much as MNTHR1)
437 *
438  maxwrk = 2*n + lwork_zgebrd_mn
439  minwrk = 2*n + m
440  IF( wntqo ) THEN
441 * Path 5o (M >> N, JOBZ='O')
442  maxwrk = max( maxwrk, 2*n + lwork_zungbr_p_nn )
443  maxwrk = max( maxwrk, 2*n + lwork_zungbr_q_mn )
444  maxwrk = maxwrk + m*n
445  minwrk = minwrk + n*n
446  ELSE IF( wntqs ) THEN
447 * Path 5s (M >> N, JOBZ='S')
448  maxwrk = max( maxwrk, 2*n + lwork_zungbr_p_nn )
449  maxwrk = max( maxwrk, 2*n + lwork_zungbr_q_mn )
450  ELSE IF( wntqa ) THEN
451 * Path 5a (M >> N, JOBZ='A')
452  maxwrk = max( maxwrk, 2*n + lwork_zungbr_p_nn )
453  maxwrk = max( maxwrk, 2*n + lwork_zungbr_q_mm )
454  END IF
455  ELSE
456 *
457 * Path 6 (M >= N, but not much larger)
458 *
459  maxwrk = 2*n + lwork_zgebrd_mn
460  minwrk = 2*n + m
461  IF( wntqo ) THEN
462 * Path 6o (M >= N, JOBZ='O')
463  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_prc_nn )
464  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_qln_mn )
465  maxwrk = maxwrk + m*n
466  minwrk = minwrk + n*n
467  ELSE IF( wntqs ) THEN
468 * Path 6s (M >= N, JOBZ='S')
469  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_qln_mn )
470  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_prc_nn )
471  ELSE IF( wntqa ) THEN
472 * Path 6a (M >= N, JOBZ='A')
473  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_qln_mm )
474  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_prc_nn )
475  END IF
476  END IF
477  ELSE IF( minmn.GT.0 ) THEN
478 *
479 * There is no complex work space needed for bidiagonal SVD
480 * The real work space needed for bidiagonal SVD (dbdsdc) is
481 * BDSPAC = 3*M*M + 4*M for singular values and vectors;
482 * BDSPAC = 4*M for singular values only;
483 * not including e, RU, and RVT matrices.
484 *
485 * Compute space preferred for each routine
486  CALL zgebrd( m, n, cdum(1), m, dum(1), dum(1), cdum(1),
487  $ cdum(1), cdum(1), -1, ierr )
488  lwork_zgebrd_mn = int( cdum(1) )
489 *
490  CALL zgebrd( m, m, cdum(1), m, dum(1), dum(1), cdum(1),
491  $ cdum(1), cdum(1), -1, ierr )
492  lwork_zgebrd_mm = int( cdum(1) )
493 *
494  CALL zgelqf( m, n, cdum(1), m, cdum(1), cdum(1), -1, ierr )
495  lwork_zgelqf_mn = int( cdum(1) )
496 *
497  CALL zungbr( 'P', m, n, m, cdum(1), m, cdum(1), cdum(1),
498  $ -1, ierr )
499  lwork_zungbr_p_mn = int( cdum(1) )
500 *
501  CALL zungbr( 'P', n, n, m, cdum(1), n, cdum(1), cdum(1),
502  $ -1, ierr )
503  lwork_zungbr_p_nn = int( cdum(1) )
504 *
505  CALL zungbr( 'Q', m, m, n, cdum(1), m, cdum(1), cdum(1),
506  $ -1, ierr )
507  lwork_zungbr_q_mm = int( cdum(1) )
508 *
509  CALL zunglq( m, n, m, cdum(1), m, cdum(1), cdum(1),
510  $ -1, ierr )
511  lwork_zunglq_mn = int( cdum(1) )
512 *
513  CALL zunglq( n, n, m, cdum(1), n, cdum(1), cdum(1),
514  $ -1, ierr )
515  lwork_zunglq_nn = int( cdum(1) )
516 *
517  CALL zunmbr( 'P', 'R', 'C', m, m, m, cdum(1), m, cdum(1),
518  $ cdum(1), m, cdum(1), -1, ierr )
519  lwork_zunmbr_prc_mm = int( cdum(1) )
520 *
521  CALL zunmbr( 'P', 'R', 'C', m, n, m, cdum(1), m, cdum(1),
522  $ cdum(1), m, cdum(1), -1, ierr )
523  lwork_zunmbr_prc_mn = int( cdum(1) )
524 *
525  CALL zunmbr( 'P', 'R', 'C', n, n, m, cdum(1), n, cdum(1),
526  $ cdum(1), n, cdum(1), -1, ierr )
527  lwork_zunmbr_prc_nn = int( cdum(1) )
528 *
529  CALL zunmbr( 'Q', 'L', 'N', m, m, m, cdum(1), m, cdum(1),
530  $ cdum(1), m, cdum(1), -1, ierr )
531  lwork_zunmbr_qln_mm = int( cdum(1) )
532 *
533  IF( n.GE.mnthr1 ) THEN
534  IF( wntqn ) THEN
535 *
536 * Path 1t (N >> M, JOBZ='N')
537 *
538  maxwrk = m + lwork_zgelqf_mn
539  maxwrk = max( maxwrk, 2*m + lwork_zgebrd_mm )
540  minwrk = 3*m
541  ELSE IF( wntqo ) THEN
542 *
543 * Path 2t (N >> M, JOBZ='O')
544 *
545  wrkbl = m + lwork_zgelqf_mn
546  wrkbl = max( wrkbl, m + lwork_zunglq_mn )
547  wrkbl = max( wrkbl, 2*m + lwork_zgebrd_mm )
548  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_qln_mm )
549  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_prc_mm )
550  maxwrk = m*n + m*m + wrkbl
551  minwrk = 2*m*m + 3*m
552  ELSE IF( wntqs ) THEN
553 *
554 * Path 3t (N >> M, JOBZ='S')
555 *
556  wrkbl = m + lwork_zgelqf_mn
557  wrkbl = max( wrkbl, m + lwork_zunglq_mn )
558  wrkbl = max( wrkbl, 2*m + lwork_zgebrd_mm )
559  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_qln_mm )
560  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_prc_mm )
561  maxwrk = m*m + wrkbl
562  minwrk = m*m + 3*m
563  ELSE IF( wntqa ) THEN
564 *
565 * Path 4t (N >> M, JOBZ='A')
566 *
567  wrkbl = m + lwork_zgelqf_mn
568  wrkbl = max( wrkbl, m + lwork_zunglq_nn )
569  wrkbl = max( wrkbl, 2*m + lwork_zgebrd_mm )
570  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_qln_mm )
571  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_prc_mm )
572  maxwrk = m*m + wrkbl
573  minwrk = m*m + max( 3*m, m + n )
574  END IF
575  ELSE IF( n.GE.mnthr2 ) THEN
576 *
577 * Path 5t (N >> M, but not as much as MNTHR1)
578 *
579  maxwrk = 2*m + lwork_zgebrd_mn
580  minwrk = 2*m + n
581  IF( wntqo ) THEN
582 * Path 5to (N >> M, JOBZ='O')
583  maxwrk = max( maxwrk, 2*m + lwork_zungbr_q_mm )
584  maxwrk = max( maxwrk, 2*m + lwork_zungbr_p_mn )
585  maxwrk = maxwrk + m*n
586  minwrk = minwrk + m*m
587  ELSE IF( wntqs ) THEN
588 * Path 5ts (N >> M, JOBZ='S')
589  maxwrk = max( maxwrk, 2*m + lwork_zungbr_q_mm )
590  maxwrk = max( maxwrk, 2*m + lwork_zungbr_p_mn )
591  ELSE IF( wntqa ) THEN
592 * Path 5ta (N >> M, JOBZ='A')
593  maxwrk = max( maxwrk, 2*m + lwork_zungbr_q_mm )
594  maxwrk = max( maxwrk, 2*m + lwork_zungbr_p_nn )
595  END IF
596  ELSE
597 *
598 * Path 6t (N > M, but not much larger)
599 *
600  maxwrk = 2*m + lwork_zgebrd_mn
601  minwrk = 2*m + n
602  IF( wntqo ) THEN
603 * Path 6to (N > M, JOBZ='O')
604  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_qln_mm )
605  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_prc_mn )
606  maxwrk = maxwrk + m*n
607  minwrk = minwrk + m*m
608  ELSE IF( wntqs ) THEN
609 * Path 6ts (N > M, JOBZ='S')
610  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_qln_mm )
611  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_prc_mn )
612  ELSE IF( wntqa ) THEN
613 * Path 6ta (N > M, JOBZ='A')
614  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_qln_mm )
615  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_prc_nn )
616  END IF
617  END IF
618  END IF
619  maxwrk = max( maxwrk, minwrk )
620  END IF
621  IF( info.EQ.0 ) THEN
622  work( 1 ) = maxwrk
623  IF( lwork.LT.minwrk .AND. .NOT. lquery ) THEN
624  info = -12
625  END IF
626  END IF
627 *
628  IF( info.NE.0 ) THEN
629  CALL xerbla( 'ZGESDD', -info )
630  RETURN
631  ELSE IF( lquery ) THEN
632  RETURN
633  END IF
634 *
635 * Quick return if possible
636 *
637  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
638  RETURN
639  END IF
640 *
641 * Get machine constants
642 *
643  eps = dlamch( 'P' )
644  smlnum = sqrt( dlamch( 'S' ) ) / eps
645  bignum = one / smlnum
646 *
647 * Scale A if max element outside range [SMLNUM,BIGNUM]
648 *
649  anrm = zlange( 'M', m, n, a, lda, dum )
650  iscl = 0
651  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
652  iscl = 1
653  CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
654  ELSE IF( anrm.GT.bignum ) THEN
655  iscl = 1
656  CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
657  END IF
658 *
659  IF( m.GE.n ) THEN
660 *
661 * A has at least as many rows as columns. If A has sufficiently
662 * more rows than columns, first reduce using the QR
663 * decomposition (if sufficient workspace available)
664 *
665  IF( m.GE.mnthr1 ) THEN
666 *
667  IF( wntqn ) THEN
668 *
669 * Path 1 (M >> N, JOBZ='N')
670 * No singular vectors to be computed
671 *
672  itau = 1
673  nwork = itau + n
674 *
675 * Compute A=Q*R
676 * CWorkspace: need N [tau] + N [work]
677 * CWorkspace: prefer N [tau] + N*NB [work]
678 * RWorkspace: need 0
679 *
680  CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
681  $ lwork-nwork+1, ierr )
682 *
683 * Zero out below R
684 *
685  CALL zlaset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
686  $ lda )
687  ie = 1
688  itauq = 1
689  itaup = itauq + n
690  nwork = itaup + n
691 *
692 * Bidiagonalize R in A
693 * CWorkspace: need 2*N [tauq, taup] + N [work]
694 * CWorkspace: prefer 2*N [tauq, taup] + 2*N*NB [work]
695 * RWorkspace: need N [e]
696 *
697  CALL zgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
698  $ work( itaup ), work( nwork ), lwork-nwork+1,
699  $ ierr )
700  nrwork = ie + n
701 *
702 * Perform bidiagonal SVD, compute singular values only
703 * CWorkspace: need 0
704 * RWorkspace: need N [e] + BDSPAC
705 *
706  CALL dbdsdc( 'U', 'N', n, s, rwork( ie ), dum,1,dum,1,
707  $ dum, idum, rwork( nrwork ), iwork, info )
708 *
709  ELSE IF( wntqo ) THEN
710 *
711 * Path 2 (M >> N, JOBZ='O')
712 * N left singular vectors to be overwritten on A and
713 * N right singular vectors to be computed in VT
714 *
715  iu = 1
716 *
717 * WORK(IU) is N by N
718 *
719  ldwrku = n
720  ir = iu + ldwrku*n
721  IF( lwork .GE. m*n + n*n + 3*n ) THEN
722 *
723 * WORK(IR) is M by N
724 *
725  ldwrkr = m
726  ELSE
727  ldwrkr = ( lwork - n*n - 3*n ) / n
728  END IF
729  itau = ir + ldwrkr*n
730  nwork = itau + n
731 *
732 * Compute A=Q*R
733 * CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work]
734 * CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
735 * RWorkspace: need 0
736 *
737  CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
738  $ lwork-nwork+1, ierr )
739 *
740 * Copy R to WORK( IR ), zeroing out below it
741 *
742  CALL zlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
743  CALL zlaset( 'L', n-1, n-1, czero, czero, work( ir+1 ),
744  $ ldwrkr )
745 *
746 * Generate Q in A
747 * CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work]
748 * CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
749 * RWorkspace: need 0
750 *
751  CALL zungqr( m, n, n, a, lda, work( itau ),
752  $ work( nwork ), lwork-nwork+1, ierr )
753  ie = 1
754  itauq = itau
755  itaup = itauq + n
756  nwork = itaup + n
757 *
758 * Bidiagonalize R in WORK(IR)
759 * CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
760 * CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
761 * RWorkspace: need N [e]
762 *
763  CALL zgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
764  $ work( itauq ), work( itaup ), work( nwork ),
765  $ lwork-nwork+1, ierr )
766 *
767 * Perform bidiagonal SVD, computing left singular vectors
768 * of R in WORK(IRU) and computing right singular vectors
769 * of R in WORK(IRVT)
770 * CWorkspace: need 0
771 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
772 *
773  iru = ie + n
774  irvt = iru + n*n
775  nrwork = irvt + n*n
776  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
777  $ n, rwork( irvt ), n, dum, idum,
778  $ rwork( nrwork ), iwork, info )
779 *
780 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
781 * Overwrite WORK(IU) by the left singular vectors of R
782 * CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
783 * CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
784 * RWorkspace: need 0
785 *
786  CALL zlacp2( 'F', n, n, rwork( iru ), n, work( iu ),
787  $ ldwrku )
788  CALL zunmbr( 'Q', 'L', 'N', n, n, n, work( ir ), ldwrkr,
789  $ work( itauq ), work( iu ), ldwrku,
790  $ work( nwork ), lwork-nwork+1, ierr )
791 *
792 * Copy real matrix RWORK(IRVT) to complex matrix VT
793 * Overwrite VT by the right singular vectors of R
794 * CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
795 * CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
796 * RWorkspace: need 0
797 *
798  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
799  CALL zunmbr( 'P', 'R', 'C', n, n, n, work( ir ), ldwrkr,
800  $ work( itaup ), vt, ldvt, work( nwork ),
801  $ lwork-nwork+1, ierr )
802 *
803 * Multiply Q in A by left singular vectors of R in
804 * WORK(IU), storing result in WORK(IR) and copying to A
805 * CWorkspace: need N*N [U] + N*N [R]
806 * CWorkspace: prefer N*N [U] + M*N [R]
807 * RWorkspace: need 0
808 *
809  DO 10 i = 1, m, ldwrkr
810  chunk = min( m-i+1, ldwrkr )
811  CALL zgemm( 'N', 'N', chunk, n, n, cone, a( i, 1 ),
812  $ lda, work( iu ), ldwrku, czero,
813  $ work( ir ), ldwrkr )
814  CALL zlacpy( 'F', chunk, n, work( ir ), ldwrkr,
815  $ a( i, 1 ), lda )
816  10 CONTINUE
817 *
818  ELSE IF( wntqs ) THEN
819 *
820 * Path 3 (M >> N, JOBZ='S')
821 * N left singular vectors to be computed in U and
822 * N right singular vectors to be computed in VT
823 *
824  ir = 1
825 *
826 * WORK(IR) is N by N
827 *
828  ldwrkr = n
829  itau = ir + ldwrkr*n
830  nwork = itau + n
831 *
832 * Compute A=Q*R
833 * CWorkspace: need N*N [R] + N [tau] + N [work]
834 * CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
835 * RWorkspace: need 0
836 *
837  CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
838  $ lwork-nwork+1, ierr )
839 *
840 * Copy R to WORK(IR), zeroing out below it
841 *
842  CALL zlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
843  CALL zlaset( 'L', n-1, n-1, czero, czero, work( ir+1 ),
844  $ ldwrkr )
845 *
846 * Generate Q in A
847 * CWorkspace: need N*N [R] + N [tau] + N [work]
848 * CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
849 * RWorkspace: need 0
850 *
851  CALL zungqr( m, n, n, a, lda, work( itau ),
852  $ work( nwork ), lwork-nwork+1, ierr )
853  ie = 1
854  itauq = itau
855  itaup = itauq + n
856  nwork = itaup + n
857 *
858 * Bidiagonalize R in WORK(IR)
859 * CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
860 * CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
861 * RWorkspace: need N [e]
862 *
863  CALL zgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
864  $ work( itauq ), work( itaup ), work( nwork ),
865  $ lwork-nwork+1, ierr )
866 *
867 * Perform bidiagonal SVD, computing left singular vectors
868 * of bidiagonal matrix in RWORK(IRU) and computing right
869 * singular vectors of bidiagonal matrix in RWORK(IRVT)
870 * CWorkspace: need 0
871 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
872 *
873  iru = ie + n
874  irvt = iru + n*n
875  nrwork = irvt + n*n
876  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
877  $ n, rwork( irvt ), n, dum, idum,
878  $ rwork( nrwork ), iwork, info )
879 *
880 * Copy real matrix RWORK(IRU) to complex matrix U
881 * Overwrite U by left singular vectors of R
882 * CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
883 * CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
884 * RWorkspace: need 0
885 *
886  CALL zlacp2( 'F', n, n, rwork( iru ), n, u, ldu )
887  CALL zunmbr( 'Q', 'L', 'N', n, n, n, work( ir ), ldwrkr,
888  $ work( itauq ), u, ldu, work( nwork ),
889  $ lwork-nwork+1, ierr )
890 *
891 * Copy real matrix RWORK(IRVT) to complex matrix VT
892 * Overwrite VT by right singular vectors of R
893 * CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
894 * CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
895 * RWorkspace: need 0
896 *
897  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
898  CALL zunmbr( 'P', 'R', 'C', n, n, n, work( ir ), ldwrkr,
899  $ work( itaup ), vt, ldvt, work( nwork ),
900  $ lwork-nwork+1, ierr )
901 *
902 * Multiply Q in A by left singular vectors of R in
903 * WORK(IR), storing result in U
904 * CWorkspace: need N*N [R]
905 * RWorkspace: need 0
906 *
907  CALL zlacpy( 'F', n, n, u, ldu, work( ir ), ldwrkr )
908  CALL zgemm( 'N', 'N', m, n, n, cone, a, lda, work( ir ),
909  $ ldwrkr, czero, u, ldu )
910 *
911  ELSE IF( wntqa ) THEN
912 *
913 * Path 4 (M >> N, JOBZ='A')
914 * M left singular vectors to be computed in U and
915 * N right singular vectors to be computed in VT
916 *
917  iu = 1
918 *
919 * WORK(IU) is N by N
920 *
921  ldwrku = n
922  itau = iu + ldwrku*n
923  nwork = itau + n
924 *
925 * Compute A=Q*R, copying result to U
926 * CWorkspace: need N*N [U] + N [tau] + N [work]
927 * CWorkspace: prefer N*N [U] + N [tau] + N*NB [work]
928 * RWorkspace: need 0
929 *
930  CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
931  $ lwork-nwork+1, ierr )
932  CALL zlacpy( 'L', m, n, a, lda, u, ldu )
933 *
934 * Generate Q in U
935 * CWorkspace: need N*N [U] + N [tau] + M [work]
936 * CWorkspace: prefer N*N [U] + N [tau] + M*NB [work]
937 * RWorkspace: need 0
938 *
939  CALL zungqr( m, m, n, u, ldu, work( itau ),
940  $ work( nwork ), lwork-nwork+1, ierr )
941 *
942 * Produce R in A, zeroing out below it
943 *
944  CALL zlaset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
945  $ lda )
946  ie = 1
947  itauq = itau
948  itaup = itauq + n
949  nwork = itaup + n
950 *
951 * Bidiagonalize R in A
952 * CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
953 * CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + 2*N*NB [work]
954 * RWorkspace: need N [e]
955 *
956  CALL zgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
957  $ work( itaup ), work( nwork ), lwork-nwork+1,
958  $ ierr )
959  iru = ie + n
960  irvt = iru + n*n
961  nrwork = irvt + n*n
962 *
963 * Perform bidiagonal SVD, computing left singular vectors
964 * of bidiagonal matrix in RWORK(IRU) and computing right
965 * singular vectors of bidiagonal matrix in RWORK(IRVT)
966 * CWorkspace: need 0
967 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
968 *
969  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
970  $ n, rwork( irvt ), n, dum, idum,
971  $ rwork( nrwork ), iwork, info )
972 *
973 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
974 * Overwrite WORK(IU) by left singular vectors of R
975 * CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
976 * CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
977 * RWorkspace: need 0
978 *
979  CALL zlacp2( 'F', n, n, rwork( iru ), n, work( iu ),
980  $ ldwrku )
981  CALL zunmbr( 'Q', 'L', 'N', n, n, n, a, lda,
982  $ work( itauq ), work( iu ), ldwrku,
983  $ work( nwork ), lwork-nwork+1, ierr )
984 *
985 * Copy real matrix RWORK(IRVT) to complex matrix VT
986 * Overwrite VT by right singular vectors of R
987 * CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
988 * CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
989 * RWorkspace: need 0
990 *
991  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
992  CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
993  $ work( itaup ), vt, ldvt, work( nwork ),
994  $ lwork-nwork+1, ierr )
995 *
996 * Multiply Q in U by left singular vectors of R in
997 * WORK(IU), storing result in A
998 * CWorkspace: need N*N [U]
999 * RWorkspace: need 0
1000 *
1001  CALL zgemm( 'N', 'N', m, n, n, cone, u, ldu, work( iu ),
1002  $ ldwrku, czero, a, lda )
1003 *
1004 * Copy left singular vectors of A from A to U
1005 *
1006  CALL zlacpy( 'F', m, n, a, lda, u, ldu )
1007 *
1008  END IF
1009 *
1010  ELSE IF( m.GE.mnthr2 ) THEN
1011 *
1012 * MNTHR2 <= M < MNTHR1
1013 *
1014 * Path 5 (M >> N, but not as much as MNTHR1)
1015 * Reduce to bidiagonal form without QR decomposition, use
1016 * ZUNGBR and matrix multiplication to compute singular vectors
1017 *
1018  ie = 1
1019  nrwork = ie + n
1020  itauq = 1
1021  itaup = itauq + n
1022  nwork = itaup + n
1023 *
1024 * Bidiagonalize A
1025 * CWorkspace: need 2*N [tauq, taup] + M [work]
1026 * CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
1027 * RWorkspace: need N [e]
1028 *
1029  CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1030  $ work( itaup ), work( nwork ), lwork-nwork+1,
1031  $ ierr )
1032  IF( wntqn ) THEN
1033 *
1034 * Path 5n (M >> N, JOBZ='N')
1035 * Compute singular values only
1036 * CWorkspace: need 0
1037 * RWorkspace: need N [e] + BDSPAC
1038 *
1039  CALL dbdsdc( 'U', 'N', n, s, rwork( ie ), dum, 1,dum,1,
1040  $ dum, idum, rwork( nrwork ), iwork, info )
1041  ELSE IF( wntqo ) THEN
1042  iu = nwork
1043  iru = nrwork
1044  irvt = iru + n*n
1045  nrwork = irvt + n*n
1046 *
1047 * Path 5o (M >> N, JOBZ='O')
1048 * Copy A to VT, generate P**H
1049 * CWorkspace: need 2*N [tauq, taup] + N [work]
1050 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1051 * RWorkspace: need 0
1052 *
1053  CALL zlacpy( 'U', n, n, a, lda, vt, ldvt )
1054  CALL zungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1055  $ work( nwork ), lwork-nwork+1, ierr )
1056 *
1057 * Generate Q in A
1058 * CWorkspace: need 2*N [tauq, taup] + N [work]
1059 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1060 * RWorkspace: need 0
1061 *
1062  CALL zungbr( 'Q', m, n, n, a, lda, work( itauq ),
1063  $ work( nwork ), lwork-nwork+1, ierr )
1064 *
1065  IF( lwork .GE. m*n + 3*n ) THEN
1066 *
1067 * WORK( IU ) is M by N
1068 *
1069  ldwrku = m
1070  ELSE
1071 *
1072 * WORK(IU) is LDWRKU by N
1073 *
1074  ldwrku = ( lwork - 3*n ) / n
1075  END IF
1076  nwork = iu + ldwrku*n
1077 *
1078 * Perform bidiagonal SVD, computing left singular vectors
1079 * of bidiagonal matrix in RWORK(IRU) and computing right
1080 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1081 * CWorkspace: need 0
1082 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1083 *
1084  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1085  $ n, rwork( irvt ), n, dum, idum,
1086  $ rwork( nrwork ), iwork, info )
1087 *
1088 * Multiply real matrix RWORK(IRVT) by P**H in VT,
1089 * storing the result in WORK(IU), copying to VT
1090 * CWorkspace: need 2*N [tauq, taup] + N*N [U]
1091 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1092 *
1093  CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt,
1094  $ work( iu ), ldwrku, rwork( nrwork ) )
1095  CALL zlacpy( 'F', n, n, work( iu ), ldwrku, vt, ldvt )
1096 *
1097 * Multiply Q in A by real matrix RWORK(IRU), storing the
1098 * result in WORK(IU), copying to A
1099 * CWorkspace: need 2*N [tauq, taup] + N*N [U]
1100 * CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
1101 * RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork]
1102 * RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1103 *
1104  nrwork = irvt
1105  DO 20 i = 1, m, ldwrku
1106  chunk = min( m-i+1, ldwrku )
1107  CALL zlacrm( chunk, n, a( i, 1 ), lda, rwork( iru ),
1108  $ n, work( iu ), ldwrku, rwork( nrwork ) )
1109  CALL zlacpy( 'F', chunk, n, work( iu ), ldwrku,
1110  $ a( i, 1 ), lda )
1111  20 CONTINUE
1112 *
1113  ELSE IF( wntqs ) THEN
1114 *
1115 * Path 5s (M >> N, JOBZ='S')
1116 * Copy A to VT, generate P**H
1117 * CWorkspace: need 2*N [tauq, taup] + N [work]
1118 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1119 * RWorkspace: need 0
1120 *
1121  CALL zlacpy( 'U', n, n, a, lda, vt, ldvt )
1122  CALL zungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1123  $ work( nwork ), lwork-nwork+1, ierr )
1124 *
1125 * Copy A to U, generate Q
1126 * CWorkspace: need 2*N [tauq, taup] + N [work]
1127 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1128 * RWorkspace: need 0
1129 *
1130  CALL zlacpy( 'L', m, n, a, lda, u, ldu )
1131  CALL zungbr( 'Q', m, n, n, u, ldu, work( itauq ),
1132  $ work( nwork ), lwork-nwork+1, ierr )
1133 *
1134 * Perform bidiagonal SVD, computing left singular vectors
1135 * of bidiagonal matrix in RWORK(IRU) and computing right
1136 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1137 * CWorkspace: need 0
1138 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1139 *
1140  iru = nrwork
1141  irvt = iru + n*n
1142  nrwork = irvt + n*n
1143  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1144  $ n, rwork( irvt ), n, dum, idum,
1145  $ rwork( nrwork ), iwork, info )
1146 *
1147 * Multiply real matrix RWORK(IRVT) by P**H in VT,
1148 * storing the result in A, copying to VT
1149 * CWorkspace: need 0
1150 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1151 *
1152  CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1153  $ rwork( nrwork ) )
1154  CALL zlacpy( 'F', n, n, a, lda, vt, ldvt )
1155 *
1156 * Multiply Q in U by real matrix RWORK(IRU), storing the
1157 * result in A, copying to U
1158 * CWorkspace: need 0
1159 * RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1160 *
1161  nrwork = irvt
1162  CALL zlacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1163  $ rwork( nrwork ) )
1164  CALL zlacpy( 'F', m, n, a, lda, u, ldu )
1165  ELSE
1166 *
1167 * Path 5a (M >> N, JOBZ='A')
1168 * Copy A to VT, generate P**H
1169 * CWorkspace: need 2*N [tauq, taup] + N [work]
1170 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1171 * RWorkspace: need 0
1172 *
1173  CALL zlacpy( 'U', n, n, a, lda, vt, ldvt )
1174  CALL zungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1175  $ work( nwork ), lwork-nwork+1, ierr )
1176 *
1177 * Copy A to U, generate Q
1178 * CWorkspace: need 2*N [tauq, taup] + M [work]
1179 * CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
1180 * RWorkspace: need 0
1181 *
1182  CALL zlacpy( 'L', m, n, a, lda, u, ldu )
1183  CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1184  $ work( nwork ), lwork-nwork+1, ierr )
1185 *
1186 * Perform bidiagonal SVD, computing left singular vectors
1187 * of bidiagonal matrix in RWORK(IRU) and computing right
1188 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1189 * CWorkspace: need 0
1190 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1191 *
1192  iru = nrwork
1193  irvt = iru + n*n
1194  nrwork = irvt + n*n
1195  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1196  $ n, rwork( irvt ), n, dum, idum,
1197  $ rwork( nrwork ), iwork, info )
1198 *
1199 * Multiply real matrix RWORK(IRVT) by P**H in VT,
1200 * storing the result in A, copying to VT
1201 * CWorkspace: need 0
1202 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1203 *
1204  CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1205  $ rwork( nrwork ) )
1206  CALL zlacpy( 'F', n, n, a, lda, vt, ldvt )
1207 *
1208 * Multiply Q in U by real matrix RWORK(IRU), storing the
1209 * result in A, copying to U
1210 * CWorkspace: need 0
1211 * RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1212 *
1213  nrwork = irvt
1214  CALL zlacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1215  $ rwork( nrwork ) )
1216  CALL zlacpy( 'F', m, n, a, lda, u, ldu )
1217  END IF
1218 *
1219  ELSE
1220 *
1221 * M .LT. MNTHR2
1222 *
1223 * Path 6 (M >= N, but not much larger)
1224 * Reduce to bidiagonal form without QR decomposition
1225 * Use ZUNMBR to compute singular vectors
1226 *
1227  ie = 1
1228  nrwork = ie + n
1229  itauq = 1
1230  itaup = itauq + n
1231  nwork = itaup + n
1232 *
1233 * Bidiagonalize A
1234 * CWorkspace: need 2*N [tauq, taup] + M [work]
1235 * CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
1236 * RWorkspace: need N [e]
1237 *
1238  CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1239  $ work( itaup ), work( nwork ), lwork-nwork+1,
1240  $ ierr )
1241  IF( wntqn ) THEN
1242 *
1243 * Path 6n (M >= N, JOBZ='N')
1244 * Compute singular values only
1245 * CWorkspace: need 0
1246 * RWorkspace: need N [e] + BDSPAC
1247 *
1248  CALL dbdsdc( 'U', 'N', n, s, rwork( ie ), dum,1,dum,1,
1249  $ dum, idum, rwork( nrwork ), iwork, info )
1250  ELSE IF( wntqo ) THEN
1251  iu = nwork
1252  iru = nrwork
1253  irvt = iru + n*n
1254  nrwork = irvt + n*n
1255  IF( lwork .GE. m*n + 3*n ) THEN
1256 *
1257 * WORK( IU ) is M by N
1258 *
1259  ldwrku = m
1260  ELSE
1261 *
1262 * WORK( IU ) is LDWRKU by N
1263 *
1264  ldwrku = ( lwork - 3*n ) / n
1265  END IF
1266  nwork = iu + ldwrku*n
1267 *
1268 * Path 6o (M >= N, JOBZ='O')
1269 * Perform bidiagonal SVD, computing left singular vectors
1270 * of bidiagonal matrix in RWORK(IRU) and computing right
1271 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1272 * CWorkspace: need 0
1273 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1274 *
1275  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1276  $ n, rwork( irvt ), n, dum, idum,
1277  $ rwork( nrwork ), iwork, info )
1278 *
1279 * Copy real matrix RWORK(IRVT) to complex matrix VT
1280 * Overwrite VT by right singular vectors of A
1281 * CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work]
1282 * CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
1283 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1284 *
1285  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
1286  CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
1287  $ work( itaup ), vt, ldvt, work( nwork ),
1288  $ lwork-nwork+1, ierr )
1289 *
1290  IF( lwork .GE. m*n + 3*n ) THEN
1291 *
1292 * Path 6o-fast
1293 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1294 * Overwrite WORK(IU) by left singular vectors of A, copying
1295 * to A
1296 * CWorkspace: need 2*N [tauq, taup] + M*N [U] + N [work]
1297 * CWorkspace: prefer 2*N [tauq, taup] + M*N [U] + N*NB [work]
1298 * RWorkspace: need N [e] + N*N [RU]
1299 *
1300  CALL zlaset( 'F', m, n, czero, czero, work( iu ),
1301  $ ldwrku )
1302  CALL zlacp2( 'F', n, n, rwork( iru ), n, work( iu ),
1303  $ ldwrku )
1304  CALL zunmbr( 'Q', 'L', 'N', m, n, n, a, lda,
1305  $ work( itauq ), work( iu ), ldwrku,
1306  $ work( nwork ), lwork-nwork+1, ierr )
1307  CALL zlacpy( 'F', m, n, work( iu ), ldwrku, a, lda )
1308  ELSE
1309 *
1310 * Path 6o-slow
1311 * Generate Q in A
1312 * CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work]
1313 * CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
1314 * RWorkspace: need 0
1315 *
1316  CALL zungbr( 'Q', m, n, n, a, lda, work( itauq ),
1317  $ work( nwork ), lwork-nwork+1, ierr )
1318 *
1319 * Multiply Q in A by real matrix RWORK(IRU), storing the
1320 * result in WORK(IU), copying to A
1321 * CWorkspace: need 2*N [tauq, taup] + N*N [U]
1322 * CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
1323 * RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork]
1324 * RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1325 *
1326  nrwork = irvt
1327  DO 30 i = 1, m, ldwrku
1328  chunk = min( m-i+1, ldwrku )
1329  CALL zlacrm( chunk, n, a( i, 1 ), lda,
1330  $ rwork( iru ), n, work( iu ), ldwrku,
1331  $ rwork( nrwork ) )
1332  CALL zlacpy( 'F', chunk, n, work( iu ), ldwrku,
1333  $ a( i, 1 ), lda )
1334  30 CONTINUE
1335  END IF
1336 *
1337  ELSE IF( wntqs ) THEN
1338 *
1339 * Path 6s (M >= N, JOBZ='S')
1340 * Perform bidiagonal SVD, computing left singular vectors
1341 * of bidiagonal matrix in RWORK(IRU) and computing right
1342 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1343 * CWorkspace: need 0
1344 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1345 *
1346  iru = nrwork
1347  irvt = iru + n*n
1348  nrwork = irvt + n*n
1349  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1350  $ n, rwork( irvt ), n, dum, idum,
1351  $ rwork( nrwork ), iwork, info )
1352 *
1353 * Copy real matrix RWORK(IRU) to complex matrix U
1354 * Overwrite U by left singular vectors of A
1355 * CWorkspace: need 2*N [tauq, taup] + N [work]
1356 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1357 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1358 *
1359  CALL zlaset( 'F', m, n, czero, czero, u, ldu )
1360  CALL zlacp2( 'F', n, n, rwork( iru ), n, u, ldu )
1361  CALL zunmbr( 'Q', 'L', 'N', m, n, n, a, lda,
1362  $ work( itauq ), u, ldu, work( nwork ),
1363  $ lwork-nwork+1, ierr )
1364 *
1365 * Copy real matrix RWORK(IRVT) to complex matrix VT
1366 * Overwrite VT by right singular vectors of A
1367 * CWorkspace: need 2*N [tauq, taup] + N [work]
1368 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1369 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1370 *
1371  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
1372  CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
1373  $ work( itaup ), vt, ldvt, work( nwork ),
1374  $ lwork-nwork+1, ierr )
1375  ELSE
1376 *
1377 * Path 6a (M >= N, JOBZ='A')
1378 * Perform bidiagonal SVD, computing left singular vectors
1379 * of bidiagonal matrix in RWORK(IRU) and computing right
1380 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1381 * CWorkspace: need 0
1382 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1383 *
1384  iru = nrwork
1385  irvt = iru + n*n
1386  nrwork = irvt + n*n
1387  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1388  $ n, rwork( irvt ), n, dum, idum,
1389  $ rwork( nrwork ), iwork, info )
1390 *
1391 * Set the right corner of U to identity matrix
1392 *
1393  CALL zlaset( 'F', m, m, czero, czero, u, ldu )
1394  IF( m.GT.n ) THEN
1395  CALL zlaset( 'F', m-n, m-n, czero, cone,
1396  $ u( n+1, n+1 ), ldu )
1397  END IF
1398 *
1399 * Copy real matrix RWORK(IRU) to complex matrix U
1400 * Overwrite U by left singular vectors of A
1401 * CWorkspace: need 2*N [tauq, taup] + M [work]
1402 * CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
1403 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1404 *
1405  CALL zlacp2( 'F', n, n, rwork( iru ), n, u, ldu )
1406  CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
1407  $ work( itauq ), u, ldu, work( nwork ),
1408  $ lwork-nwork+1, ierr )
1409 *
1410 * Copy real matrix RWORK(IRVT) to complex matrix VT
1411 * Overwrite VT by right singular vectors of A
1412 * CWorkspace: need 2*N [tauq, taup] + N [work]
1413 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1414 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1415 *
1416  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
1417  CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
1418  $ work( itaup ), vt, ldvt, work( nwork ),
1419  $ lwork-nwork+1, ierr )
1420  END IF
1421 *
1422  END IF
1423 *
1424  ELSE
1425 *
1426 * A has more columns than rows. If A has sufficiently more
1427 * columns than rows, first reduce using the LQ decomposition (if
1428 * sufficient workspace available)
1429 *
1430  IF( n.GE.mnthr1 ) THEN
1431 *
1432  IF( wntqn ) THEN
1433 *
1434 * Path 1t (N >> M, JOBZ='N')
1435 * No singular vectors to be computed
1436 *
1437  itau = 1
1438  nwork = itau + m
1439 *
1440 * Compute A=L*Q
1441 * CWorkspace: need M [tau] + M [work]
1442 * CWorkspace: prefer M [tau] + M*NB [work]
1443 * RWorkspace: need 0
1444 *
1445  CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1446  $ lwork-nwork+1, ierr )
1447 *
1448 * Zero out above L
1449 *
1450  CALL zlaset( 'U', m-1, m-1, czero, czero, a( 1, 2 ),
1451  $ lda )
1452  ie = 1
1453  itauq = 1
1454  itaup = itauq + m
1455  nwork = itaup + m
1456 *
1457 * Bidiagonalize L in A
1458 * CWorkspace: need 2*M [tauq, taup] + M [work]
1459 * CWorkspace: prefer 2*M [tauq, taup] + 2*M*NB [work]
1460 * RWorkspace: need M [e]
1461 *
1462  CALL zgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1463  $ work( itaup ), work( nwork ), lwork-nwork+1,
1464  $ ierr )
1465  nrwork = ie + m
1466 *
1467 * Perform bidiagonal SVD, compute singular values only
1468 * CWorkspace: need 0
1469 * RWorkspace: need M [e] + BDSPAC
1470 *
1471  CALL dbdsdc( 'U', 'N', m, s, rwork( ie ), dum,1,dum,1,
1472  $ dum, idum, rwork( nrwork ), iwork, info )
1473 *
1474  ELSE IF( wntqo ) THEN
1475 *
1476 * Path 2t (N >> M, JOBZ='O')
1477 * M right singular vectors to be overwritten on A and
1478 * M left singular vectors to be computed in U
1479 *
1480  ivt = 1
1481  ldwkvt = m
1482 *
1483 * WORK(IVT) is M by M
1484 *
1485  il = ivt + ldwkvt*m
1486  IF( lwork .GE. m*n + m*m + 3*m ) THEN
1487 *
1488 * WORK(IL) M by N
1489 *
1490  ldwrkl = m
1491  chunk = n
1492  ELSE
1493 *
1494 * WORK(IL) is M by CHUNK
1495 *
1496  ldwrkl = m
1497  chunk = ( lwork - m*m - 3*m ) / m
1498  END IF
1499  itau = il + ldwrkl*chunk
1500  nwork = itau + m
1501 *
1502 * Compute A=L*Q
1503 * CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1504 * CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1505 * RWorkspace: need 0
1506 *
1507  CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1508  $ lwork-nwork+1, ierr )
1509 *
1510 * Copy L to WORK(IL), zeroing about above it
1511 *
1512  CALL zlacpy( 'L', m, m, a, lda, work( il ), ldwrkl )
1513  CALL zlaset( 'U', m-1, m-1, czero, czero,
1514  $ work( il+ldwrkl ), ldwrkl )
1515 *
1516 * Generate Q in A
1517 * CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1518 * CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1519 * RWorkspace: need 0
1520 *
1521  CALL zunglq( m, n, m, a, lda, work( itau ),
1522  $ work( nwork ), lwork-nwork+1, ierr )
1523  ie = 1
1524  itauq = itau
1525  itaup = itauq + m
1526  nwork = itaup + m
1527 *
1528 * Bidiagonalize L in WORK(IL)
1529 * CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1530 * CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
1531 * RWorkspace: need M [e]
1532 *
1533  CALL zgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1534  $ work( itauq ), work( itaup ), work( nwork ),
1535  $ lwork-nwork+1, ierr )
1536 *
1537 * Perform bidiagonal SVD, computing left singular vectors
1538 * of bidiagonal matrix in RWORK(IRU) and computing right
1539 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1540 * CWorkspace: need 0
1541 * RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1542 *
1543  iru = ie + m
1544  irvt = iru + m*m
1545  nrwork = irvt + m*m
1546  CALL dbdsdc( 'U', 'I', m, s, rwork( ie ), rwork( iru ),
1547  $ m, rwork( irvt ), m, dum, idum,
1548  $ rwork( nrwork ), iwork, info )
1549 *
1550 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1551 * Overwrite WORK(IU) by the left singular vectors of L
1552 * CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1553 * CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
1554 * RWorkspace: need 0
1555 *
1556  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
1557  CALL zunmbr( 'Q', 'L', 'N', m, m, m, work( il ), ldwrkl,
1558  $ work( itauq ), u, ldu, work( nwork ),
1559  $ lwork-nwork+1, ierr )
1560 *
1561 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1562 * Overwrite WORK(IVT) by the right singular vectors of L
1563 * CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1564 * CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
1565 * RWorkspace: need 0
1566 *
1567  CALL zlacp2( 'F', m, m, rwork( irvt ), m, work( ivt ),
1568  $ ldwkvt )
1569  CALL zunmbr( 'P', 'R', 'C', m, m, m, work( il ), ldwrkl,
1570  $ work( itaup ), work( ivt ), ldwkvt,
1571  $ work( nwork ), lwork-nwork+1, ierr )
1572 *
1573 * Multiply right singular vectors of L in WORK(IL) by Q
1574 * in A, storing result in WORK(IL) and copying to A
1575 * CWorkspace: need M*M [VT] + M*M [L]
1576 * CWorkspace: prefer M*M [VT] + M*N [L]
1577 * RWorkspace: need 0
1578 *
1579  DO 40 i = 1, n, chunk
1580  blk = min( n-i+1, chunk )
1581  CALL zgemm( 'N', 'N', m, blk, m, cone, work( ivt ), m,
1582  $ a( 1, i ), lda, czero, work( il ),
1583  $ ldwrkl )
1584  CALL zlacpy( 'F', m, blk, work( il ), ldwrkl,
1585  $ a( 1, i ), lda )
1586  40 CONTINUE
1587 *
1588  ELSE IF( wntqs ) THEN
1589 *
1590 * Path 3t (N >> M, JOBZ='S')
1591 * M right singular vectors to be computed in VT and
1592 * M left singular vectors to be computed in U
1593 *
1594  il = 1
1595 *
1596 * WORK(IL) is M by M
1597 *
1598  ldwrkl = m
1599  itau = il + ldwrkl*m
1600  nwork = itau + m
1601 *
1602 * Compute A=L*Q
1603 * CWorkspace: need M*M [L] + M [tau] + M [work]
1604 * CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
1605 * RWorkspace: need 0
1606 *
1607  CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1608  $ lwork-nwork+1, ierr )
1609 *
1610 * Copy L to WORK(IL), zeroing out above it
1611 *
1612  CALL zlacpy( 'L', m, m, a, lda, work( il ), ldwrkl )
1613  CALL zlaset( 'U', m-1, m-1, czero, czero,
1614  $ work( il+ldwrkl ), ldwrkl )
1615 *
1616 * Generate Q in A
1617 * CWorkspace: need M*M [L] + M [tau] + M [work]
1618 * CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
1619 * RWorkspace: need 0
1620 *
1621  CALL zunglq( m, n, m, a, lda, work( itau ),
1622  $ work( nwork ), lwork-nwork+1, ierr )
1623  ie = 1
1624  itauq = itau
1625  itaup = itauq + m
1626  nwork = itaup + m
1627 *
1628 * Bidiagonalize L in WORK(IL)
1629 * CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1630 * CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
1631 * RWorkspace: need M [e]
1632 *
1633  CALL zgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1634  $ work( itauq ), work( itaup ), work( nwork ),
1635  $ lwork-nwork+1, ierr )
1636 *
1637 * Perform bidiagonal SVD, computing left singular vectors
1638 * of bidiagonal matrix in RWORK(IRU) and computing right
1639 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1640 * CWorkspace: need 0
1641 * RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1642 *
1643  iru = ie + m
1644  irvt = iru + m*m
1645  nrwork = irvt + m*m
1646  CALL dbdsdc( 'U', 'I', m, s, rwork( ie ), rwork( iru ),
1647  $ m, rwork( irvt ), m, dum, idum,
1648  $ rwork( nrwork ), iwork, info )
1649 *
1650 * Copy real matrix RWORK(IRU) to complex matrix U
1651 * Overwrite U by left singular vectors of L
1652 * CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1653 * CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
1654 * RWorkspace: need 0
1655 *
1656  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
1657  CALL zunmbr( 'Q', 'L', 'N', m, m, m, work( il ), ldwrkl,
1658  $ work( itauq ), u, ldu, work( nwork ),
1659  $ lwork-nwork+1, ierr )
1660 *
1661 * Copy real matrix RWORK(IRVT) to complex matrix VT
1662 * Overwrite VT by left singular vectors of L
1663 * CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1664 * CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
1665 * RWorkspace: need 0
1666 *
1667  CALL zlacp2( 'F', m, m, rwork( irvt ), m, vt, ldvt )
1668  CALL zunmbr( 'P', 'R', 'C', m, m, m, work( il ), ldwrkl,
1669  $ work( itaup ), vt, ldvt, work( nwork ),
1670  $ lwork-nwork+1, ierr )
1671 *
1672 * Copy VT to WORK(IL), multiply right singular vectors of L
1673 * in WORK(IL) by Q in A, storing result in VT
1674 * CWorkspace: need M*M [L]
1675 * RWorkspace: need 0
1676 *
1677  CALL zlacpy( 'F', m, m, vt, ldvt, work( il ), ldwrkl )
1678  CALL zgemm( 'N', 'N', m, n, m, cone, work( il ), ldwrkl,
1679  $ a, lda, czero, vt, ldvt )
1680 *
1681  ELSE IF( wntqa ) THEN
1682 *
1683 * Path 4t (N >> M, JOBZ='A')
1684 * N right singular vectors to be computed in VT and
1685 * M left singular vectors to be computed in U
1686 *
1687  ivt = 1
1688 *
1689 * WORK(IVT) is M by M
1690 *
1691  ldwkvt = m
1692  itau = ivt + ldwkvt*m
1693  nwork = itau + m
1694 *
1695 * Compute A=L*Q, copying result to VT
1696 * CWorkspace: need M*M [VT] + M [tau] + M [work]
1697 * CWorkspace: prefer M*M [VT] + M [tau] + M*NB [work]
1698 * RWorkspace: need 0
1699 *
1700  CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1701  $ lwork-nwork+1, ierr )
1702  CALL zlacpy( 'U', m, n, a, lda, vt, ldvt )
1703 *
1704 * Generate Q in VT
1705 * CWorkspace: need M*M [VT] + M [tau] + N [work]
1706 * CWorkspace: prefer M*M [VT] + M [tau] + N*NB [work]
1707 * RWorkspace: need 0
1708 *
1709  CALL zunglq( n, n, m, vt, ldvt, work( itau ),
1710  $ work( nwork ), lwork-nwork+1, ierr )
1711 *
1712 * Produce L in A, zeroing out above it
1713 *
1714  CALL zlaset( 'U', m-1, m-1, czero, czero, a( 1, 2 ),
1715  $ lda )
1716  ie = 1
1717  itauq = itau
1718  itaup = itauq + m
1719  nwork = itaup + m
1720 *
1721 * Bidiagonalize L in A
1722 * CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1723 * CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + 2*M*NB [work]
1724 * RWorkspace: need M [e]
1725 *
1726  CALL zgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1727  $ work( itaup ), work( nwork ), lwork-nwork+1,
1728  $ ierr )
1729 *
1730 * Perform bidiagonal SVD, computing left singular vectors
1731 * of bidiagonal matrix in RWORK(IRU) and computing right
1732 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1733 * CWorkspace: need 0
1734 * RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1735 *
1736  iru = ie + m
1737  irvt = iru + m*m
1738  nrwork = irvt + m*m
1739  CALL dbdsdc( 'U', 'I', m, s, rwork( ie ), rwork( iru ),
1740  $ m, rwork( irvt ), m, dum, idum,
1741  $ rwork( nrwork ), iwork, info )
1742 *
1743 * Copy real matrix RWORK(IRU) to complex matrix U
1744 * Overwrite U by left singular vectors of L
1745 * CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1746 * CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
1747 * RWorkspace: need 0
1748 *
1749  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
1750  CALL zunmbr( 'Q', 'L', 'N', m, m, m, a, lda,
1751  $ work( itauq ), u, ldu, work( nwork ),
1752  $ lwork-nwork+1, ierr )
1753 *
1754 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1755 * Overwrite WORK(IVT) by right singular vectors of L
1756 * CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1757 * CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
1758 * RWorkspace: need 0
1759 *
1760  CALL zlacp2( 'F', m, m, rwork( irvt ), m, work( ivt ),
1761  $ ldwkvt )
1762  CALL zunmbr( 'P', 'R', 'C', m, m, m, a, lda,
1763  $ work( itaup ), work( ivt ), ldwkvt,
1764  $ work( nwork ), lwork-nwork+1, ierr )
1765 *
1766 * Multiply right singular vectors of L in WORK(IVT) by
1767 * Q in VT, storing result in A
1768 * CWorkspace: need M*M [VT]
1769 * RWorkspace: need 0
1770 *
1771  CALL zgemm( 'N', 'N', m, n, m, cone, work( ivt ), ldwkvt,
1772  $ vt, ldvt, czero, a, lda )
1773 *
1774 * Copy right singular vectors of A from A to VT
1775 *
1776  CALL zlacpy( 'F', m, n, a, lda, vt, ldvt )
1777 *
1778  END IF
1779 *
1780  ELSE IF( n.GE.mnthr2 ) THEN
1781 *
1782 * MNTHR2 <= N < MNTHR1
1783 *
1784 * Path 5t (N >> M, but not as much as MNTHR1)
1785 * Reduce to bidiagonal form without QR decomposition, use
1786 * ZUNGBR and matrix multiplication to compute singular vectors
1787 *
1788  ie = 1
1789  nrwork = ie + m
1790  itauq = 1
1791  itaup = itauq + m
1792  nwork = itaup + m
1793 *
1794 * Bidiagonalize A
1795 * CWorkspace: need 2*M [tauq, taup] + N [work]
1796 * CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
1797 * RWorkspace: need M [e]
1798 *
1799  CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1800  $ work( itaup ), work( nwork ), lwork-nwork+1,
1801  $ ierr )
1802 *
1803  IF( wntqn ) THEN
1804 *
1805 * Path 5tn (N >> M, JOBZ='N')
1806 * Compute singular values only
1807 * CWorkspace: need 0
1808 * RWorkspace: need M [e] + BDSPAC
1809 *
1810  CALL dbdsdc( 'L', 'N', m, s, rwork( ie ), dum,1,dum,1,
1811  $ dum, idum, rwork( nrwork ), iwork, info )
1812  ELSE IF( wntqo ) THEN
1813  irvt = nrwork
1814  iru = irvt + m*m
1815  nrwork = iru + m*m
1816  ivt = nwork
1817 *
1818 * Path 5to (N >> M, JOBZ='O')
1819 * Copy A to U, generate Q
1820 * CWorkspace: need 2*M [tauq, taup] + M [work]
1821 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1822 * RWorkspace: need 0
1823 *
1824  CALL zlacpy( 'L', m, m, a, lda, u, ldu )
1825  CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1826  $ work( nwork ), lwork-nwork+1, ierr )
1827 *
1828 * Generate P**H in A
1829 * CWorkspace: need 2*M [tauq, taup] + M [work]
1830 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1831 * RWorkspace: need 0
1832 *
1833  CALL zungbr( 'P', m, n, m, a, lda, work( itaup ),
1834  $ work( nwork ), lwork-nwork+1, ierr )
1835 *
1836  ldwkvt = m
1837  IF( lwork .GE. m*n + 3*m ) THEN
1838 *
1839 * WORK( IVT ) is M by N
1840 *
1841  nwork = ivt + ldwkvt*n
1842  chunk = n
1843  ELSE
1844 *
1845 * WORK( IVT ) is M by CHUNK
1846 *
1847  chunk = ( lwork - 3*m ) / m
1848  nwork = ivt + ldwkvt*chunk
1849  END IF
1850 *
1851 * Perform bidiagonal SVD, computing left singular vectors
1852 * of bidiagonal matrix in RWORK(IRU) and computing right
1853 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1854 * CWorkspace: need 0
1855 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1856 *
1857  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
1858  $ m, rwork( irvt ), m, dum, idum,
1859  $ rwork( nrwork ), iwork, info )
1860 *
1861 * Multiply Q in U by real matrix RWORK(IRVT)
1862 * storing the result in WORK(IVT), copying to U
1863 * CWorkspace: need 2*M [tauq, taup] + M*M [VT]
1864 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1865 *
1866  CALL zlacrm( m, m, u, ldu, rwork( iru ), m, work( ivt ),
1867  $ ldwkvt, rwork( nrwork ) )
1868  CALL zlacpy( 'F', m, m, work( ivt ), ldwkvt, u, ldu )
1869 *
1870 * Multiply RWORK(IRVT) by P**H in A, storing the
1871 * result in WORK(IVT), copying to A
1872 * CWorkspace: need 2*M [tauq, taup] + M*M [VT]
1873 * CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
1874 * RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork]
1875 * RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1876 *
1877  nrwork = iru
1878  DO 50 i = 1, n, chunk
1879  blk = min( n-i+1, chunk )
1880  CALL zlarcm( m, blk, rwork( irvt ), m, a( 1, i ), lda,
1881  $ work( ivt ), ldwkvt, rwork( nrwork ) )
1882  CALL zlacpy( 'F', m, blk, work( ivt ), ldwkvt,
1883  $ a( 1, i ), lda )
1884  50 CONTINUE
1885  ELSE IF( wntqs ) THEN
1886 *
1887 * Path 5ts (N >> M, JOBZ='S')
1888 * Copy A to U, generate Q
1889 * CWorkspace: need 2*M [tauq, taup] + M [work]
1890 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1891 * RWorkspace: need 0
1892 *
1893  CALL zlacpy( 'L', m, m, a, lda, u, ldu )
1894  CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1895  $ work( nwork ), lwork-nwork+1, ierr )
1896 *
1897 * Copy A to VT, generate P**H
1898 * CWorkspace: need 2*M [tauq, taup] + M [work]
1899 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1900 * RWorkspace: need 0
1901 *
1902  CALL zlacpy( 'U', m, n, a, lda, vt, ldvt )
1903  CALL zungbr( 'P', m, n, m, vt, ldvt, work( itaup ),
1904  $ work( nwork ), lwork-nwork+1, ierr )
1905 *
1906 * Perform bidiagonal SVD, computing left singular vectors
1907 * of bidiagonal matrix in RWORK(IRU) and computing right
1908 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1909 * CWorkspace: need 0
1910 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1911 *
1912  irvt = nrwork
1913  iru = irvt + m*m
1914  nrwork = iru + m*m
1915  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
1916  $ m, rwork( irvt ), m, dum, idum,
1917  $ rwork( nrwork ), iwork, info )
1918 *
1919 * Multiply Q in U by real matrix RWORK(IRU), storing the
1920 * result in A, copying to U
1921 * CWorkspace: need 0
1922 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1923 *
1924  CALL zlacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1925  $ rwork( nrwork ) )
1926  CALL zlacpy( 'F', m, m, a, lda, u, ldu )
1927 *
1928 * Multiply real matrix RWORK(IRVT) by P**H in VT,
1929 * storing the result in A, copying to VT
1930 * CWorkspace: need 0
1931 * RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1932 *
1933  nrwork = iru
1934  CALL zlarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1935  $ rwork( nrwork ) )
1936  CALL zlacpy( 'F', m, n, a, lda, vt, ldvt )
1937  ELSE
1938 *
1939 * Path 5ta (N >> M, JOBZ='A')
1940 * Copy A to U, generate Q
1941 * CWorkspace: need 2*M [tauq, taup] + M [work]
1942 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1943 * RWorkspace: need 0
1944 *
1945  CALL zlacpy( 'L', m, m, a, lda, u, ldu )
1946  CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1947  $ work( nwork ), lwork-nwork+1, ierr )
1948 *
1949 * Copy A to VT, generate P**H
1950 * CWorkspace: need 2*M [tauq, taup] + N [work]
1951 * CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
1952 * RWorkspace: need 0
1953 *
1954  CALL zlacpy( 'U', m, n, a, lda, vt, ldvt )
1955  CALL zungbr( 'P', n, n, m, vt, ldvt, work( itaup ),
1956  $ work( nwork ), lwork-nwork+1, ierr )
1957 *
1958 * Perform bidiagonal SVD, computing left singular vectors
1959 * of bidiagonal matrix in RWORK(IRU) and computing right
1960 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1961 * CWorkspace: need 0
1962 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1963 *
1964  irvt = nrwork
1965  iru = irvt + m*m
1966  nrwork = iru + m*m
1967  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
1968  $ m, rwork( irvt ), m, dum, idum,
1969  $ rwork( nrwork ), iwork, info )
1970 *
1971 * Multiply Q in U by real matrix RWORK(IRU), storing the
1972 * result in A, copying to U
1973 * CWorkspace: need 0
1974 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1975 *
1976  CALL zlacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1977  $ rwork( nrwork ) )
1978  CALL zlacpy( 'F', m, m, a, lda, u, ldu )
1979 *
1980 * Multiply real matrix RWORK(IRVT) by P**H in VT,
1981 * storing the result in A, copying to VT
1982 * CWorkspace: need 0
1983 * RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1984 *
1985  nrwork = iru
1986  CALL zlarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1987  $ rwork( nrwork ) )
1988  CALL zlacpy( 'F', m, n, a, lda, vt, ldvt )
1989  END IF
1990 *
1991  ELSE
1992 *
1993 * N .LT. MNTHR2
1994 *
1995 * Path 6t (N > M, but not much larger)
1996 * Reduce to bidiagonal form without LQ decomposition
1997 * Use ZUNMBR to compute singular vectors
1998 *
1999  ie = 1
2000  nrwork = ie + m
2001  itauq = 1
2002  itaup = itauq + m
2003  nwork = itaup + m
2004 *
2005 * Bidiagonalize A
2006 * CWorkspace: need 2*M [tauq, taup] + N [work]
2007 * CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
2008 * RWorkspace: need M [e]
2009 *
2010  CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2011  $ work( itaup ), work( nwork ), lwork-nwork+1,
2012  $ ierr )
2013  IF( wntqn ) THEN
2014 *
2015 * Path 6tn (N > M, JOBZ='N')
2016 * Compute singular values only
2017 * CWorkspace: need 0
2018 * RWorkspace: need M [e] + BDSPAC
2019 *
2020  CALL dbdsdc( 'L', 'N', m, s, rwork( ie ), dum,1,dum,1,
2021  $ dum, idum, rwork( nrwork ), iwork, info )
2022  ELSE IF( wntqo ) THEN
2023 * Path 6to (N > M, JOBZ='O')
2024  ldwkvt = m
2025  ivt = nwork
2026  IF( lwork .GE. m*n + 3*m ) THEN
2027 *
2028 * WORK( IVT ) is M by N
2029 *
2030  CALL zlaset( 'F', m, n, czero, czero, work( ivt ),
2031  $ ldwkvt )
2032  nwork = ivt + ldwkvt*n
2033  ELSE
2034 *
2035 * WORK( IVT ) is M by CHUNK
2036 *
2037  chunk = ( lwork - 3*m ) / m
2038  nwork = ivt + ldwkvt*chunk
2039  END IF
2040 *
2041 * Perform bidiagonal SVD, computing left singular vectors
2042 * of bidiagonal matrix in RWORK(IRU) and computing right
2043 * singular vectors of bidiagonal matrix in RWORK(IRVT)
2044 * CWorkspace: need 0
2045 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2046 *
2047  irvt = nrwork
2048  iru = irvt + m*m
2049  nrwork = iru + m*m
2050  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
2051  $ m, rwork( irvt ), m, dum, idum,
2052  $ rwork( nrwork ), iwork, info )
2053 *
2054 * Copy real matrix RWORK(IRU) to complex matrix U
2055 * Overwrite U by left singular vectors of A
2056 * CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work]
2057 * CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
2058 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2059 *
2060  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
2061  CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
2062  $ work( itauq ), u, ldu, work( nwork ),
2063  $ lwork-nwork+1, ierr )
2064 *
2065  IF( lwork .GE. m*n + 3*m ) THEN
2066 *
2067 * Path 6to-fast
2068 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
2069 * Overwrite WORK(IVT) by right singular vectors of A,
2070 * copying to A
2071 * CWorkspace: need 2*M [tauq, taup] + M*N [VT] + M [work]
2072 * CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] + M*NB [work]
2073 * RWorkspace: need M [e] + M*M [RVT]
2074 *
2075  CALL zlacp2( 'F', m, m, rwork( irvt ), m, work( ivt ),
2076  $ ldwkvt )
2077  CALL zunmbr( 'P', 'R', 'C', m, n, m, a, lda,
2078  $ work( itaup ), work( ivt ), ldwkvt,
2079  $ work( nwork ), lwork-nwork+1, ierr )
2080  CALL zlacpy( 'F', m, n, work( ivt ), ldwkvt, a, lda )
2081  ELSE
2082 *
2083 * Path 6to-slow
2084 * Generate P**H in A
2085 * CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work]
2086 * CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
2087 * RWorkspace: need 0
2088 *
2089  CALL zungbr( 'P', m, n, m, a, lda, work( itaup ),
2090  $ work( nwork ), lwork-nwork+1, ierr )
2091 *
2092 * Multiply Q in A by real matrix RWORK(IRU), storing the
2093 * result in WORK(IU), copying to A
2094 * CWorkspace: need 2*M [tauq, taup] + M*M [VT]
2095 * CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
2096 * RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork]
2097 * RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
2098 *
2099  nrwork = iru
2100  DO 60 i = 1, n, chunk
2101  blk = min( n-i+1, chunk )
2102  CALL zlarcm( m, blk, rwork( irvt ), m, a( 1, i ),
2103  $ lda, work( ivt ), ldwkvt,
2104  $ rwork( nrwork ) )
2105  CALL zlacpy( 'F', m, blk, work( ivt ), ldwkvt,
2106  $ a( 1, i ), lda )
2107  60 CONTINUE
2108  END IF
2109  ELSE IF( wntqs ) THEN
2110 *
2111 * Path 6ts (N > M, JOBZ='S')
2112 * Perform bidiagonal SVD, computing left singular vectors
2113 * of bidiagonal matrix in RWORK(IRU) and computing right
2114 * singular vectors of bidiagonal matrix in RWORK(IRVT)
2115 * CWorkspace: need 0
2116 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2117 *
2118  irvt = nrwork
2119  iru = irvt + m*m
2120  nrwork = iru + m*m
2121  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
2122  $ m, rwork( irvt ), m, dum, idum,
2123  $ rwork( nrwork ), iwork, info )
2124 *
2125 * Copy real matrix RWORK(IRU) to complex matrix U
2126 * Overwrite U by left singular vectors of A
2127 * CWorkspace: need 2*M [tauq, taup] + M [work]
2128 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2129 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2130 *
2131  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
2132  CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
2133  $ work( itauq ), u, ldu, work( nwork ),
2134  $ lwork-nwork+1, ierr )
2135 *
2136 * Copy real matrix RWORK(IRVT) to complex matrix VT
2137 * Overwrite VT by right singular vectors of A
2138 * CWorkspace: need 2*M [tauq, taup] + M [work]
2139 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2140 * RWorkspace: need M [e] + M*M [RVT]
2141 *
2142  CALL zlaset( 'F', m, n, czero, czero, vt, ldvt )
2143  CALL zlacp2( 'F', m, m, rwork( irvt ), m, vt, ldvt )
2144  CALL zunmbr( 'P', 'R', 'C', m, n, m, a, lda,
2145  $ work( itaup ), vt, ldvt, work( nwork ),
2146  $ lwork-nwork+1, ierr )
2147  ELSE
2148 *
2149 * Path 6ta (N > M, JOBZ='A')
2150 * Perform bidiagonal SVD, computing left singular vectors
2151 * of bidiagonal matrix in RWORK(IRU) and computing right
2152 * singular vectors of bidiagonal matrix in RWORK(IRVT)
2153 * CWorkspace: need 0
2154 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2155 *
2156  irvt = nrwork
2157  iru = irvt + m*m
2158  nrwork = iru + m*m
2159 *
2160  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
2161  $ m, rwork( irvt ), m, dum, idum,
2162  $ rwork( nrwork ), iwork, info )
2163 *
2164 * Copy real matrix RWORK(IRU) to complex matrix U
2165 * Overwrite U by left singular vectors of A
2166 * CWorkspace: need 2*M [tauq, taup] + M [work]
2167 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2168 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2169 *
2170  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
2171  CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
2172  $ work( itauq ), u, ldu, work( nwork ),
2173  $ lwork-nwork+1, ierr )
2174 *
2175 * Set all of VT to identity matrix
2176 *
2177  CALL zlaset( 'F', n, n, czero, cone, vt, ldvt )
2178 *
2179 * Copy real matrix RWORK(IRVT) to complex matrix VT
2180 * Overwrite VT by right singular vectors of A
2181 * CWorkspace: need 2*M [tauq, taup] + N [work]
2182 * CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
2183 * RWorkspace: need M [e] + M*M [RVT]
2184 *
2185  CALL zlacp2( 'F', m, m, rwork( irvt ), m, vt, ldvt )
2186  CALL zunmbr( 'P', 'R', 'C', n, n, m, a, lda,
2187  $ work( itaup ), vt, ldvt, work( nwork ),
2188  $ lwork-nwork+1, ierr )
2189  END IF
2190 *
2191  END IF
2192 *
2193  END IF
2194 *
2195 * Undo scaling if necessary
2196 *
2197  IF( iscl.EQ.1 ) THEN
2198  IF( anrm.GT.bignum )
2199  $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
2200  $ ierr )
2201  IF( info.NE.0 .AND. anrm.GT.bignum )
2202  $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn-1, 1,
2203  $ rwork( ie ), minmn, ierr )
2204  IF( anrm.LT.smlnum )
2205  $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
2206  $ ierr )
2207  IF( info.NE.0 .AND. anrm.LT.smlnum )
2208  $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn-1, 1,
2209  $ rwork( ie ), minmn, ierr )
2210  END IF
2211 *
2212 * Return optimal workspace in WORK(1)
2213 *
2214  work( 1 ) = maxwrk
2215 *
2216  RETURN
2217 *
2218 * End of ZGESDD
2219 *
2220  END
zlacp2
subroutine zlacp2(UPLO, M, N, A, LDA, B, LDB)
ZLACP2 copies all or part of a real two-dimensional array to a complex array.
Definition: zlacp2.f:106
dbdsdc
subroutine dbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
DBDSDC
Definition: dbdsdc.f:207
zgesdd
subroutine zgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, IWORK, INFO)
ZGESDD
Definition: zgesdd.f:228
zungbr
subroutine zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
Definition: zungbr.f:159
dlascl
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:145
zgelqf
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
Definition: zgelqf.f:145
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:108
zlascl
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:145
zungqr
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zunmbr
subroutine zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
Definition: zunmbr.f:198
zunglq
subroutine zunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGLQ
Definition: zunglq.f:129
zgeqrf
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
zlarcm
subroutine zlarcm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
ZLARCM copies all or part of a real two-dimensional array to a complex array.
Definition: zlarcm.f:116
zgebrd
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD
Definition: zgebrd.f:207
zlacrm
subroutine zlacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
ZLACRM multiplies a complex matrix by a square real matrix.
Definition: zlacrm.f:116