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