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