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