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