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