LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
cgesvd.f
Go to the documentation of this file.
1*> \brief <b> CGESVD computes the singular value decomposition (SVD) for GE matrices</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download CGESVD + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesvd.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesvd.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesvd.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
20* WORK, LWORK, RWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER JOBU, JOBVT
24* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
25* ..
26* .. Array Arguments ..
27* REAL RWORK( * ), S( * )
28* COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
29* $ WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CGESVD computes the singular value decomposition (SVD) of a complex
39*> M-by-N matrix A, optionally computing the left and/or right singular
40*> vectors. The SVD is written
41*>
42*> A = U * SIGMA * conjugate-transpose(V)
43*>
44*> where SIGMA is an M-by-N matrix which is zero except for its
45*> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
46*> V is an N-by-N unitary matrix. The diagonal elements of SIGMA
47*> are the singular values of A; they are real and non-negative, and
48*> are returned in descending order. The first min(m,n) columns of
49*> U and V are the left and right singular vectors of A.
50*>
51*> Note that the routine returns V**H, not V.
52*> \endverbatim
53*
54* Arguments:
55* ==========
56*
57*> \param[in] JOBU
58*> \verbatim
59*> JOBU is CHARACTER*1
60*> Specifies options for computing all or part of the matrix U:
61*> = 'A': all M columns of U are returned in array U:
62*> = 'S': the first min(m,n) columns of U (the left singular
63*> vectors) are returned in the array U;
64*> = 'O': the first min(m,n) columns of U (the left singular
65*> vectors) are overwritten on the array A;
66*> = 'N': no columns of U (no left singular vectors) are
67*> computed.
68*> \endverbatim
69*>
70*> \param[in] JOBVT
71*> \verbatim
72*> JOBVT is CHARACTER*1
73*> Specifies options for computing all or part of the matrix
74*> V**H:
75*> = 'A': all N rows of V**H are returned in the array VT;
76*> = 'S': the first min(m,n) rows of V**H (the right singular
77*> vectors) are returned in the array VT;
78*> = 'O': the first min(m,n) rows of V**H (the right singular
79*> vectors) are overwritten on the array A;
80*> = 'N': no rows of V**H (no right singular vectors) are
81*> computed.
82*>
83*> JOBVT and JOBU cannot both be 'O'.
84*> \endverbatim
85*>
86*> \param[in] M
87*> \verbatim
88*> M is INTEGER
89*> The number of rows of the input matrix A. M >= 0.
90*> \endverbatim
91*>
92*> \param[in] N
93*> \verbatim
94*> N is INTEGER
95*> The number of columns of the input matrix A. N >= 0.
96*> \endverbatim
97*>
98*> \param[in,out] A
99*> \verbatim
100*> A is COMPLEX array, dimension (LDA,N)
101*> On entry, the M-by-N matrix A.
102*> On exit,
103*> if JOBU = 'O', A is overwritten with the first min(m,n)
104*> columns of U (the left singular vectors,
105*> stored columnwise);
106*> if JOBVT = 'O', A is overwritten with the first min(m,n)
107*> rows of V**H (the right singular vectors,
108*> stored rowwise);
109*> if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
110*> are destroyed.
111*> \endverbatim
112*>
113*> \param[in] LDA
114*> \verbatim
115*> LDA is INTEGER
116*> The leading dimension of the array A. LDA >= max(1,M).
117*> \endverbatim
118*>
119*> \param[out] S
120*> \verbatim
121*> S is REAL array, dimension (min(M,N))
122*> The singular values of A, sorted so that S(i) >= S(i+1).
123*> \endverbatim
124*>
125*> \param[out] U
126*> \verbatim
127*> U is COMPLEX array, dimension (LDU,UCOL)
128*> (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
129*> If JOBU = 'A', U contains the M-by-M unitary matrix U;
130*> if JOBU = 'S', U contains the first min(m,n) columns of U
131*> (the left singular vectors, stored columnwise);
132*> if JOBU = 'N' or 'O', U is not referenced.
133*> \endverbatim
134*>
135*> \param[in] LDU
136*> \verbatim
137*> LDU is INTEGER
138*> The leading dimension of the array U. LDU >= 1; if
139*> JOBU = 'S' or 'A', LDU >= M.
140*> \endverbatim
141*>
142*> \param[out] VT
143*> \verbatim
144*> VT is COMPLEX array, dimension (LDVT,N)
145*> If JOBVT = 'A', VT contains the N-by-N unitary matrix
146*> V**H;
147*> if JOBVT = 'S', VT contains the first min(m,n) rows of
148*> V**H (the right singular vectors, stored rowwise);
149*> if JOBVT = 'N' or 'O', VT is not referenced.
150*> \endverbatim
151*>
152*> \param[in] LDVT
153*> \verbatim
154*> LDVT is INTEGER
155*> The leading dimension of the array VT. LDVT >= 1; if
156*> JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
157*> \endverbatim
158*>
159*> \param[out] WORK
160*> \verbatim
161*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
162*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
163*> \endverbatim
164*>
165*> \param[in] LWORK
166*> \verbatim
167*> LWORK is INTEGER
168*> The dimension of the array WORK.
169*> LWORK >= MAX(1,2*MIN(M,N)+MAX(M,N)).
170*> For good performance, LWORK should generally be larger.
171*>
172*> If LWORK = -1, then a workspace query is assumed; the routine
173*> only calculates the optimal size of the WORK array, returns
174*> this value as the first entry of the WORK array, and no error
175*> message related to LWORK is issued by XERBLA.
176*> \endverbatim
177*>
178*> \param[out] RWORK
179*> \verbatim
180*> RWORK is REAL array, dimension (5*min(M,N))
181*> On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the
182*> unconverged superdiagonal elements of an upper bidiagonal
183*> matrix B whose diagonal is in S (not necessarily sorted).
184*> B satisfies A = U * B * VT, so it has the same singular
185*> values as A, and singular vectors related by U and VT.
186*> \endverbatim
187*>
188*> \param[out] INFO
189*> \verbatim
190*> INFO is INTEGER
191*> = 0: successful exit.
192*> < 0: if INFO = -i, the i-th argument had an illegal value.
193*> > 0: if CBDSQR did not converge, INFO specifies how many
194*> superdiagonals of an intermediate bidiagonal form B
195*> did not converge to zero. See the description of RWORK
196*> above for details.
197*> \endverbatim
198*
199* Authors:
200* ========
201*
202*> \author Univ. of Tennessee
203*> \author Univ. of California Berkeley
204*> \author Univ. of Colorado Denver
205*> \author NAG Ltd.
206*
207*> \ingroup gesvd
208*
209* =====================================================================
210 SUBROUTINE cgesvd( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT,
211 $ LDVT,
212 $ WORK, LWORK, RWORK, INFO )
213*
214* -- LAPACK driver routine --
215* -- LAPACK is a software package provided by Univ. of Tennessee, --
216* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217*
218* .. Scalar Arguments ..
219 CHARACTER JOBU, JOBVT
220 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
221* ..
222* .. Array Arguments ..
223 REAL RWORK( * ), S( * )
224 COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
225 $ work( * )
226* ..
227*
228* =====================================================================
229*
230* .. Parameters ..
231 COMPLEX CZERO, CONE
232 PARAMETER ( CZERO = ( 0.0e0, 0.0e0 ),
233 $ cone = ( 1.0e0, 0.0e0 ) )
234 REAL ZERO, ONE
235 PARAMETER ( ZERO = 0.0e0, one = 1.0e0 )
236* ..
237* .. Local Scalars ..
238 LOGICAL LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS,
239 $ WNTVA, WNTVAS, WNTVN, WNTVO, WNTVS
240 INTEGER BLK, CHUNK, I, IE, IERR, IR, IRWORK, ISCL,
241 $ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,
242 $ maxwrk, minmn, minwrk, mnthr, ncu, ncvt, nru,
243 $ nrvt, wrkbl
244 INTEGER LWORK_CGEQRF, LWORK_CUNGQR_N, LWORK_CUNGQR_M,
245 $ LWORK_CGEBRD, LWORK_CUNGBR_P, LWORK_CUNGBR_Q,
246 $ lwork_cgelqf, lwork_cunglq_n, lwork_cunglq_m
247 REAL ANRM, BIGNUM, EPS, SMLNUM
248* ..
249* .. Local Arrays ..
250 REAL DUM( 1 )
251 COMPLEX CDUM( 1 )
252* ..
253* .. External Subroutines ..
254 EXTERNAL cbdsqr, cgebrd, cgelqf, cgemm, cgeqrf,
255 $ clacpy,
256 $ clascl, claset, cungbr, cunglq, cungqr, cunmbr,
257 $ slascl, xerbla
258* ..
259* .. External Functions ..
260 LOGICAL LSAME
261 INTEGER ILAENV
262 REAL CLANGE, SLAMCH, SROUNDUP_LWORK
263 EXTERNAL lsame, ilaenv, clange, slamch,
264 $ sroundup_lwork
265* ..
266* .. Intrinsic Functions ..
267 INTRINSIC max, min, sqrt
268* ..
269* .. Executable Statements ..
270*
271* Test the input arguments
272*
273 info = 0
274 minmn = min( m, n )
275 wntua = lsame( jobu, 'A' )
276 wntus = lsame( jobu, 'S' )
277 wntuas = wntua .OR. wntus
278 wntuo = lsame( jobu, 'O' )
279 wntun = lsame( jobu, 'N' )
280 wntva = lsame( jobvt, 'A' )
281 wntvs = lsame( jobvt, 'S' )
282 wntvas = wntva .OR. wntvs
283 wntvo = lsame( jobvt, 'O' )
284 wntvn = lsame( jobvt, 'N' )
285 lquery = ( lwork.EQ.-1 )
286*
287 IF( .NOT.( wntua .OR. wntus .OR. wntuo .OR. wntun ) ) THEN
288 info = -1
289 ELSE IF( .NOT.( wntva .OR. wntvs .OR. wntvo .OR. wntvn ) .OR.
290 $ ( wntvo .AND. wntuo ) ) THEN
291 info = -2
292 ELSE IF( m.LT.0 ) THEN
293 info = -3
294 ELSE IF( n.LT.0 ) THEN
295 info = -4
296 ELSE IF( lda.LT.max( 1, m ) ) THEN
297 info = -6
298 ELSE IF( ldu.LT.1 .OR. ( wntuas .AND. ldu.LT.m ) ) THEN
299 info = -9
300 ELSE IF( ldvt.LT.1 .OR. ( wntva .AND. ldvt.LT.n ) .OR.
301 $ ( wntvs .AND. ldvt.LT.minmn ) ) THEN
302 info = -11
303 END IF
304*
305* Compute workspace
306* (Note: Comments in the code beginning "Workspace:" describe the
307* minimal amount of workspace needed at that point in the code,
308* as well as the preferred amount for good performance.
309* CWorkspace refers to complex workspace, and RWorkspace to
310* real workspace. NB refers to the optimal block size for the
311* immediately following subroutine, as returned by ILAENV.)
312*
313 IF( info.EQ.0 ) THEN
314 minwrk = 1
315 maxwrk = 1
316 IF( m.GE.n .AND. minmn.GT.0 ) THEN
317*
318* Space needed for ZBDSQR is BDSPAC = 5*N
319*
320 mnthr = ilaenv( 6, 'CGESVD', jobu // jobvt, m, n, 0, 0 )
321* Compute space needed for CGEQRF
322 CALL cgeqrf( m, n, a, lda, cdum(1), cdum(1), -1, ierr )
323 lwork_cgeqrf = int( cdum(1) )
324* Compute space needed for CUNGQR
325 CALL cungqr( m, n, n, a, lda, cdum(1), cdum(1), -1,
326 $ ierr )
327 lwork_cungqr_n = int( cdum(1) )
328 CALL cungqr( m, m, n, a, lda, cdum(1), cdum(1), -1,
329 $ ierr )
330 lwork_cungqr_m = int( cdum(1) )
331* Compute space needed for CGEBRD
332 CALL cgebrd( n, n, a, lda, s, dum(1), cdum(1),
333 $ cdum(1), cdum(1), -1, ierr )
334 lwork_cgebrd = int( cdum(1) )
335* Compute space needed for CUNGBR
336 CALL cungbr( 'P', n, n, n, a, lda, cdum(1),
337 $ cdum(1), -1, ierr )
338 lwork_cungbr_p = int( cdum(1) )
339 CALL cungbr( 'Q', n, n, n, a, lda, cdum(1),
340 $ cdum(1), -1, ierr )
341 lwork_cungbr_q = int( cdum(1) )
342*
343 mnthr = ilaenv( 6, 'CGESVD', jobu // jobvt, m, n, 0, 0 )
344 IF( m.GE.mnthr ) THEN
345 IF( wntun ) THEN
346*
347* Path 1 (M much larger than N, JOBU='N')
348*
349 maxwrk = n + lwork_cgeqrf
350 maxwrk = max( maxwrk, 2*n+lwork_cgebrd )
351 IF( wntvo .OR. wntvas )
352 $ maxwrk = max( maxwrk, 2*n+lwork_cungbr_p )
353 minwrk = 3*n
354 ELSE IF( wntuo .AND. wntvn ) THEN
355*
356* Path 2 (M much larger than N, JOBU='O', JOBVT='N')
357*
358 wrkbl = n + lwork_cgeqrf
359 wrkbl = max( wrkbl, n+lwork_cungqr_n )
360 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
361 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
362 maxwrk = max( n*n+wrkbl, n*n+m*n )
363 minwrk = 2*n + m
364 ELSE IF( wntuo .AND. wntvas ) THEN
365*
366* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
367* 'A')
368*
369 wrkbl = n + lwork_cgeqrf
370 wrkbl = max( wrkbl, n+lwork_cungqr_n )
371 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
372 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
373 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
374 maxwrk = max( n*n+wrkbl, n*n+m*n )
375 minwrk = 2*n + m
376 ELSE IF( wntus .AND. wntvn ) THEN
377*
378* Path 4 (M much larger than N, JOBU='S', JOBVT='N')
379*
380 wrkbl = n + lwork_cgeqrf
381 wrkbl = max( wrkbl, n+lwork_cungqr_n )
382 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
383 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
384 maxwrk = n*n + wrkbl
385 minwrk = 2*n + m
386 ELSE IF( wntus .AND. wntvo ) THEN
387*
388* Path 5 (M much larger than N, JOBU='S', JOBVT='O')
389*
390 wrkbl = n + lwork_cgeqrf
391 wrkbl = max( wrkbl, n+lwork_cungqr_n )
392 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
393 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
394 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
395 maxwrk = 2*n*n + wrkbl
396 minwrk = 2*n + m
397 ELSE IF( wntus .AND. wntvas ) THEN
398*
399* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
400* 'A')
401*
402 wrkbl = n + lwork_cgeqrf
403 wrkbl = max( wrkbl, n+lwork_cungqr_n )
404 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
405 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
406 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
407 maxwrk = n*n + wrkbl
408 minwrk = 2*n + m
409 ELSE IF( wntua .AND. wntvn ) THEN
410*
411* Path 7 (M much larger than N, JOBU='A', JOBVT='N')
412*
413 wrkbl = n + lwork_cgeqrf
414 wrkbl = max( wrkbl, n+lwork_cungqr_m )
415 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
416 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
417 maxwrk = n*n + wrkbl
418 minwrk = 2*n + m
419 ELSE IF( wntua .AND. wntvo ) THEN
420*
421* Path 8 (M much larger than N, JOBU='A', JOBVT='O')
422*
423 wrkbl = n + lwork_cgeqrf
424 wrkbl = max( wrkbl, n+lwork_cungqr_m )
425 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
426 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
427 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
428 maxwrk = 2*n*n + wrkbl
429 minwrk = 2*n + m
430 ELSE IF( wntua .AND. wntvas ) THEN
431*
432* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
433* 'A')
434*
435 wrkbl = n + lwork_cgeqrf
436 wrkbl = max( wrkbl, n+lwork_cungqr_m )
437 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
438 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
439 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
440 maxwrk = n*n + wrkbl
441 minwrk = 2*n + m
442 END IF
443 ELSE
444*
445* Path 10 (M at least N, but not much larger)
446*
447 CALL cgebrd( m, n, a, lda, s, dum(1), cdum(1),
448 $ cdum(1), cdum(1), -1, ierr )
449 lwork_cgebrd = int( cdum(1) )
450 maxwrk = 2*n + lwork_cgebrd
451 IF( wntus .OR. wntuo ) THEN
452 CALL cungbr( 'Q', m, n, n, a, lda, cdum(1),
453 $ cdum(1), -1, ierr )
454 lwork_cungbr_q = int( cdum(1) )
455 maxwrk = max( maxwrk, 2*n+lwork_cungbr_q )
456 END IF
457 IF( wntua ) THEN
458 CALL cungbr( 'Q', m, m, n, a, lda, cdum(1),
459 $ cdum(1), -1, ierr )
460 lwork_cungbr_q = int( cdum(1) )
461 maxwrk = max( maxwrk, 2*n+lwork_cungbr_q )
462 END IF
463 IF( .NOT.wntvn ) THEN
464 maxwrk = max( maxwrk, 2*n+lwork_cungbr_p )
465 END IF
466 minwrk = 2*n + m
467 END IF
468 ELSE IF( minmn.GT.0 ) THEN
469*
470* Space needed for CBDSQR is BDSPAC = 5*M
471*
472 mnthr = ilaenv( 6, 'CGESVD', jobu // jobvt, m, n, 0, 0 )
473* Compute space needed for CGELQF
474 CALL cgelqf( m, n, a, lda, cdum(1), cdum(1), -1, ierr )
475 lwork_cgelqf = int( cdum(1) )
476* Compute space needed for CUNGLQ
477 CALL cunglq( n, n, m, cdum(1), n, cdum(1), cdum(1), -1,
478 $ ierr )
479 lwork_cunglq_n = int( cdum(1) )
480 CALL cunglq( m, n, m, a, lda, cdum(1), cdum(1), -1,
481 $ ierr )
482 lwork_cunglq_m = int( cdum(1) )
483* Compute space needed for CGEBRD
484 CALL cgebrd( m, m, a, lda, s, dum(1), cdum(1),
485 $ cdum(1), cdum(1), -1, ierr )
486 lwork_cgebrd = int( cdum(1) )
487* Compute space needed for CUNGBR P
488 CALL cungbr( 'P', m, m, m, a, n, cdum(1),
489 $ cdum(1), -1, ierr )
490 lwork_cungbr_p = int( cdum(1) )
491* Compute space needed for CUNGBR Q
492 CALL cungbr( 'Q', m, m, m, a, n, cdum(1),
493 $ cdum(1), -1, ierr )
494 lwork_cungbr_q = int( cdum(1) )
495 IF( n.GE.mnthr ) THEN
496 IF( wntvn ) THEN
497*
498* Path 1t(N much larger than M, JOBVT='N')
499*
500 maxwrk = m + lwork_cgelqf
501 maxwrk = max( maxwrk, 2*m+lwork_cgebrd )
502 IF( wntuo .OR. wntuas )
503 $ maxwrk = max( maxwrk, 2*m+lwork_cungbr_q )
504 minwrk = 3*m
505 ELSE IF( wntvo .AND. wntun ) THEN
506*
507* Path 2t(N much larger than M, JOBU='N', JOBVT='O')
508*
509 wrkbl = m + lwork_cgelqf
510 wrkbl = max( wrkbl, m+lwork_cunglq_m )
511 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
512 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
513 maxwrk = max( m*m+wrkbl, m*m+m*n )
514 minwrk = 2*m + n
515 ELSE IF( wntvo .AND. wntuas ) THEN
516*
517* Path 3t(N much larger than M, JOBU='S' or 'A',
518* JOBVT='O')
519*
520 wrkbl = m + lwork_cgelqf
521 wrkbl = max( wrkbl, m+lwork_cunglq_m )
522 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
523 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
524 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
525 maxwrk = max( m*m+wrkbl, m*m+m*n )
526 minwrk = 2*m + n
527 ELSE IF( wntvs .AND. wntun ) THEN
528*
529* Path 4t(N much larger than M, JOBU='N', JOBVT='S')
530*
531 wrkbl = m + lwork_cgelqf
532 wrkbl = max( wrkbl, m+lwork_cunglq_m )
533 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
534 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
535 maxwrk = m*m + wrkbl
536 minwrk = 2*m + n
537 ELSE IF( wntvs .AND. wntuo ) THEN
538*
539* Path 5t(N much larger than M, JOBU='O', JOBVT='S')
540*
541 wrkbl = m + lwork_cgelqf
542 wrkbl = max( wrkbl, m+lwork_cunglq_m )
543 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
544 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
545 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
546 maxwrk = 2*m*m + wrkbl
547 minwrk = 2*m + n
548 ELSE IF( wntvs .AND. wntuas ) THEN
549*
550* Path 6t(N much larger than M, JOBU='S' or 'A',
551* JOBVT='S')
552*
553 wrkbl = m + lwork_cgelqf
554 wrkbl = max( wrkbl, m+lwork_cunglq_m )
555 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
556 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
557 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
558 maxwrk = m*m + wrkbl
559 minwrk = 2*m + n
560 ELSE IF( wntva .AND. wntun ) THEN
561*
562* Path 7t(N much larger than M, JOBU='N', JOBVT='A')
563*
564 wrkbl = m + lwork_cgelqf
565 wrkbl = max( wrkbl, m+lwork_cunglq_n )
566 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
567 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
568 maxwrk = m*m + wrkbl
569 minwrk = 2*m + n
570 ELSE IF( wntva .AND. wntuo ) THEN
571*
572* Path 8t(N much larger than M, JOBU='O', JOBVT='A')
573*
574 wrkbl = m + lwork_cgelqf
575 wrkbl = max( wrkbl, m+lwork_cunglq_n )
576 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
577 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
578 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
579 maxwrk = 2*m*m + wrkbl
580 minwrk = 2*m + n
581 ELSE IF( wntva .AND. wntuas ) THEN
582*
583* Path 9t(N much larger than M, JOBU='S' or 'A',
584* JOBVT='A')
585*
586 wrkbl = m + lwork_cgelqf
587 wrkbl = max( wrkbl, m+lwork_cunglq_n )
588 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
589 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
590 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
591 maxwrk = m*m + wrkbl
592 minwrk = 2*m + n
593 END IF
594 ELSE
595*
596* Path 10t(N greater than M, but not much larger)
597*
598 CALL cgebrd( m, n, a, lda, s, dum(1), cdum(1),
599 $ cdum(1), cdum(1), -1, ierr )
600 lwork_cgebrd = int( cdum(1) )
601 maxwrk = 2*m + lwork_cgebrd
602 IF( wntvs .OR. wntvo ) THEN
603* Compute space needed for CUNGBR P
604 CALL cungbr( 'P', m, n, m, a, n, cdum(1),
605 $ cdum(1), -1, ierr )
606 lwork_cungbr_p = int( cdum(1) )
607 maxwrk = max( maxwrk, 2*m+lwork_cungbr_p )
608 END IF
609 IF( wntva ) THEN
610 CALL cungbr( 'P', n, n, m, a, n, cdum(1),
611 $ cdum(1), -1, ierr )
612 lwork_cungbr_p = int( cdum(1) )
613 maxwrk = max( maxwrk, 2*m+lwork_cungbr_p )
614 END IF
615 IF( .NOT.wntun ) THEN
616 maxwrk = max( maxwrk, 2*m+lwork_cungbr_q )
617 END IF
618 minwrk = 2*m + n
619 END IF
620 END IF
621 maxwrk = max( minwrk, maxwrk )
622 work( 1 ) = sroundup_lwork(maxwrk)
623*
624 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
625 info = -13
626 END IF
627 END IF
628*
629 IF( info.NE.0 ) THEN
630 CALL xerbla( 'CGESVD', -info )
631 RETURN
632 ELSE IF( lquery ) THEN
633 RETURN
634 END IF
635*
636* Quick return if possible
637*
638 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
639 RETURN
640 END IF
641*
642* Get machine constants
643*
644 eps = slamch( 'P' )
645 smlnum = sqrt( slamch( 'S' ) ) / eps
646 bignum = one / smlnum
647*
648* Scale A if max element outside range [SMLNUM,BIGNUM]
649*
650 anrm = clange( 'M', m, n, a, lda, dum )
651 iscl = 0
652 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
653 iscl = 1
654 CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
655 ELSE IF( anrm.GT.bignum ) THEN
656 iscl = 1
657 CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
658 END IF
659*
660 IF( m.GE.n ) THEN
661*
662* A has at least as many rows as columns. If A has sufficiently
663* more rows than columns, first reduce using the QR
664* decomposition (if sufficient workspace available)
665*
666 IF( m.GE.mnthr ) THEN
667*
668 IF( wntun ) THEN
669*
670* Path 1 (M much larger than N, JOBU='N')
671* No left singular vectors to be computed
672*
673 itau = 1
674 iwork = itau + n
675*
676* Compute A=Q*R
677* (CWorkspace: need 2*N, prefer N+N*NB)
678* (RWorkspace: need 0)
679*
680 CALL cgeqrf( m, n, a, lda, work( itau ),
681 $ work( iwork ),
682 $ lwork-iwork+1, ierr )
683*
684* Zero out below R
685*
686 IF( n .GT. 1 ) THEN
687 CALL claset( 'L', n-1, n-1, czero, czero, a( 2,
688 $ 1 ),
689 $ lda )
690 END IF
691 ie = 1
692 itauq = 1
693 itaup = itauq + n
694 iwork = itaup + n
695*
696* Bidiagonalize R in A
697* (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
698* (RWorkspace: need N)
699*
700 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
701 $ work( itauq ),
702 $ work( itaup ), work( iwork ), lwork-iwork+1,
703 $ ierr )
704 ncvt = 0
705 IF( wntvo .OR. wntvas ) THEN
706*
707* If right singular vectors desired, generate P'.
708* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
709* (RWorkspace: 0)
710*
711 CALL cungbr( 'P', n, n, n, a, lda, work( itaup ),
712 $ work( iwork ), lwork-iwork+1, ierr )
713 ncvt = n
714 END IF
715 irwork = ie + n
716*
717* Perform bidiagonal QR iteration, computing right
718* singular vectors of A in A if desired
719* (CWorkspace: 0)
720* (RWorkspace: need BDSPAC)
721*
722 CALL cbdsqr( 'U', n, ncvt, 0, 0, s, rwork( ie ), a,
723 $ lda,
724 $ cdum, 1, cdum, 1, rwork( irwork ), info )
725*
726* If right singular vectors desired in VT, copy them there
727*
728 IF( wntvas )
729 $ CALL clacpy( 'F', n, n, a, lda, vt, ldvt )
730*
731 ELSE IF( wntuo .AND. wntvn ) THEN
732*
733* Path 2 (M much larger than N, JOBU='O', JOBVT='N')
734* N left singular vectors to be overwritten on A and
735* no right singular vectors to be computed
736*
737 IF( lwork.GE.n*n+3*n ) THEN
738*
739* Sufficient workspace for a fast algorithm
740*
741 ir = 1
742 IF( lwork.GE.max( wrkbl, lda*n )+lda*n ) THEN
743*
744* WORK(IU) is LDA by N, WORK(IR) is LDA by N
745*
746 ldwrku = lda
747 ldwrkr = lda
748 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n ) THEN
749*
750* WORK(IU) is LDA by N, WORK(IR) is N by N
751*
752 ldwrku = lda
753 ldwrkr = n
754 ELSE
755*
756* WORK(IU) is LDWRKU by N, WORK(IR) is N by N
757*
758 ldwrku = ( lwork-n*n ) / n
759 ldwrkr = n
760 END IF
761 itau = ir + ldwrkr*n
762 iwork = itau + n
763*
764* Compute A=Q*R
765* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
766* (RWorkspace: 0)
767*
768 CALL cgeqrf( m, n, a, lda, work( itau ),
769 $ work( iwork ), lwork-iwork+1, ierr )
770*
771* Copy R to WORK(IR) and zero out below it
772*
773 CALL clacpy( 'U', n, n, a, lda, work( ir ),
774 $ ldwrkr )
775 CALL claset( 'L', n-1, n-1, czero, czero,
776 $ work( ir+1 ), ldwrkr )
777*
778* Generate Q in A
779* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
780* (RWorkspace: 0)
781*
782 CALL cungqr( m, n, n, a, lda, work( itau ),
783 $ work( iwork ), lwork-iwork+1, ierr )
784 ie = 1
785 itauq = itau
786 itaup = itauq + n
787 iwork = itaup + n
788*
789* Bidiagonalize R in WORK(IR)
790* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
791* (RWorkspace: need N)
792*
793 CALL cgebrd( n, n, work( ir ), ldwrkr, s,
794 $ rwork( ie ),
795 $ work( itauq ), work( itaup ),
796 $ work( iwork ), lwork-iwork+1, ierr )
797*
798* Generate left vectors bidiagonalizing R
799* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
800* (RWorkspace: need 0)
801*
802 CALL cungbr( 'Q', n, n, n, work( ir ), ldwrkr,
803 $ work( itauq ), work( iwork ),
804 $ lwork-iwork+1, ierr )
805 irwork = ie + n
806*
807* Perform bidiagonal QR iteration, computing left
808* singular vectors of R in WORK(IR)
809* (CWorkspace: need N*N)
810* (RWorkspace: need BDSPAC)
811*
812 CALL cbdsqr( 'U', n, 0, n, 0, s, rwork( ie ), cdum,
813 $ 1,
814 $ work( ir ), ldwrkr, cdum, 1,
815 $ rwork( irwork ), info )
816 iu = itauq
817*
818* Multiply Q in A by left singular vectors of R in
819* WORK(IR), storing result in WORK(IU) and copying to A
820* (CWorkspace: need N*N+N, prefer N*N+M*N)
821* (RWorkspace: 0)
822*
823 DO 10 i = 1, m, ldwrku
824 chunk = min( m-i+1, ldwrku )
825 CALL cgemm( 'N', 'N', chunk, n, n, cone, a( i,
826 $ 1 ),
827 $ lda, work( ir ), ldwrkr, czero,
828 $ work( iu ), ldwrku )
829 CALL clacpy( 'F', chunk, n, work( iu ), ldwrku,
830 $ a( i, 1 ), lda )
831 10 CONTINUE
832*
833 ELSE
834*
835* Insufficient workspace for a fast algorithm
836*
837 ie = 1
838 itauq = 1
839 itaup = itauq + n
840 iwork = itaup + n
841*
842* Bidiagonalize A
843* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
844* (RWorkspace: N)
845*
846 CALL cgebrd( m, n, a, lda, s, rwork( ie ),
847 $ work( itauq ), work( itaup ),
848 $ work( iwork ), lwork-iwork+1, ierr )
849*
850* Generate left vectors bidiagonalizing A
851* (CWorkspace: need 3*N, prefer 2*N+N*NB)
852* (RWorkspace: 0)
853*
854 CALL cungbr( 'Q', m, n, n, a, lda, work( itauq ),
855 $ work( iwork ), lwork-iwork+1, ierr )
856 irwork = ie + n
857*
858* Perform bidiagonal QR iteration, computing left
859* singular vectors of A in A
860* (CWorkspace: need 0)
861* (RWorkspace: need BDSPAC)
862*
863 CALL cbdsqr( 'U', n, 0, m, 0, s, rwork( ie ), cdum,
864 $ 1,
865 $ a, lda, cdum, 1, rwork( irwork ), info )
866*
867 END IF
868*
869 ELSE IF( wntuo .AND. wntvas ) THEN
870*
871* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
872* N left singular vectors to be overwritten on A and
873* N right singular vectors to be computed in VT
874*
875 IF( lwork.GE.n*n+3*n ) THEN
876*
877* Sufficient workspace for a fast algorithm
878*
879 ir = 1
880 IF( lwork.GE.max( wrkbl, lda*n )+lda*n ) THEN
881*
882* WORK(IU) is LDA by N and WORK(IR) is LDA by N
883*
884 ldwrku = lda
885 ldwrkr = lda
886 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n ) THEN
887*
888* WORK(IU) is LDA by N and WORK(IR) is N by N
889*
890 ldwrku = lda
891 ldwrkr = n
892 ELSE
893*
894* WORK(IU) is LDWRKU by N and WORK(IR) is N by N
895*
896 ldwrku = ( lwork-n*n ) / n
897 ldwrkr = n
898 END IF
899 itau = ir + ldwrkr*n
900 iwork = itau + n
901*
902* Compute A=Q*R
903* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
904* (RWorkspace: 0)
905*
906 CALL cgeqrf( m, n, a, lda, work( itau ),
907 $ work( iwork ), lwork-iwork+1, ierr )
908*
909* Copy R to VT, zeroing out below it
910*
911 CALL clacpy( 'U', n, n, a, lda, vt, ldvt )
912 IF( n.GT.1 )
913 $ CALL claset( 'L', n-1, n-1, czero, czero,
914 $ vt( 2, 1 ), ldvt )
915*
916* Generate Q in A
917* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
918* (RWorkspace: 0)
919*
920 CALL cungqr( m, n, n, a, lda, work( itau ),
921 $ work( iwork ), lwork-iwork+1, ierr )
922 ie = 1
923 itauq = itau
924 itaup = itauq + n
925 iwork = itaup + n
926*
927* Bidiagonalize R in VT, copying result to WORK(IR)
928* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
929* (RWorkspace: need N)
930*
931 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
932 $ work( itauq ), work( itaup ),
933 $ work( iwork ), lwork-iwork+1, ierr )
934 CALL clacpy( 'L', n, n, vt, ldvt, work( ir ),
935 $ ldwrkr )
936*
937* Generate left vectors bidiagonalizing R in WORK(IR)
938* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
939* (RWorkspace: 0)
940*
941 CALL cungbr( 'Q', n, n, n, work( ir ), ldwrkr,
942 $ work( itauq ), work( iwork ),
943 $ lwork-iwork+1, ierr )
944*
945* Generate right vectors bidiagonalizing R in VT
946* (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB)
947* (RWorkspace: 0)
948*
949 CALL cungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
950 $ work( iwork ), lwork-iwork+1, ierr )
951 irwork = ie + n
952*
953* Perform bidiagonal QR iteration, computing left
954* singular vectors of R in WORK(IR) and computing right
955* singular vectors of R in VT
956* (CWorkspace: need N*N)
957* (RWorkspace: need BDSPAC)
958*
959 CALL cbdsqr( 'U', n, n, n, 0, s, rwork( ie ), vt,
960 $ ldvt, work( ir ), ldwrkr, cdum, 1,
961 $ rwork( irwork ), info )
962 iu = itauq
963*
964* Multiply Q in A by left singular vectors of R in
965* WORK(IR), storing result in WORK(IU) and copying to A
966* (CWorkspace: need N*N+N, prefer N*N+M*N)
967* (RWorkspace: 0)
968*
969 DO 20 i = 1, m, ldwrku
970 chunk = min( m-i+1, ldwrku )
971 CALL cgemm( 'N', 'N', chunk, n, n, cone, a( i,
972 $ 1 ),
973 $ lda, work( ir ), ldwrkr, czero,
974 $ work( iu ), ldwrku )
975 CALL clacpy( 'F', chunk, n, work( iu ), ldwrku,
976 $ a( i, 1 ), lda )
977 20 CONTINUE
978*
979 ELSE
980*
981* Insufficient workspace for a fast algorithm
982*
983 itau = 1
984 iwork = itau + n
985*
986* Compute A=Q*R
987* (CWorkspace: need 2*N, prefer N+N*NB)
988* (RWorkspace: 0)
989*
990 CALL cgeqrf( m, n, a, lda, work( itau ),
991 $ work( iwork ), lwork-iwork+1, ierr )
992*
993* Copy R to VT, zeroing out below it
994*
995 CALL clacpy( 'U', n, n, a, lda, vt, ldvt )
996 IF( n.GT.1 )
997 $ CALL claset( 'L', n-1, n-1, czero, czero,
998 $ vt( 2, 1 ), ldvt )
999*
1000* Generate Q in A
1001* (CWorkspace: need 2*N, prefer N+N*NB)
1002* (RWorkspace: 0)
1003*
1004 CALL cungqr( m, n, n, a, lda, work( itau ),
1005 $ work( iwork ), lwork-iwork+1, ierr )
1006 ie = 1
1007 itauq = itau
1008 itaup = itauq + n
1009 iwork = itaup + n
1010*
1011* Bidiagonalize R in VT
1012* (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1013* (RWorkspace: N)
1014*
1015 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1016 $ work( itauq ), work( itaup ),
1017 $ work( iwork ), lwork-iwork+1, ierr )
1018*
1019* Multiply Q in A by left vectors bidiagonalizing R
1020* (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1021* (RWorkspace: 0)
1022*
1023 CALL cunmbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
1024 $ work( itauq ), a, lda, work( iwork ),
1025 $ lwork-iwork+1, ierr )
1026*
1027* Generate right vectors bidiagonalizing R in VT
1028* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1029* (RWorkspace: 0)
1030*
1031 CALL cungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1032 $ work( iwork ), lwork-iwork+1, ierr )
1033 irwork = ie + n
1034*
1035* Perform bidiagonal QR iteration, computing left
1036* singular vectors of A in A and computing right
1037* singular vectors of A in VT
1038* (CWorkspace: 0)
1039* (RWorkspace: need BDSPAC)
1040*
1041 CALL cbdsqr( 'U', n, n, m, 0, s, rwork( ie ), vt,
1042 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
1043 $ info )
1044*
1045 END IF
1046*
1047 ELSE IF( wntus ) THEN
1048*
1049 IF( wntvn ) THEN
1050*
1051* Path 4 (M much larger than N, JOBU='S', JOBVT='N')
1052* N left singular vectors to be computed in U and
1053* no right singular vectors to be computed
1054*
1055 IF( lwork.GE.n*n+3*n ) THEN
1056*
1057* Sufficient workspace for a fast algorithm
1058*
1059 ir = 1
1060 IF( lwork.GE.wrkbl+lda*n ) THEN
1061*
1062* WORK(IR) is LDA by N
1063*
1064 ldwrkr = lda
1065 ELSE
1066*
1067* WORK(IR) is N by N
1068*
1069 ldwrkr = n
1070 END IF
1071 itau = ir + ldwrkr*n
1072 iwork = itau + n
1073*
1074* Compute A=Q*R
1075* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1076* (RWorkspace: 0)
1077*
1078 CALL cgeqrf( m, n, a, lda, work( itau ),
1079 $ work( iwork ), lwork-iwork+1, ierr )
1080*
1081* Copy R to WORK(IR), zeroing out below it
1082*
1083 CALL clacpy( 'U', n, n, a, lda, work( ir ),
1084 $ ldwrkr )
1085 CALL claset( 'L', n-1, n-1, czero, czero,
1086 $ work( ir+1 ), ldwrkr )
1087*
1088* Generate Q in A
1089* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1090* (RWorkspace: 0)
1091*
1092 CALL cungqr( m, n, n, a, lda, work( itau ),
1093 $ work( iwork ), lwork-iwork+1, ierr )
1094 ie = 1
1095 itauq = itau
1096 itaup = itauq + n
1097 iwork = itaup + n
1098*
1099* Bidiagonalize R in WORK(IR)
1100* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1101* (RWorkspace: need N)
1102*
1103 CALL cgebrd( n, n, work( ir ), ldwrkr, s,
1104 $ rwork( ie ), work( itauq ),
1105 $ work( itaup ), work( iwork ),
1106 $ lwork-iwork+1, ierr )
1107*
1108* Generate left vectors bidiagonalizing R in WORK(IR)
1109* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1110* (RWorkspace: 0)
1111*
1112 CALL cungbr( 'Q', n, n, n, work( ir ), ldwrkr,
1113 $ work( itauq ), work( iwork ),
1114 $ lwork-iwork+1, ierr )
1115 irwork = ie + n
1116*
1117* Perform bidiagonal QR iteration, computing left
1118* singular vectors of R in WORK(IR)
1119* (CWorkspace: need N*N)
1120* (RWorkspace: need BDSPAC)
1121*
1122 CALL cbdsqr( 'U', n, 0, n, 0, s, rwork( ie ),
1123 $ cdum,
1124 $ 1, work( ir ), ldwrkr, cdum, 1,
1125 $ rwork( irwork ), info )
1126*
1127* Multiply Q in A by left singular vectors of R in
1128* WORK(IR), storing result in U
1129* (CWorkspace: need N*N)
1130* (RWorkspace: 0)
1131*
1132 CALL cgemm( 'N', 'N', m, n, n, cone, a, lda,
1133 $ work( ir ), ldwrkr, czero, u, ldu )
1134*
1135 ELSE
1136*
1137* Insufficient workspace for a fast algorithm
1138*
1139 itau = 1
1140 iwork = itau + n
1141*
1142* Compute A=Q*R, copying result to U
1143* (CWorkspace: need 2*N, prefer N+N*NB)
1144* (RWorkspace: 0)
1145*
1146 CALL cgeqrf( m, n, a, lda, work( itau ),
1147 $ work( iwork ), lwork-iwork+1, ierr )
1148 CALL clacpy( 'L', m, n, a, lda, u, ldu )
1149*
1150* Generate Q in U
1151* (CWorkspace: need 2*N, prefer N+N*NB)
1152* (RWorkspace: 0)
1153*
1154 CALL cungqr( m, n, n, u, ldu, work( itau ),
1155 $ work( iwork ), lwork-iwork+1, ierr )
1156 ie = 1
1157 itauq = itau
1158 itaup = itauq + n
1159 iwork = itaup + n
1160*
1161* Zero out below R in A
1162*
1163 IF( n .GT. 1 ) THEN
1164 CALL claset( 'L', n-1, n-1, czero, czero,
1165 $ a( 2, 1 ), lda )
1166 END IF
1167*
1168* Bidiagonalize R in A
1169* (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1170* (RWorkspace: need N)
1171*
1172 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1173 $ work( itauq ), work( itaup ),
1174 $ work( iwork ), lwork-iwork+1, ierr )
1175*
1176* Multiply Q in U by left vectors bidiagonalizing R
1177* (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1178* (RWorkspace: 0)
1179*
1180 CALL cunmbr( 'Q', 'R', 'N', m, n, n, a, lda,
1181 $ work( itauq ), u, ldu, work( iwork ),
1182 $ lwork-iwork+1, ierr )
1183 irwork = ie + n
1184*
1185* Perform bidiagonal QR iteration, computing left
1186* singular vectors of A in U
1187* (CWorkspace: 0)
1188* (RWorkspace: need BDSPAC)
1189*
1190 CALL cbdsqr( 'U', n, 0, m, 0, s, rwork( ie ),
1191 $ cdum,
1192 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1193 $ info )
1194*
1195 END IF
1196*
1197 ELSE IF( wntvo ) THEN
1198*
1199* Path 5 (M much larger than N, JOBU='S', JOBVT='O')
1200* N left singular vectors to be computed in U and
1201* N right singular vectors to be overwritten on A
1202*
1203 IF( lwork.GE.2*n*n+3*n ) THEN
1204*
1205* Sufficient workspace for a fast algorithm
1206*
1207 iu = 1
1208 IF( lwork.GE.wrkbl+2*lda*n ) THEN
1209*
1210* WORK(IU) is LDA by N and WORK(IR) is LDA by N
1211*
1212 ldwrku = lda
1213 ir = iu + ldwrku*n
1214 ldwrkr = lda
1215 ELSE IF( lwork.GE.wrkbl+( lda+n )*n ) THEN
1216*
1217* WORK(IU) is LDA by N and WORK(IR) is N by N
1218*
1219 ldwrku = lda
1220 ir = iu + ldwrku*n
1221 ldwrkr = n
1222 ELSE
1223*
1224* WORK(IU) is N by N and WORK(IR) is N by N
1225*
1226 ldwrku = n
1227 ir = iu + ldwrku*n
1228 ldwrkr = n
1229 END IF
1230 itau = ir + ldwrkr*n
1231 iwork = itau + n
1232*
1233* Compute A=Q*R
1234* (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1235* (RWorkspace: 0)
1236*
1237 CALL cgeqrf( m, n, a, lda, work( itau ),
1238 $ work( iwork ), lwork-iwork+1, ierr )
1239*
1240* Copy R to WORK(IU), zeroing out below it
1241*
1242 CALL clacpy( 'U', n, n, a, lda, work( iu ),
1243 $ ldwrku )
1244 CALL claset( 'L', n-1, n-1, czero, czero,
1245 $ work( iu+1 ), ldwrku )
1246*
1247* Generate Q in A
1248* (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1249* (RWorkspace: 0)
1250*
1251 CALL cungqr( m, n, n, a, lda, work( itau ),
1252 $ work( iwork ), lwork-iwork+1, ierr )
1253 ie = 1
1254 itauq = itau
1255 itaup = itauq + n
1256 iwork = itaup + n
1257*
1258* Bidiagonalize R in WORK(IU), copying result to
1259* WORK(IR)
1260* (CWorkspace: need 2*N*N+3*N,
1261* prefer 2*N*N+2*N+2*N*NB)
1262* (RWorkspace: need N)
1263*
1264 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1265 $ rwork( ie ), work( itauq ),
1266 $ work( itaup ), work( iwork ),
1267 $ lwork-iwork+1, ierr )
1268 CALL clacpy( 'U', n, n, work( iu ), ldwrku,
1269 $ work( ir ), ldwrkr )
1270*
1271* Generate left bidiagonalizing vectors in WORK(IU)
1272* (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
1273* (RWorkspace: 0)
1274*
1275 CALL cungbr( 'Q', n, n, n, work( iu ), ldwrku,
1276 $ work( itauq ), work( iwork ),
1277 $ lwork-iwork+1, ierr )
1278*
1279* Generate right bidiagonalizing vectors in WORK(IR)
1280* (CWorkspace: need 2*N*N+3*N-1,
1281* prefer 2*N*N+2*N+(N-1)*NB)
1282* (RWorkspace: 0)
1283*
1284 CALL cungbr( 'P', n, n, n, work( ir ), ldwrkr,
1285 $ work( itaup ), work( iwork ),
1286 $ lwork-iwork+1, ierr )
1287 irwork = ie + n
1288*
1289* Perform bidiagonal QR iteration, computing left
1290* singular vectors of R in WORK(IU) and computing
1291* right singular vectors of R in WORK(IR)
1292* (CWorkspace: need 2*N*N)
1293* (RWorkspace: need BDSPAC)
1294*
1295 CALL cbdsqr( 'U', n, n, n, 0, s, rwork( ie ),
1296 $ work( ir ), ldwrkr, work( iu ),
1297 $ ldwrku, cdum, 1, rwork( irwork ),
1298 $ info )
1299*
1300* Multiply Q in A by left singular vectors of R in
1301* WORK(IU), storing result in U
1302* (CWorkspace: need N*N)
1303* (RWorkspace: 0)
1304*
1305 CALL cgemm( 'N', 'N', m, n, n, cone, a, lda,
1306 $ work( iu ), ldwrku, czero, u, ldu )
1307*
1308* Copy right singular vectors of R to A
1309* (CWorkspace: need N*N)
1310* (RWorkspace: 0)
1311*
1312 CALL clacpy( 'F', n, n, work( ir ), ldwrkr, a,
1313 $ lda )
1314*
1315 ELSE
1316*
1317* Insufficient workspace for a fast algorithm
1318*
1319 itau = 1
1320 iwork = itau + n
1321*
1322* Compute A=Q*R, copying result to U
1323* (CWorkspace: need 2*N, prefer N+N*NB)
1324* (RWorkspace: 0)
1325*
1326 CALL cgeqrf( m, n, a, lda, work( itau ),
1327 $ work( iwork ), lwork-iwork+1, ierr )
1328 CALL clacpy( 'L', m, n, a, lda, u, ldu )
1329*
1330* Generate Q in U
1331* (CWorkspace: need 2*N, prefer N+N*NB)
1332* (RWorkspace: 0)
1333*
1334 CALL cungqr( m, n, n, u, ldu, work( itau ),
1335 $ work( iwork ), lwork-iwork+1, ierr )
1336 ie = 1
1337 itauq = itau
1338 itaup = itauq + n
1339 iwork = itaup + n
1340*
1341* Zero out below R in A
1342*
1343 IF( n .GT. 1 ) THEN
1344 CALL claset( 'L', n-1, n-1, czero, czero,
1345 $ a( 2, 1 ), lda )
1346 END IF
1347*
1348* Bidiagonalize R in A
1349* (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1350* (RWorkspace: need N)
1351*
1352 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1353 $ work( itauq ), work( itaup ),
1354 $ work( iwork ), lwork-iwork+1, ierr )
1355*
1356* Multiply Q in U by left vectors bidiagonalizing R
1357* (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1358* (RWorkspace: 0)
1359*
1360 CALL cunmbr( 'Q', 'R', 'N', m, n, n, a, lda,
1361 $ work( itauq ), u, ldu, work( iwork ),
1362 $ lwork-iwork+1, ierr )
1363*
1364* Generate right vectors bidiagonalizing R in A
1365* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1366* (RWorkspace: 0)
1367*
1368 CALL cungbr( 'P', n, n, n, a, lda,
1369 $ work( itaup ),
1370 $ work( iwork ), lwork-iwork+1, ierr )
1371 irwork = ie + n
1372*
1373* Perform bidiagonal QR iteration, computing left
1374* singular vectors of A in U and computing right
1375* singular vectors of A in A
1376* (CWorkspace: 0)
1377* (RWorkspace: need BDSPAC)
1378*
1379 CALL cbdsqr( 'U', n, n, m, 0, s, rwork( ie ), a,
1380 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1381 $ info )
1382*
1383 END IF
1384*
1385 ELSE IF( wntvas ) THEN
1386*
1387* Path 6 (M much larger than N, JOBU='S', JOBVT='S'
1388* or 'A')
1389* N left singular vectors to be computed in U and
1390* N right singular vectors to be computed in VT
1391*
1392 IF( lwork.GE.n*n+3*n ) THEN
1393*
1394* Sufficient workspace for a fast algorithm
1395*
1396 iu = 1
1397 IF( lwork.GE.wrkbl+lda*n ) THEN
1398*
1399* WORK(IU) is LDA by N
1400*
1401 ldwrku = lda
1402 ELSE
1403*
1404* WORK(IU) is N by N
1405*
1406 ldwrku = n
1407 END IF
1408 itau = iu + ldwrku*n
1409 iwork = itau + n
1410*
1411* Compute A=Q*R
1412* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1413* (RWorkspace: 0)
1414*
1415 CALL cgeqrf( m, n, a, lda, work( itau ),
1416 $ work( iwork ), lwork-iwork+1, ierr )
1417*
1418* Copy R to WORK(IU), zeroing out below it
1419*
1420 CALL clacpy( 'U', n, n, a, lda, work( iu ),
1421 $ ldwrku )
1422 CALL claset( 'L', n-1, n-1, czero, czero,
1423 $ work( iu+1 ), ldwrku )
1424*
1425* Generate Q in A
1426* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1427* (RWorkspace: 0)
1428*
1429 CALL cungqr( m, n, n, a, lda, work( itau ),
1430 $ work( iwork ), lwork-iwork+1, ierr )
1431 ie = 1
1432 itauq = itau
1433 itaup = itauq + n
1434 iwork = itaup + n
1435*
1436* Bidiagonalize R in WORK(IU), copying result to VT
1437* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1438* (RWorkspace: need N)
1439*
1440 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1441 $ rwork( ie ), work( itauq ),
1442 $ work( itaup ), work( iwork ),
1443 $ lwork-iwork+1, ierr )
1444 CALL clacpy( 'U', n, n, work( iu ), ldwrku, vt,
1445 $ ldvt )
1446*
1447* Generate left bidiagonalizing vectors in WORK(IU)
1448* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1449* (RWorkspace: 0)
1450*
1451 CALL cungbr( 'Q', n, n, n, work( iu ), ldwrku,
1452 $ work( itauq ), work( iwork ),
1453 $ lwork-iwork+1, ierr )
1454*
1455* Generate right bidiagonalizing vectors in VT
1456* (CWorkspace: need N*N+3*N-1,
1457* prefer N*N+2*N+(N-1)*NB)
1458* (RWorkspace: 0)
1459*
1460 CALL cungbr( 'P', n, n, n, vt, ldvt,
1461 $ work( itaup ),
1462 $ work( iwork ), lwork-iwork+1, ierr )
1463 irwork = ie + n
1464*
1465* Perform bidiagonal QR iteration, computing left
1466* singular vectors of R in WORK(IU) and computing
1467* right singular vectors of R in VT
1468* (CWorkspace: need N*N)
1469* (RWorkspace: need BDSPAC)
1470*
1471 CALL cbdsqr( 'U', n, n, n, 0, s, rwork( ie ),
1472 $ vt,
1473 $ ldvt, work( iu ), ldwrku, cdum, 1,
1474 $ rwork( irwork ), info )
1475*
1476* Multiply Q in A by left singular vectors of R in
1477* WORK(IU), storing result in U
1478* (CWorkspace: need N*N)
1479* (RWorkspace: 0)
1480*
1481 CALL cgemm( 'N', 'N', m, n, n, cone, a, lda,
1482 $ work( iu ), ldwrku, czero, u, ldu )
1483*
1484 ELSE
1485*
1486* Insufficient workspace for a fast algorithm
1487*
1488 itau = 1
1489 iwork = itau + n
1490*
1491* Compute A=Q*R, copying result to U
1492* (CWorkspace: need 2*N, prefer N+N*NB)
1493* (RWorkspace: 0)
1494*
1495 CALL cgeqrf( m, n, a, lda, work( itau ),
1496 $ work( iwork ), lwork-iwork+1, ierr )
1497 CALL clacpy( 'L', m, n, a, lda, u, ldu )
1498*
1499* Generate Q in U
1500* (CWorkspace: need 2*N, prefer N+N*NB)
1501* (RWorkspace: 0)
1502*
1503 CALL cungqr( m, n, n, u, ldu, work( itau ),
1504 $ work( iwork ), lwork-iwork+1, ierr )
1505*
1506* Copy R to VT, zeroing out below it
1507*
1508 CALL clacpy( 'U', n, n, a, lda, vt, ldvt )
1509 IF( n.GT.1 )
1510 $ CALL claset( 'L', n-1, n-1, czero, czero,
1511 $ vt( 2, 1 ), ldvt )
1512 ie = 1
1513 itauq = itau
1514 itaup = itauq + n
1515 iwork = itaup + n
1516*
1517* Bidiagonalize R in VT
1518* (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1519* (RWorkspace: need N)
1520*
1521 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1522 $ work( itauq ), work( itaup ),
1523 $ work( iwork ), lwork-iwork+1, ierr )
1524*
1525* Multiply Q in U by left bidiagonalizing vectors
1526* in VT
1527* (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1528* (RWorkspace: 0)
1529*
1530 CALL cunmbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
1531 $ work( itauq ), u, ldu, work( iwork ),
1532 $ lwork-iwork+1, ierr )
1533*
1534* Generate right bidiagonalizing vectors in VT
1535* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1536* (RWorkspace: 0)
1537*
1538 CALL cungbr( 'P', n, n, n, vt, ldvt,
1539 $ work( itaup ),
1540 $ work( iwork ), lwork-iwork+1, ierr )
1541 irwork = ie + n
1542*
1543* Perform bidiagonal QR iteration, computing left
1544* singular vectors of A in U and computing right
1545* singular vectors of A in VT
1546* (CWorkspace: 0)
1547* (RWorkspace: need BDSPAC)
1548*
1549 CALL cbdsqr( 'U', n, n, m, 0, s, rwork( ie ),
1550 $ vt,
1551 $ ldvt, u, ldu, cdum, 1,
1552 $ rwork( irwork ), info )
1553*
1554 END IF
1555*
1556 END IF
1557*
1558 ELSE IF( wntua ) THEN
1559*
1560 IF( wntvn ) THEN
1561*
1562* Path 7 (M much larger than N, JOBU='A', JOBVT='N')
1563* M left singular vectors to be computed in U and
1564* no right singular vectors to be computed
1565*
1566 IF( lwork.GE.n*n+max( n+m, 3*n ) ) THEN
1567*
1568* Sufficient workspace for a fast algorithm
1569*
1570 ir = 1
1571 IF( lwork.GE.wrkbl+lda*n ) THEN
1572*
1573* WORK(IR) is LDA by N
1574*
1575 ldwrkr = lda
1576 ELSE
1577*
1578* WORK(IR) is N by N
1579*
1580 ldwrkr = n
1581 END IF
1582 itau = ir + ldwrkr*n
1583 iwork = itau + n
1584*
1585* Compute A=Q*R, copying result to U
1586* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1587* (RWorkspace: 0)
1588*
1589 CALL cgeqrf( m, n, a, lda, work( itau ),
1590 $ work( iwork ), lwork-iwork+1, ierr )
1591 CALL clacpy( 'L', m, n, a, lda, u, ldu )
1592*
1593* Copy R to WORK(IR), zeroing out below it
1594*
1595 CALL clacpy( 'U', n, n, a, lda, work( ir ),
1596 $ ldwrkr )
1597 CALL claset( 'L', n-1, n-1, czero, czero,
1598 $ work( ir+1 ), ldwrkr )
1599*
1600* Generate Q in U
1601* (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
1602* (RWorkspace: 0)
1603*
1604 CALL cungqr( m, m, n, u, ldu, work( itau ),
1605 $ work( iwork ), lwork-iwork+1, ierr )
1606 ie = 1
1607 itauq = itau
1608 itaup = itauq + n
1609 iwork = itaup + n
1610*
1611* Bidiagonalize R in WORK(IR)
1612* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1613* (RWorkspace: need N)
1614*
1615 CALL cgebrd( n, n, work( ir ), ldwrkr, s,
1616 $ rwork( ie ), work( itauq ),
1617 $ work( itaup ), work( iwork ),
1618 $ lwork-iwork+1, ierr )
1619*
1620* Generate left bidiagonalizing vectors in WORK(IR)
1621* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1622* (RWorkspace: 0)
1623*
1624 CALL cungbr( 'Q', n, n, n, work( ir ), ldwrkr,
1625 $ work( itauq ), work( iwork ),
1626 $ lwork-iwork+1, ierr )
1627 irwork = ie + n
1628*
1629* Perform bidiagonal QR iteration, computing left
1630* singular vectors of R in WORK(IR)
1631* (CWorkspace: need N*N)
1632* (RWorkspace: need BDSPAC)
1633*
1634 CALL cbdsqr( 'U', n, 0, n, 0, s, rwork( ie ),
1635 $ cdum,
1636 $ 1, work( ir ), ldwrkr, cdum, 1,
1637 $ rwork( irwork ), info )
1638*
1639* Multiply Q in U by left singular vectors of R in
1640* WORK(IR), storing result in A
1641* (CWorkspace: need N*N)
1642* (RWorkspace: 0)
1643*
1644 CALL cgemm( 'N', 'N', m, n, n, cone, u, ldu,
1645 $ work( ir ), ldwrkr, czero, a, lda )
1646*
1647* Copy left singular vectors of A from A to U
1648*
1649 CALL clacpy( 'F', m, n, a, lda, u, ldu )
1650*
1651 ELSE
1652*
1653* Insufficient workspace for a fast algorithm
1654*
1655 itau = 1
1656 iwork = itau + n
1657*
1658* Compute A=Q*R, copying result to U
1659* (CWorkspace: need 2*N, prefer N+N*NB)
1660* (RWorkspace: 0)
1661*
1662 CALL cgeqrf( m, n, a, lda, work( itau ),
1663 $ work( iwork ), lwork-iwork+1, ierr )
1664 CALL clacpy( 'L', m, n, a, lda, u, ldu )
1665*
1666* Generate Q in U
1667* (CWorkspace: need N+M, prefer N+M*NB)
1668* (RWorkspace: 0)
1669*
1670 CALL cungqr( m, m, n, u, ldu, work( itau ),
1671 $ work( iwork ), lwork-iwork+1, ierr )
1672 ie = 1
1673 itauq = itau
1674 itaup = itauq + n
1675 iwork = itaup + n
1676*
1677* Zero out below R in A
1678*
1679 IF( n .GT. 1 ) THEN
1680 CALL claset( 'L', n-1, n-1, czero, czero,
1681 $ a( 2, 1 ), lda )
1682 END IF
1683*
1684* Bidiagonalize R in A
1685* (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1686* (RWorkspace: need N)
1687*
1688 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1689 $ work( itauq ), work( itaup ),
1690 $ work( iwork ), lwork-iwork+1, ierr )
1691*
1692* Multiply Q in U by left bidiagonalizing vectors
1693* in A
1694* (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1695* (RWorkspace: 0)
1696*
1697 CALL cunmbr( 'Q', 'R', 'N', m, n, n, a, lda,
1698 $ work( itauq ), u, ldu, work( iwork ),
1699 $ lwork-iwork+1, ierr )
1700 irwork = ie + n
1701*
1702* Perform bidiagonal QR iteration, computing left
1703* singular vectors of A in U
1704* (CWorkspace: 0)
1705* (RWorkspace: need BDSPAC)
1706*
1707 CALL cbdsqr( 'U', n, 0, m, 0, s, rwork( ie ),
1708 $ cdum,
1709 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1710 $ info )
1711*
1712 END IF
1713*
1714 ELSE IF( wntvo ) THEN
1715*
1716* Path 8 (M much larger than N, JOBU='A', JOBVT='O')
1717* M left singular vectors to be computed in U and
1718* N right singular vectors to be overwritten on A
1719*
1720 IF( lwork.GE.2*n*n+max( n+m, 3*n ) ) THEN
1721*
1722* Sufficient workspace for a fast algorithm
1723*
1724 iu = 1
1725 IF( lwork.GE.wrkbl+2*lda*n ) THEN
1726*
1727* WORK(IU) is LDA by N and WORK(IR) is LDA by N
1728*
1729 ldwrku = lda
1730 ir = iu + ldwrku*n
1731 ldwrkr = lda
1732 ELSE IF( lwork.GE.wrkbl+( lda+n )*n ) THEN
1733*
1734* WORK(IU) is LDA by N and WORK(IR) is N by N
1735*
1736 ldwrku = lda
1737 ir = iu + ldwrku*n
1738 ldwrkr = n
1739 ELSE
1740*
1741* WORK(IU) is N by N and WORK(IR) is N by N
1742*
1743 ldwrku = n
1744 ir = iu + ldwrku*n
1745 ldwrkr = n
1746 END IF
1747 itau = ir + ldwrkr*n
1748 iwork = itau + n
1749*
1750* Compute A=Q*R, copying result to U
1751* (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1752* (RWorkspace: 0)
1753*
1754 CALL cgeqrf( m, n, a, lda, work( itau ),
1755 $ work( iwork ), lwork-iwork+1, ierr )
1756 CALL clacpy( 'L', m, n, a, lda, u, ldu )
1757*
1758* Generate Q in U
1759* (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
1760* (RWorkspace: 0)
1761*
1762 CALL cungqr( m, m, n, u, ldu, work( itau ),
1763 $ work( iwork ), lwork-iwork+1, ierr )
1764*
1765* Copy R to WORK(IU), zeroing out below it
1766*
1767 CALL clacpy( 'U', n, n, a, lda, work( iu ),
1768 $ ldwrku )
1769 CALL claset( 'L', n-1, n-1, czero, czero,
1770 $ work( iu+1 ), ldwrku )
1771 ie = 1
1772 itauq = itau
1773 itaup = itauq + n
1774 iwork = itaup + n
1775*
1776* Bidiagonalize R in WORK(IU), copying result to
1777* WORK(IR)
1778* (CWorkspace: need 2*N*N+3*N,
1779* prefer 2*N*N+2*N+2*N*NB)
1780* (RWorkspace: need N)
1781*
1782 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1783 $ rwork( ie ), work( itauq ),
1784 $ work( itaup ), work( iwork ),
1785 $ lwork-iwork+1, ierr )
1786 CALL clacpy( 'U', n, n, work( iu ), ldwrku,
1787 $ work( ir ), ldwrkr )
1788*
1789* Generate left bidiagonalizing vectors in WORK(IU)
1790* (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
1791* (RWorkspace: 0)
1792*
1793 CALL cungbr( 'Q', n, n, n, work( iu ), ldwrku,
1794 $ work( itauq ), work( iwork ),
1795 $ lwork-iwork+1, ierr )
1796*
1797* Generate right bidiagonalizing vectors in WORK(IR)
1798* (CWorkspace: need 2*N*N+3*N-1,
1799* prefer 2*N*N+2*N+(N-1)*NB)
1800* (RWorkspace: 0)
1801*
1802 CALL cungbr( 'P', n, n, n, work( ir ), ldwrkr,
1803 $ work( itaup ), work( iwork ),
1804 $ lwork-iwork+1, ierr )
1805 irwork = ie + n
1806*
1807* Perform bidiagonal QR iteration, computing left
1808* singular vectors of R in WORK(IU) and computing
1809* right singular vectors of R in WORK(IR)
1810* (CWorkspace: need 2*N*N)
1811* (RWorkspace: need BDSPAC)
1812*
1813 CALL cbdsqr( 'U', n, n, n, 0, s, rwork( ie ),
1814 $ work( ir ), ldwrkr, work( iu ),
1815 $ ldwrku, cdum, 1, rwork( irwork ),
1816 $ info )
1817*
1818* Multiply Q in U by left singular vectors of R in
1819* WORK(IU), storing result in A
1820* (CWorkspace: need N*N)
1821* (RWorkspace: 0)
1822*
1823 CALL cgemm( 'N', 'N', m, n, n, cone, u, ldu,
1824 $ work( iu ), ldwrku, czero, a, lda )
1825*
1826* Copy left singular vectors of A from A to U
1827*
1828 CALL clacpy( 'F', m, n, a, lda, u, ldu )
1829*
1830* Copy right singular vectors of R from WORK(IR) to A
1831*
1832 CALL clacpy( 'F', n, n, work( ir ), ldwrkr, a,
1833 $ lda )
1834*
1835 ELSE
1836*
1837* Insufficient workspace for a fast algorithm
1838*
1839 itau = 1
1840 iwork = itau + n
1841*
1842* Compute A=Q*R, copying result to U
1843* (CWorkspace: need 2*N, prefer N+N*NB)
1844* (RWorkspace: 0)
1845*
1846 CALL cgeqrf( m, n, a, lda, work( itau ),
1847 $ work( iwork ), lwork-iwork+1, ierr )
1848 CALL clacpy( 'L', m, n, a, lda, u, ldu )
1849*
1850* Generate Q in U
1851* (CWorkspace: need N+M, prefer N+M*NB)
1852* (RWorkspace: 0)
1853*
1854 CALL cungqr( m, m, n, u, ldu, work( itau ),
1855 $ work( iwork ), lwork-iwork+1, ierr )
1856 ie = 1
1857 itauq = itau
1858 itaup = itauq + n
1859 iwork = itaup + n
1860*
1861* Zero out below R in A
1862*
1863 IF( n .GT. 1 ) THEN
1864 CALL claset( 'L', n-1, n-1, czero, czero,
1865 $ a( 2, 1 ), lda )
1866 END IF
1867*
1868* Bidiagonalize R in A
1869* (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1870* (RWorkspace: need N)
1871*
1872 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1873 $ work( itauq ), work( itaup ),
1874 $ work( iwork ), lwork-iwork+1, ierr )
1875*
1876* Multiply Q in U by left bidiagonalizing vectors
1877* in A
1878* (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1879* (RWorkspace: 0)
1880*
1881 CALL cunmbr( 'Q', 'R', 'N', m, n, n, a, lda,
1882 $ work( itauq ), u, ldu, work( iwork ),
1883 $ lwork-iwork+1, ierr )
1884*
1885* Generate right bidiagonalizing vectors in A
1886* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1887* (RWorkspace: 0)
1888*
1889 CALL cungbr( 'P', n, n, n, a, lda,
1890 $ work( itaup ),
1891 $ work( iwork ), lwork-iwork+1, ierr )
1892 irwork = ie + n
1893*
1894* Perform bidiagonal QR iteration, computing left
1895* singular vectors of A in U and computing right
1896* singular vectors of A in A
1897* (CWorkspace: 0)
1898* (RWorkspace: need BDSPAC)
1899*
1900 CALL cbdsqr( 'U', n, n, m, 0, s, rwork( ie ), a,
1901 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1902 $ info )
1903*
1904 END IF
1905*
1906 ELSE IF( wntvas ) THEN
1907*
1908* Path 9 (M much larger than N, JOBU='A', JOBVT='S'
1909* or 'A')
1910* M left singular vectors to be computed in U and
1911* N right singular vectors to be computed in VT
1912*
1913 IF( lwork.GE.n*n+max( n+m, 3*n ) ) THEN
1914*
1915* Sufficient workspace for a fast algorithm
1916*
1917 iu = 1
1918 IF( lwork.GE.wrkbl+lda*n ) THEN
1919*
1920* WORK(IU) is LDA by N
1921*
1922 ldwrku = lda
1923 ELSE
1924*
1925* WORK(IU) is N by N
1926*
1927 ldwrku = n
1928 END IF
1929 itau = iu + ldwrku*n
1930 iwork = itau + n
1931*
1932* Compute A=Q*R, copying result to U
1933* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1934* (RWorkspace: 0)
1935*
1936 CALL cgeqrf( m, n, a, lda, work( itau ),
1937 $ work( iwork ), lwork-iwork+1, ierr )
1938 CALL clacpy( 'L', m, n, a, lda, u, ldu )
1939*
1940* Generate Q in U
1941* (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
1942* (RWorkspace: 0)
1943*
1944 CALL cungqr( m, m, n, u, ldu, work( itau ),
1945 $ work( iwork ), lwork-iwork+1, ierr )
1946*
1947* Copy R to WORK(IU), zeroing out below it
1948*
1949 CALL clacpy( 'U', n, n, a, lda, work( iu ),
1950 $ ldwrku )
1951 CALL claset( 'L', n-1, n-1, czero, czero,
1952 $ work( iu+1 ), ldwrku )
1953 ie = 1
1954 itauq = itau
1955 itaup = itauq + n
1956 iwork = itaup + n
1957*
1958* Bidiagonalize R in WORK(IU), copying result to VT
1959* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1960* (RWorkspace: need N)
1961*
1962 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1963 $ rwork( ie ), work( itauq ),
1964 $ work( itaup ), work( iwork ),
1965 $ lwork-iwork+1, ierr )
1966 CALL clacpy( 'U', n, n, work( iu ), ldwrku, vt,
1967 $ ldvt )
1968*
1969* Generate left bidiagonalizing vectors in WORK(IU)
1970* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1971* (RWorkspace: 0)
1972*
1973 CALL cungbr( 'Q', n, n, n, work( iu ), ldwrku,
1974 $ work( itauq ), work( iwork ),
1975 $ lwork-iwork+1, ierr )
1976*
1977* Generate right bidiagonalizing vectors in VT
1978* (CWorkspace: need N*N+3*N-1,
1979* prefer N*N+2*N+(N-1)*NB)
1980* (RWorkspace: need 0)
1981*
1982 CALL cungbr( 'P', n, n, n, vt, ldvt,
1983 $ work( itaup ),
1984 $ work( iwork ), lwork-iwork+1, ierr )
1985 irwork = ie + n
1986*
1987* Perform bidiagonal QR iteration, computing left
1988* singular vectors of R in WORK(IU) and computing
1989* right singular vectors of R in VT
1990* (CWorkspace: need N*N)
1991* (RWorkspace: need BDSPAC)
1992*
1993 CALL cbdsqr( 'U', n, n, n, 0, s, rwork( ie ),
1994 $ vt,
1995 $ ldvt, work( iu ), ldwrku, cdum, 1,
1996 $ rwork( irwork ), info )
1997*
1998* Multiply Q in U by left singular vectors of R in
1999* WORK(IU), storing result in A
2000* (CWorkspace: need N*N)
2001* (RWorkspace: 0)
2002*
2003 CALL cgemm( 'N', 'N', m, n, n, cone, u, ldu,
2004 $ work( iu ), ldwrku, czero, a, lda )
2005*
2006* Copy left singular vectors of A from A to U
2007*
2008 CALL clacpy( 'F', m, n, a, lda, u, ldu )
2009*
2010 ELSE
2011*
2012* Insufficient workspace for a fast algorithm
2013*
2014 itau = 1
2015 iwork = itau + n
2016*
2017* Compute A=Q*R, copying result to U
2018* (CWorkspace: need 2*N, prefer N+N*NB)
2019* (RWorkspace: 0)
2020*
2021 CALL cgeqrf( m, n, a, lda, work( itau ),
2022 $ work( iwork ), lwork-iwork+1, ierr )
2023 CALL clacpy( 'L', m, n, a, lda, u, ldu )
2024*
2025* Generate Q in U
2026* (CWorkspace: need N+M, prefer N+M*NB)
2027* (RWorkspace: 0)
2028*
2029 CALL cungqr( m, m, n, u, ldu, work( itau ),
2030 $ work( iwork ), lwork-iwork+1, ierr )
2031*
2032* Copy R from A to VT, zeroing out below it
2033*
2034 CALL clacpy( 'U', n, n, a, lda, vt, ldvt )
2035 IF( n.GT.1 )
2036 $ CALL claset( 'L', n-1, n-1, czero, czero,
2037 $ vt( 2, 1 ), ldvt )
2038 ie = 1
2039 itauq = itau
2040 itaup = itauq + n
2041 iwork = itaup + n
2042*
2043* Bidiagonalize R in VT
2044* (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
2045* (RWorkspace: need N)
2046*
2047 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
2048 $ work( itauq ), work( itaup ),
2049 $ work( iwork ), lwork-iwork+1, ierr )
2050*
2051* Multiply Q in U by left bidiagonalizing vectors
2052* in VT
2053* (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
2054* (RWorkspace: 0)
2055*
2056 CALL cunmbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
2057 $ work( itauq ), u, ldu, work( iwork ),
2058 $ lwork-iwork+1, ierr )
2059*
2060* Generate right bidiagonalizing vectors in VT
2061* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2062* (RWorkspace: 0)
2063*
2064 CALL cungbr( 'P', n, n, n, vt, ldvt,
2065 $ work( itaup ),
2066 $ work( iwork ), lwork-iwork+1, ierr )
2067 irwork = ie + n
2068*
2069* Perform bidiagonal QR iteration, computing left
2070* singular vectors of A in U and computing right
2071* singular vectors of A in VT
2072* (CWorkspace: 0)
2073* (RWorkspace: need BDSPAC)
2074*
2075 CALL cbdsqr( 'U', n, n, m, 0, s, rwork( ie ),
2076 $ vt,
2077 $ ldvt, u, ldu, cdum, 1,
2078 $ rwork( irwork ), info )
2079*
2080 END IF
2081*
2082 END IF
2083*
2084 END IF
2085*
2086 ELSE
2087*
2088* M .LT. MNTHR
2089*
2090* Path 10 (M at least N, but not much larger)
2091* Reduce to bidiagonal form without QR decomposition
2092*
2093 ie = 1
2094 itauq = 1
2095 itaup = itauq + n
2096 iwork = itaup + n
2097*
2098* Bidiagonalize A
2099* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
2100* (RWorkspace: need N)
2101*
2102 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2103 $ work( itaup ), work( iwork ), lwork-iwork+1,
2104 $ ierr )
2105 IF( wntuas ) THEN
2106*
2107* If left singular vectors desired in U, copy result to U
2108* and generate left bidiagonalizing vectors in U
2109* (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB)
2110* (RWorkspace: 0)
2111*
2112 CALL clacpy( 'L', m, n, a, lda, u, ldu )
2113 IF( wntus )
2114 $ ncu = n
2115 IF( wntua )
2116 $ ncu = m
2117 CALL cungbr( 'Q', m, ncu, n, u, ldu, work( itauq ),
2118 $ work( iwork ), lwork-iwork+1, ierr )
2119 END IF
2120 IF( wntvas ) THEN
2121*
2122* If right singular vectors desired in VT, copy result to
2123* VT and generate right bidiagonalizing vectors in VT
2124* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2125* (RWorkspace: 0)
2126*
2127 CALL clacpy( 'U', n, n, a, lda, vt, ldvt )
2128 CALL cungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
2129 $ work( iwork ), lwork-iwork+1, ierr )
2130 END IF
2131 IF( wntuo ) THEN
2132*
2133* If left singular vectors desired in A, generate left
2134* bidiagonalizing vectors in A
2135* (CWorkspace: need 3*N, prefer 2*N+N*NB)
2136* (RWorkspace: 0)
2137*
2138 CALL cungbr( 'Q', m, n, n, a, lda, work( itauq ),
2139 $ work( iwork ), lwork-iwork+1, ierr )
2140 END IF
2141 IF( wntvo ) THEN
2142*
2143* If right singular vectors desired in A, generate right
2144* bidiagonalizing vectors in A
2145* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2146* (RWorkspace: 0)
2147*
2148 CALL cungbr( 'P', n, n, n, a, lda, work( itaup ),
2149 $ work( iwork ), lwork-iwork+1, ierr )
2150 END IF
2151 irwork = ie + n
2152 IF( wntuas .OR. wntuo )
2153 $ nru = m
2154 IF( wntun )
2155 $ nru = 0
2156 IF( wntvas .OR. wntvo )
2157 $ ncvt = n
2158 IF( wntvn )
2159 $ ncvt = 0
2160 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) ) THEN
2161*
2162* Perform bidiagonal QR iteration, if desired, computing
2163* left singular vectors in U and computing right singular
2164* vectors in VT
2165* (CWorkspace: 0)
2166* (RWorkspace: need BDSPAC)
2167*
2168 CALL cbdsqr( 'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2169 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
2170 $ info )
2171 ELSE IF( ( .NOT.wntuo ) .AND. wntvo ) THEN
2172*
2173* Perform bidiagonal QR iteration, if desired, computing
2174* left singular vectors in U and computing right singular
2175* vectors in A
2176* (CWorkspace: 0)
2177* (RWorkspace: need BDSPAC)
2178*
2179 CALL cbdsqr( 'U', n, ncvt, nru, 0, s, rwork( ie ), a,
2180 $ lda, u, ldu, cdum, 1, rwork( irwork ),
2181 $ info )
2182 ELSE
2183*
2184* Perform bidiagonal QR iteration, if desired, computing
2185* left singular vectors in A and computing right singular
2186* vectors in VT
2187* (CWorkspace: 0)
2188* (RWorkspace: need BDSPAC)
2189*
2190 CALL cbdsqr( 'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2191 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
2192 $ info )
2193 END IF
2194*
2195 END IF
2196*
2197 ELSE
2198*
2199* A has more columns than rows. If A has sufficiently more
2200* columns than rows, first reduce using the LQ decomposition (if
2201* sufficient workspace available)
2202*
2203 IF( n.GE.mnthr ) THEN
2204*
2205 IF( wntvn ) THEN
2206*
2207* Path 1t(N much larger than M, JOBVT='N')
2208* No right singular vectors to be computed
2209*
2210 itau = 1
2211 iwork = itau + m
2212*
2213* Compute A=L*Q
2214* (CWorkspace: need 2*M, prefer M+M*NB)
2215* (RWorkspace: 0)
2216*
2217 CALL cgelqf( m, n, a, lda, work( itau ),
2218 $ work( iwork ),
2219 $ lwork-iwork+1, ierr )
2220*
2221* Zero out above L
2222*
2223 CALL claset( 'U', m-1, m-1, czero, czero, a( 1, 2 ),
2224 $ lda )
2225 ie = 1
2226 itauq = 1
2227 itaup = itauq + m
2228 iwork = itaup + m
2229*
2230* Bidiagonalize L in A
2231* (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2232* (RWorkspace: need M)
2233*
2234 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2235 $ work( itauq ),
2236 $ work( itaup ), work( iwork ), lwork-iwork+1,
2237 $ ierr )
2238 IF( wntuo .OR. wntuas ) THEN
2239*
2240* If left singular vectors desired, generate Q
2241* (CWorkspace: need 3*M, prefer 2*M+M*NB)
2242* (RWorkspace: 0)
2243*
2244 CALL cungbr( 'Q', m, m, m, a, lda, work( itauq ),
2245 $ work( iwork ), lwork-iwork+1, ierr )
2246 END IF
2247 irwork = ie + m
2248 nru = 0
2249 IF( wntuo .OR. wntuas )
2250 $ nru = m
2251*
2252* Perform bidiagonal QR iteration, computing left singular
2253* vectors of A in A if desired
2254* (CWorkspace: 0)
2255* (RWorkspace: need BDSPAC)
2256*
2257 CALL cbdsqr( 'U', m, 0, nru, 0, s, rwork( ie ), cdum,
2258 $ 1,
2259 $ a, lda, cdum, 1, rwork( irwork ), info )
2260*
2261* If left singular vectors desired in U, copy them there
2262*
2263 IF( wntuas )
2264 $ CALL clacpy( 'F', m, m, a, lda, u, ldu )
2265*
2266 ELSE IF( wntvo .AND. wntun ) THEN
2267*
2268* Path 2t(N much larger than M, JOBU='N', JOBVT='O')
2269* M right singular vectors to be overwritten on A and
2270* no left singular vectors to be computed
2271*
2272 IF( lwork.GE.m*m+3*m ) THEN
2273*
2274* Sufficient workspace for a fast algorithm
2275*
2276 ir = 1
2277 IF( lwork.GE.max( wrkbl, lda*n )+lda*m ) THEN
2278*
2279* WORK(IU) is LDA by N and WORK(IR) is LDA by M
2280*
2281 ldwrku = lda
2282 chunk = n
2283 ldwrkr = lda
2284 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m ) THEN
2285*
2286* WORK(IU) is LDA by N and WORK(IR) is M by M
2287*
2288 ldwrku = lda
2289 chunk = n
2290 ldwrkr = m
2291 ELSE
2292*
2293* WORK(IU) is M by CHUNK and WORK(IR) is M by M
2294*
2295 ldwrku = m
2296 chunk = ( lwork-m*m ) / m
2297 ldwrkr = m
2298 END IF
2299 itau = ir + ldwrkr*m
2300 iwork = itau + m
2301*
2302* Compute A=L*Q
2303* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2304* (RWorkspace: 0)
2305*
2306 CALL cgelqf( m, n, a, lda, work( itau ),
2307 $ work( iwork ), lwork-iwork+1, ierr )
2308*
2309* Copy L to WORK(IR) and zero out above it
2310*
2311 CALL clacpy( 'L', m, m, a, lda, work( ir ),
2312 $ ldwrkr )
2313 CALL claset( 'U', m-1, m-1, czero, czero,
2314 $ work( ir+ldwrkr ), ldwrkr )
2315*
2316* Generate Q in A
2317* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2318* (RWorkspace: 0)
2319*
2320 CALL cunglq( m, n, m, a, lda, work( itau ),
2321 $ work( iwork ), lwork-iwork+1, ierr )
2322 ie = 1
2323 itauq = itau
2324 itaup = itauq + m
2325 iwork = itaup + m
2326*
2327* Bidiagonalize L in WORK(IR)
2328* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2329* (RWorkspace: need M)
2330*
2331 CALL cgebrd( m, m, work( ir ), ldwrkr, s,
2332 $ rwork( ie ),
2333 $ work( itauq ), work( itaup ),
2334 $ work( iwork ), lwork-iwork+1, ierr )
2335*
2336* Generate right vectors bidiagonalizing L
2337* (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2338* (RWorkspace: 0)
2339*
2340 CALL cungbr( 'P', m, m, m, work( ir ), ldwrkr,
2341 $ work( itaup ), work( iwork ),
2342 $ lwork-iwork+1, ierr )
2343 irwork = ie + m
2344*
2345* Perform bidiagonal QR iteration, computing right
2346* singular vectors of L in WORK(IR)
2347* (CWorkspace: need M*M)
2348* (RWorkspace: need BDSPAC)
2349*
2350 CALL cbdsqr( 'U', m, m, 0, 0, s, rwork( ie ),
2351 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2352 $ rwork( irwork ), info )
2353 iu = itauq
2354*
2355* Multiply right singular vectors of L in WORK(IR) by Q
2356* in A, storing result in WORK(IU) and copying to A
2357* (CWorkspace: need M*M+M, prefer M*M+M*N)
2358* (RWorkspace: 0)
2359*
2360 DO 30 i = 1, n, chunk
2361 blk = min( n-i+1, chunk )
2362 CALL cgemm( 'N', 'N', m, blk, m, cone,
2363 $ work( ir ),
2364 $ ldwrkr, a( 1, i ), lda, czero,
2365 $ work( iu ), ldwrku )
2366 CALL clacpy( 'F', m, blk, work( iu ), ldwrku,
2367 $ a( 1, i ), lda )
2368 30 CONTINUE
2369*
2370 ELSE
2371*
2372* Insufficient workspace for a fast algorithm
2373*
2374 ie = 1
2375 itauq = 1
2376 itaup = itauq + m
2377 iwork = itaup + m
2378*
2379* Bidiagonalize A
2380* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
2381* (RWorkspace: need M)
2382*
2383 CALL cgebrd( m, n, a, lda, s, rwork( ie ),
2384 $ work( itauq ), work( itaup ),
2385 $ work( iwork ), lwork-iwork+1, ierr )
2386*
2387* Generate right vectors bidiagonalizing A
2388* (CWorkspace: need 3*M, prefer 2*M+M*NB)
2389* (RWorkspace: 0)
2390*
2391 CALL cungbr( 'P', m, n, m, a, lda, work( itaup ),
2392 $ work( iwork ), lwork-iwork+1, ierr )
2393 irwork = ie + m
2394*
2395* Perform bidiagonal QR iteration, computing right
2396* singular vectors of A in A
2397* (CWorkspace: 0)
2398* (RWorkspace: need BDSPAC)
2399*
2400 CALL cbdsqr( 'L', m, n, 0, 0, s, rwork( ie ), a,
2401 $ lda,
2402 $ cdum, 1, cdum, 1, rwork( irwork ), info )
2403*
2404 END IF
2405*
2406 ELSE IF( wntvo .AND. wntuas ) THEN
2407*
2408* Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
2409* M right singular vectors to be overwritten on A and
2410* M left singular vectors to be computed in U
2411*
2412 IF( lwork.GE.m*m+3*m ) THEN
2413*
2414* Sufficient workspace for a fast algorithm
2415*
2416 ir = 1
2417 IF( lwork.GE.max( wrkbl, lda*n )+lda*m ) THEN
2418*
2419* WORK(IU) is LDA by N and WORK(IR) is LDA by M
2420*
2421 ldwrku = lda
2422 chunk = n
2423 ldwrkr = lda
2424 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m ) THEN
2425*
2426* WORK(IU) is LDA by N and WORK(IR) is M by M
2427*
2428 ldwrku = lda
2429 chunk = n
2430 ldwrkr = m
2431 ELSE
2432*
2433* WORK(IU) is M by CHUNK and WORK(IR) is M by M
2434*
2435 ldwrku = m
2436 chunk = ( lwork-m*m ) / m
2437 ldwrkr = m
2438 END IF
2439 itau = ir + ldwrkr*m
2440 iwork = itau + m
2441*
2442* Compute A=L*Q
2443* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2444* (RWorkspace: 0)
2445*
2446 CALL cgelqf( m, n, a, lda, work( itau ),
2447 $ work( iwork ), lwork-iwork+1, ierr )
2448*
2449* Copy L to U, zeroing about above it
2450*
2451 CALL clacpy( 'L', m, m, a, lda, u, ldu )
2452 CALL claset( 'U', m-1, m-1, czero, czero, u( 1,
2453 $ 2 ),
2454 $ ldu )
2455*
2456* Generate Q in A
2457* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2458* (RWorkspace: 0)
2459*
2460 CALL cunglq( m, n, m, a, lda, work( itau ),
2461 $ work( iwork ), lwork-iwork+1, ierr )
2462 ie = 1
2463 itauq = itau
2464 itaup = itauq + m
2465 iwork = itaup + m
2466*
2467* Bidiagonalize L in U, copying result to WORK(IR)
2468* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2469* (RWorkspace: need M)
2470*
2471 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
2472 $ work( itauq ), work( itaup ),
2473 $ work( iwork ), lwork-iwork+1, ierr )
2474 CALL clacpy( 'U', m, m, u, ldu, work( ir ),
2475 $ ldwrkr )
2476*
2477* Generate right vectors bidiagonalizing L in WORK(IR)
2478* (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2479* (RWorkspace: 0)
2480*
2481 CALL cungbr( 'P', m, m, m, work( ir ), ldwrkr,
2482 $ work( itaup ), work( iwork ),
2483 $ lwork-iwork+1, ierr )
2484*
2485* Generate left vectors bidiagonalizing L in U
2486* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
2487* (RWorkspace: 0)
2488*
2489 CALL cungbr( 'Q', m, m, m, u, ldu, work( itauq ),
2490 $ work( iwork ), lwork-iwork+1, ierr )
2491 irwork = ie + m
2492*
2493* Perform bidiagonal QR iteration, computing left
2494* singular vectors of L in U, and computing right
2495* singular vectors of L in WORK(IR)
2496* (CWorkspace: need M*M)
2497* (RWorkspace: need BDSPAC)
2498*
2499 CALL cbdsqr( 'U', m, m, m, 0, s, rwork( ie ),
2500 $ work( ir ), ldwrkr, u, ldu, cdum, 1,
2501 $ rwork( irwork ), info )
2502 iu = itauq
2503*
2504* Multiply right singular vectors of L in WORK(IR) by Q
2505* in A, storing result in WORK(IU) and copying to A
2506* (CWorkspace: need M*M+M, prefer M*M+M*N))
2507* (RWorkspace: 0)
2508*
2509 DO 40 i = 1, n, chunk
2510 blk = min( n-i+1, chunk )
2511 CALL cgemm( 'N', 'N', m, blk, m, cone,
2512 $ work( ir ),
2513 $ ldwrkr, a( 1, i ), lda, czero,
2514 $ work( iu ), ldwrku )
2515 CALL clacpy( 'F', m, blk, work( iu ), ldwrku,
2516 $ a( 1, i ), lda )
2517 40 CONTINUE
2518*
2519 ELSE
2520*
2521* Insufficient workspace for a fast algorithm
2522*
2523 itau = 1
2524 iwork = itau + m
2525*
2526* Compute A=L*Q
2527* (CWorkspace: need 2*M, prefer M+M*NB)
2528* (RWorkspace: 0)
2529*
2530 CALL cgelqf( m, n, a, lda, work( itau ),
2531 $ work( iwork ), lwork-iwork+1, ierr )
2532*
2533* Copy L to U, zeroing out above it
2534*
2535 CALL clacpy( 'L', m, m, a, lda, u, ldu )
2536 CALL claset( 'U', m-1, m-1, czero, czero, u( 1,
2537 $ 2 ),
2538 $ ldu )
2539*
2540* Generate Q in A
2541* (CWorkspace: need 2*M, prefer M+M*NB)
2542* (RWorkspace: 0)
2543*
2544 CALL cunglq( m, n, m, a, lda, work( itau ),
2545 $ work( iwork ), lwork-iwork+1, ierr )
2546 ie = 1
2547 itauq = itau
2548 itaup = itauq + m
2549 iwork = itaup + m
2550*
2551* Bidiagonalize L in U
2552* (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2553* (RWorkspace: need M)
2554*
2555 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
2556 $ work( itauq ), work( itaup ),
2557 $ work( iwork ), lwork-iwork+1, ierr )
2558*
2559* Multiply right vectors bidiagonalizing L by Q in A
2560* (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2561* (RWorkspace: 0)
2562*
2563 CALL cunmbr( 'P', 'L', 'C', m, n, m, u, ldu,
2564 $ work( itaup ), a, lda, work( iwork ),
2565 $ lwork-iwork+1, ierr )
2566*
2567* Generate left vectors bidiagonalizing L in U
2568* (CWorkspace: need 3*M, prefer 2*M+M*NB)
2569* (RWorkspace: 0)
2570*
2571 CALL cungbr( 'Q', m, m, m, u, ldu, work( itauq ),
2572 $ work( iwork ), lwork-iwork+1, ierr )
2573 irwork = ie + m
2574*
2575* Perform bidiagonal QR iteration, computing left
2576* singular vectors of A in U and computing right
2577* singular vectors of A in A
2578* (CWorkspace: 0)
2579* (RWorkspace: need BDSPAC)
2580*
2581 CALL cbdsqr( 'U', m, n, m, 0, s, rwork( ie ), a,
2582 $ lda,
2583 $ u, ldu, cdum, 1, rwork( irwork ), info )
2584*
2585 END IF
2586*
2587 ELSE IF( wntvs ) THEN
2588*
2589 IF( wntun ) THEN
2590*
2591* Path 4t(N much larger than M, JOBU='N', JOBVT='S')
2592* M right singular vectors to be computed in VT and
2593* no left singular vectors to be computed
2594*
2595 IF( lwork.GE.m*m+3*m ) THEN
2596*
2597* Sufficient workspace for a fast algorithm
2598*
2599 ir = 1
2600 IF( lwork.GE.wrkbl+lda*m ) THEN
2601*
2602* WORK(IR) is LDA by M
2603*
2604 ldwrkr = lda
2605 ELSE
2606*
2607* WORK(IR) is M by M
2608*
2609 ldwrkr = m
2610 END IF
2611 itau = ir + ldwrkr*m
2612 iwork = itau + m
2613*
2614* Compute A=L*Q
2615* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2616* (RWorkspace: 0)
2617*
2618 CALL cgelqf( m, n, a, lda, work( itau ),
2619 $ work( iwork ), lwork-iwork+1, ierr )
2620*
2621* Copy L to WORK(IR), zeroing out above it
2622*
2623 CALL clacpy( 'L', m, m, a, lda, work( ir ),
2624 $ ldwrkr )
2625 CALL claset( 'U', m-1, m-1, czero, czero,
2626 $ work( ir+ldwrkr ), ldwrkr )
2627*
2628* Generate Q in A
2629* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2630* (RWorkspace: 0)
2631*
2632 CALL cunglq( m, n, m, a, lda, work( itau ),
2633 $ work( iwork ), lwork-iwork+1, ierr )
2634 ie = 1
2635 itauq = itau
2636 itaup = itauq + m
2637 iwork = itaup + m
2638*
2639* Bidiagonalize L in WORK(IR)
2640* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2641* (RWorkspace: need M)
2642*
2643 CALL cgebrd( m, m, work( ir ), ldwrkr, s,
2644 $ rwork( ie ), work( itauq ),
2645 $ work( itaup ), work( iwork ),
2646 $ lwork-iwork+1, ierr )
2647*
2648* Generate right vectors bidiagonalizing L in
2649* WORK(IR)
2650* (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
2651* (RWorkspace: 0)
2652*
2653 CALL cungbr( 'P', m, m, m, work( ir ), ldwrkr,
2654 $ work( itaup ), work( iwork ),
2655 $ lwork-iwork+1, ierr )
2656 irwork = ie + m
2657*
2658* Perform bidiagonal QR iteration, computing right
2659* singular vectors of L in WORK(IR)
2660* (CWorkspace: need M*M)
2661* (RWorkspace: need BDSPAC)
2662*
2663 CALL cbdsqr( 'U', m, m, 0, 0, s, rwork( ie ),
2664 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2665 $ rwork( irwork ), info )
2666*
2667* Multiply right singular vectors of L in WORK(IR) by
2668* Q in A, storing result in VT
2669* (CWorkspace: need M*M)
2670* (RWorkspace: 0)
2671*
2672 CALL cgemm( 'N', 'N', m, n, m, cone, work( ir ),
2673 $ ldwrkr, a, lda, czero, vt, ldvt )
2674*
2675 ELSE
2676*
2677* Insufficient workspace for a fast algorithm
2678*
2679 itau = 1
2680 iwork = itau + m
2681*
2682* Compute A=L*Q
2683* (CWorkspace: need 2*M, prefer M+M*NB)
2684* (RWorkspace: 0)
2685*
2686 CALL cgelqf( m, n, a, lda, work( itau ),
2687 $ work( iwork ), lwork-iwork+1, ierr )
2688*
2689* Copy result to VT
2690*
2691 CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
2692*
2693* Generate Q in VT
2694* (CWorkspace: need 2*M, prefer M+M*NB)
2695* (RWorkspace: 0)
2696*
2697 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
2698 $ work( iwork ), lwork-iwork+1, ierr )
2699 ie = 1
2700 itauq = itau
2701 itaup = itauq + m
2702 iwork = itaup + m
2703*
2704* Zero out above L in A
2705*
2706 CALL claset( 'U', m-1, m-1, czero, czero,
2707 $ a( 1, 2 ), lda )
2708*
2709* Bidiagonalize L in A
2710* (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2711* (RWorkspace: need M)
2712*
2713 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2714 $ work( itauq ), work( itaup ),
2715 $ work( iwork ), lwork-iwork+1, ierr )
2716*
2717* Multiply right vectors bidiagonalizing L by Q in VT
2718* (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2719* (RWorkspace: 0)
2720*
2721 CALL cunmbr( 'P', 'L', 'C', m, n, m, a, lda,
2722 $ work( itaup ), vt, ldvt,
2723 $ work( iwork ), lwork-iwork+1, ierr )
2724 irwork = ie + m
2725*
2726* Perform bidiagonal QR iteration, computing right
2727* singular vectors of A in VT
2728* (CWorkspace: 0)
2729* (RWorkspace: need BDSPAC)
2730*
2731 CALL cbdsqr( 'U', m, n, 0, 0, s, rwork( ie ),
2732 $ vt,
2733 $ ldvt, cdum, 1, cdum, 1,
2734 $ rwork( irwork ), info )
2735*
2736 END IF
2737*
2738 ELSE IF( wntuo ) THEN
2739*
2740* Path 5t(N much larger than M, JOBU='O', JOBVT='S')
2741* M right singular vectors to be computed in VT and
2742* M left singular vectors to be overwritten on A
2743*
2744 IF( lwork.GE.2*m*m+3*m ) THEN
2745*
2746* Sufficient workspace for a fast algorithm
2747*
2748 iu = 1
2749 IF( lwork.GE.wrkbl+2*lda*m ) THEN
2750*
2751* WORK(IU) is LDA by M and WORK(IR) is LDA by M
2752*
2753 ldwrku = lda
2754 ir = iu + ldwrku*m
2755 ldwrkr = lda
2756 ELSE IF( lwork.GE.wrkbl+( lda+m )*m ) THEN
2757*
2758* WORK(IU) is LDA by M and WORK(IR) is M by M
2759*
2760 ldwrku = lda
2761 ir = iu + ldwrku*m
2762 ldwrkr = m
2763 ELSE
2764*
2765* WORK(IU) is M by M and WORK(IR) is M by M
2766*
2767 ldwrku = m
2768 ir = iu + ldwrku*m
2769 ldwrkr = m
2770 END IF
2771 itau = ir + ldwrkr*m
2772 iwork = itau + m
2773*
2774* Compute A=L*Q
2775* (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
2776* (RWorkspace: 0)
2777*
2778 CALL cgelqf( m, n, a, lda, work( itau ),
2779 $ work( iwork ), lwork-iwork+1, ierr )
2780*
2781* Copy L to WORK(IU), zeroing out below it
2782*
2783 CALL clacpy( 'L', m, m, a, lda, work( iu ),
2784 $ ldwrku )
2785 CALL claset( 'U', m-1, m-1, czero, czero,
2786 $ work( iu+ldwrku ), ldwrku )
2787*
2788* Generate Q in A
2789* (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
2790* (RWorkspace: 0)
2791*
2792 CALL cunglq( m, n, m, a, lda, work( itau ),
2793 $ work( iwork ), lwork-iwork+1, ierr )
2794 ie = 1
2795 itauq = itau
2796 itaup = itauq + m
2797 iwork = itaup + m
2798*
2799* Bidiagonalize L in WORK(IU), copying result to
2800* WORK(IR)
2801* (CWorkspace: need 2*M*M+3*M,
2802* prefer 2*M*M+2*M+2*M*NB)
2803* (RWorkspace: need M)
2804*
2805 CALL cgebrd( m, m, work( iu ), ldwrku, s,
2806 $ rwork( ie ), work( itauq ),
2807 $ work( itaup ), work( iwork ),
2808 $ lwork-iwork+1, ierr )
2809 CALL clacpy( 'L', m, m, work( iu ), ldwrku,
2810 $ work( ir ), ldwrkr )
2811*
2812* Generate right bidiagonalizing vectors in WORK(IU)
2813* (CWorkspace: need 2*M*M+3*M-1,
2814* prefer 2*M*M+2*M+(M-1)*NB)
2815* (RWorkspace: 0)
2816*
2817 CALL cungbr( 'P', m, m, m, work( iu ), ldwrku,
2818 $ work( itaup ), work( iwork ),
2819 $ lwork-iwork+1, ierr )
2820*
2821* Generate left bidiagonalizing vectors in WORK(IR)
2822* (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
2823* (RWorkspace: 0)
2824*
2825 CALL cungbr( 'Q', m, m, m, work( ir ), ldwrkr,
2826 $ work( itauq ), work( iwork ),
2827 $ lwork-iwork+1, ierr )
2828 irwork = ie + m
2829*
2830* Perform bidiagonal QR iteration, computing left
2831* singular vectors of L in WORK(IR) and computing
2832* right singular vectors of L in WORK(IU)
2833* (CWorkspace: need 2*M*M)
2834* (RWorkspace: need BDSPAC)
2835*
2836 CALL cbdsqr( 'U', m, m, m, 0, s, rwork( ie ),
2837 $ work( iu ), ldwrku, work( ir ),
2838 $ ldwrkr, cdum, 1, rwork( irwork ),
2839 $ info )
2840*
2841* Multiply right singular vectors of L in WORK(IU) by
2842* Q in A, storing result in VT
2843* (CWorkspace: need M*M)
2844* (RWorkspace: 0)
2845*
2846 CALL cgemm( 'N', 'N', m, n, m, cone, work( iu ),
2847 $ ldwrku, a, lda, czero, vt, ldvt )
2848*
2849* Copy left singular vectors of L to A
2850* (CWorkspace: need M*M)
2851* (RWorkspace: 0)
2852*
2853 CALL clacpy( 'F', m, m, work( ir ), ldwrkr, a,
2854 $ lda )
2855*
2856 ELSE
2857*
2858* Insufficient workspace for a fast algorithm
2859*
2860 itau = 1
2861 iwork = itau + m
2862*
2863* Compute A=L*Q, copying result to VT
2864* (CWorkspace: need 2*M, prefer M+M*NB)
2865* (RWorkspace: 0)
2866*
2867 CALL cgelqf( m, n, a, lda, work( itau ),
2868 $ work( iwork ), lwork-iwork+1, ierr )
2869 CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
2870*
2871* Generate Q in VT
2872* (CWorkspace: need 2*M, prefer M+M*NB)
2873* (RWorkspace: 0)
2874*
2875 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
2876 $ work( iwork ), lwork-iwork+1, ierr )
2877 ie = 1
2878 itauq = itau
2879 itaup = itauq + m
2880 iwork = itaup + m
2881*
2882* Zero out above L in A
2883*
2884 CALL claset( 'U', m-1, m-1, czero, czero,
2885 $ a( 1, 2 ), lda )
2886*
2887* Bidiagonalize L in A
2888* (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2889* (RWorkspace: need M)
2890*
2891 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2892 $ work( itauq ), work( itaup ),
2893 $ work( iwork ), lwork-iwork+1, ierr )
2894*
2895* Multiply right vectors bidiagonalizing L by Q in VT
2896* (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2897* (RWorkspace: 0)
2898*
2899 CALL cunmbr( 'P', 'L', 'C', m, n, m, a, lda,
2900 $ work( itaup ), vt, ldvt,
2901 $ work( iwork ), lwork-iwork+1, ierr )
2902*
2903* Generate left bidiagonalizing vectors of L in A
2904* (CWorkspace: need 3*M, prefer 2*M+M*NB)
2905* (RWorkspace: 0)
2906*
2907 CALL cungbr( 'Q', m, m, m, a, lda,
2908 $ work( itauq ),
2909 $ work( iwork ), lwork-iwork+1, ierr )
2910 irwork = ie + m
2911*
2912* Perform bidiagonal QR iteration, computing left
2913* singular vectors of A in A and computing right
2914* singular vectors of A in VT
2915* (CWorkspace: 0)
2916* (RWorkspace: need BDSPAC)
2917*
2918 CALL cbdsqr( 'U', m, n, m, 0, s, rwork( ie ),
2919 $ vt,
2920 $ ldvt, a, lda, cdum, 1,
2921 $ rwork( irwork ), info )
2922*
2923 END IF
2924*
2925 ELSE IF( wntuas ) THEN
2926*
2927* Path 6t(N much larger than M, JOBU='S' or 'A',
2928* JOBVT='S')
2929* M right singular vectors to be computed in VT and
2930* M left singular vectors to be computed in U
2931*
2932 IF( lwork.GE.m*m+3*m ) THEN
2933*
2934* Sufficient workspace for a fast algorithm
2935*
2936 iu = 1
2937 IF( lwork.GE.wrkbl+lda*m ) THEN
2938*
2939* WORK(IU) is LDA by N
2940*
2941 ldwrku = lda
2942 ELSE
2943*
2944* WORK(IU) is LDA by M
2945*
2946 ldwrku = m
2947 END IF
2948 itau = iu + ldwrku*m
2949 iwork = itau + m
2950*
2951* Compute A=L*Q
2952* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2953* (RWorkspace: 0)
2954*
2955 CALL cgelqf( m, n, a, lda, work( itau ),
2956 $ work( iwork ), lwork-iwork+1, ierr )
2957*
2958* Copy L to WORK(IU), zeroing out above it
2959*
2960 CALL clacpy( 'L', m, m, a, lda, work( iu ),
2961 $ ldwrku )
2962 CALL claset( 'U', m-1, m-1, czero, czero,
2963 $ work( iu+ldwrku ), ldwrku )
2964*
2965* Generate Q in A
2966* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2967* (RWorkspace: 0)
2968*
2969 CALL cunglq( m, n, m, a, lda, work( itau ),
2970 $ work( iwork ), lwork-iwork+1, ierr )
2971 ie = 1
2972 itauq = itau
2973 itaup = itauq + m
2974 iwork = itaup + m
2975*
2976* Bidiagonalize L in WORK(IU), copying result to U
2977* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2978* (RWorkspace: need M)
2979*
2980 CALL cgebrd( m, m, work( iu ), ldwrku, s,
2981 $ rwork( ie ), work( itauq ),
2982 $ work( itaup ), work( iwork ),
2983 $ lwork-iwork+1, ierr )
2984 CALL clacpy( 'L', m, m, work( iu ), ldwrku, u,
2985 $ ldu )
2986*
2987* Generate right bidiagonalizing vectors in WORK(IU)
2988* (CWorkspace: need M*M+3*M-1,
2989* prefer M*M+2*M+(M-1)*NB)
2990* (RWorkspace: 0)
2991*
2992 CALL cungbr( 'P', m, m, m, work( iu ), ldwrku,
2993 $ work( itaup ), work( iwork ),
2994 $ lwork-iwork+1, ierr )
2995*
2996* Generate left bidiagonalizing vectors in U
2997* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
2998* (RWorkspace: 0)
2999*
3000 CALL cungbr( 'Q', m, m, m, u, ldu,
3001 $ work( itauq ),
3002 $ work( iwork ), lwork-iwork+1, ierr )
3003 irwork = ie + m
3004*
3005* Perform bidiagonal QR iteration, computing left
3006* singular vectors of L in U and computing right
3007* singular vectors of L in WORK(IU)
3008* (CWorkspace: need M*M)
3009* (RWorkspace: need BDSPAC)
3010*
3011 CALL cbdsqr( 'U', m, m, m, 0, s, rwork( ie ),
3012 $ work( iu ), ldwrku, u, ldu, cdum, 1,
3013 $ rwork( irwork ), info )
3014*
3015* Multiply right singular vectors of L in WORK(IU) by
3016* Q in A, storing result in VT
3017* (CWorkspace: need M*M)
3018* (RWorkspace: 0)
3019*
3020 CALL cgemm( 'N', 'N', m, n, m, cone, work( iu ),
3021 $ ldwrku, a, lda, czero, vt, ldvt )
3022*
3023 ELSE
3024*
3025* Insufficient workspace for a fast algorithm
3026*
3027 itau = 1
3028 iwork = itau + m
3029*
3030* Compute A=L*Q, copying result to VT
3031* (CWorkspace: need 2*M, prefer M+M*NB)
3032* (RWorkspace: 0)
3033*
3034 CALL cgelqf( m, n, a, lda, work( itau ),
3035 $ work( iwork ), lwork-iwork+1, ierr )
3036 CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3037*
3038* Generate Q in VT
3039* (CWorkspace: need 2*M, prefer M+M*NB)
3040* (RWorkspace: 0)
3041*
3042 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
3043 $ work( iwork ), lwork-iwork+1, ierr )
3044*
3045* Copy L to U, zeroing out above it
3046*
3047 CALL clacpy( 'L', m, m, a, lda, u, ldu )
3048 CALL claset( 'U', m-1, m-1, czero, czero,
3049 $ u( 1, 2 ), ldu )
3050 ie = 1
3051 itauq = itau
3052 itaup = itauq + m
3053 iwork = itaup + m
3054*
3055* Bidiagonalize L in U
3056* (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3057* (RWorkspace: need M)
3058*
3059 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
3060 $ work( itauq ), work( itaup ),
3061 $ work( iwork ), lwork-iwork+1, ierr )
3062*
3063* Multiply right bidiagonalizing vectors in U by Q
3064* in VT
3065* (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3066* (RWorkspace: 0)
3067*
3068 CALL cunmbr( 'P', 'L', 'C', m, n, m, u, ldu,
3069 $ work( itaup ), vt, ldvt,
3070 $ work( iwork ), lwork-iwork+1, ierr )
3071*
3072* Generate left bidiagonalizing vectors in U
3073* (CWorkspace: need 3*M, prefer 2*M+M*NB)
3074* (RWorkspace: 0)
3075*
3076 CALL cungbr( 'Q', m, m, m, u, ldu,
3077 $ work( itauq ),
3078 $ work( iwork ), lwork-iwork+1, ierr )
3079 irwork = ie + m
3080*
3081* Perform bidiagonal QR iteration, computing left
3082* singular vectors of A in U and computing right
3083* singular vectors of A in VT
3084* (CWorkspace: 0)
3085* (RWorkspace: need BDSPAC)
3086*
3087 CALL cbdsqr( 'U', m, n, m, 0, s, rwork( ie ),
3088 $ vt,
3089 $ ldvt, u, ldu, cdum, 1,
3090 $ rwork( irwork ), info )
3091*
3092 END IF
3093*
3094 END IF
3095*
3096 ELSE IF( wntva ) THEN
3097*
3098 IF( wntun ) THEN
3099*
3100* Path 7t(N much larger than M, JOBU='N', JOBVT='A')
3101* N right singular vectors to be computed in VT and
3102* no left singular vectors to be computed
3103*
3104 IF( lwork.GE.m*m+max( n+m, 3*m ) ) THEN
3105*
3106* Sufficient workspace for a fast algorithm
3107*
3108 ir = 1
3109 IF( lwork.GE.wrkbl+lda*m ) THEN
3110*
3111* WORK(IR) is LDA by M
3112*
3113 ldwrkr = lda
3114 ELSE
3115*
3116* WORK(IR) is M by M
3117*
3118 ldwrkr = m
3119 END IF
3120 itau = ir + ldwrkr*m
3121 iwork = itau + m
3122*
3123* Compute A=L*Q, copying result to VT
3124* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3125* (RWorkspace: 0)
3126*
3127 CALL cgelqf( m, n, a, lda, work( itau ),
3128 $ work( iwork ), lwork-iwork+1, ierr )
3129 CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3130*
3131* Copy L to WORK(IR), zeroing out above it
3132*
3133 CALL clacpy( 'L', m, m, a, lda, work( ir ),
3134 $ ldwrkr )
3135 CALL claset( 'U', m-1, m-1, czero, czero,
3136 $ work( ir+ldwrkr ), ldwrkr )
3137*
3138* Generate Q in VT
3139* (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
3140* (RWorkspace: 0)
3141*
3142 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3143 $ work( iwork ), lwork-iwork+1, ierr )
3144 ie = 1
3145 itauq = itau
3146 itaup = itauq + m
3147 iwork = itaup + m
3148*
3149* Bidiagonalize L in WORK(IR)
3150* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3151* (RWorkspace: need M)
3152*
3153 CALL cgebrd( m, m, work( ir ), ldwrkr, s,
3154 $ rwork( ie ), work( itauq ),
3155 $ work( itaup ), work( iwork ),
3156 $ lwork-iwork+1, ierr )
3157*
3158* Generate right bidiagonalizing vectors in WORK(IR)
3159* (CWorkspace: need M*M+3*M-1,
3160* prefer M*M+2*M+(M-1)*NB)
3161* (RWorkspace: 0)
3162*
3163 CALL cungbr( 'P', m, m, m, work( ir ), ldwrkr,
3164 $ work( itaup ), work( iwork ),
3165 $ lwork-iwork+1, ierr )
3166 irwork = ie + m
3167*
3168* Perform bidiagonal QR iteration, computing right
3169* singular vectors of L in WORK(IR)
3170* (CWorkspace: need M*M)
3171* (RWorkspace: need BDSPAC)
3172*
3173 CALL cbdsqr( 'U', m, m, 0, 0, s, rwork( ie ),
3174 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
3175 $ rwork( irwork ), info )
3176*
3177* Multiply right singular vectors of L in WORK(IR) by
3178* Q in VT, storing result in A
3179* (CWorkspace: need M*M)
3180* (RWorkspace: 0)
3181*
3182 CALL cgemm( 'N', 'N', m, n, m, cone, work( ir ),
3183 $ ldwrkr, vt, ldvt, czero, a, lda )
3184*
3185* Copy right singular vectors of A from A to VT
3186*
3187 CALL clacpy( 'F', m, n, a, lda, vt, ldvt )
3188*
3189 ELSE
3190*
3191* Insufficient workspace for a fast algorithm
3192*
3193 itau = 1
3194 iwork = itau + m
3195*
3196* Compute A=L*Q, copying result to VT
3197* (CWorkspace: need 2*M, prefer M+M*NB)
3198* (RWorkspace: 0)
3199*
3200 CALL cgelqf( m, n, a, lda, work( itau ),
3201 $ work( iwork ), lwork-iwork+1, ierr )
3202 CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3203*
3204* Generate Q in VT
3205* (CWorkspace: need M+N, prefer M+N*NB)
3206* (RWorkspace: 0)
3207*
3208 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3209 $ work( iwork ), lwork-iwork+1, ierr )
3210 ie = 1
3211 itauq = itau
3212 itaup = itauq + m
3213 iwork = itaup + m
3214*
3215* Zero out above L in A
3216*
3217 CALL claset( 'U', m-1, m-1, czero, czero,
3218 $ a( 1, 2 ), lda )
3219*
3220* Bidiagonalize L in A
3221* (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3222* (RWorkspace: need M)
3223*
3224 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
3225 $ work( itauq ), work( itaup ),
3226 $ work( iwork ), lwork-iwork+1, ierr )
3227*
3228* Multiply right bidiagonalizing vectors in A by Q
3229* in VT
3230* (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3231* (RWorkspace: 0)
3232*
3233 CALL cunmbr( 'P', 'L', 'C', m, n, m, a, lda,
3234 $ work( itaup ), vt, ldvt,
3235 $ work( iwork ), lwork-iwork+1, ierr )
3236 irwork = ie + m
3237*
3238* Perform bidiagonal QR iteration, computing right
3239* singular vectors of A in VT
3240* (CWorkspace: 0)
3241* (RWorkspace: need BDSPAC)
3242*
3243 CALL cbdsqr( 'U', m, n, 0, 0, s, rwork( ie ),
3244 $ vt,
3245 $ ldvt, cdum, 1, cdum, 1,
3246 $ rwork( irwork ), info )
3247*
3248 END IF
3249*
3250 ELSE IF( wntuo ) THEN
3251*
3252* Path 8t(N much larger than M, JOBU='O', JOBVT='A')
3253* N right singular vectors to be computed in VT and
3254* M left singular vectors to be overwritten on A
3255*
3256 IF( lwork.GE.2*m*m+max( n+m, 3*m ) ) THEN
3257*
3258* Sufficient workspace for a fast algorithm
3259*
3260 iu = 1
3261 IF( lwork.GE.wrkbl+2*lda*m ) THEN
3262*
3263* WORK(IU) is LDA by M and WORK(IR) is LDA by M
3264*
3265 ldwrku = lda
3266 ir = iu + ldwrku*m
3267 ldwrkr = lda
3268 ELSE IF( lwork.GE.wrkbl+( lda+m )*m ) THEN
3269*
3270* WORK(IU) is LDA by M and WORK(IR) is M by M
3271*
3272 ldwrku = lda
3273 ir = iu + ldwrku*m
3274 ldwrkr = m
3275 ELSE
3276*
3277* WORK(IU) is M by M and WORK(IR) is M by M
3278*
3279 ldwrku = m
3280 ir = iu + ldwrku*m
3281 ldwrkr = m
3282 END IF
3283 itau = ir + ldwrkr*m
3284 iwork = itau + m
3285*
3286* Compute A=L*Q, copying result to VT
3287* (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3288* (RWorkspace: 0)
3289*
3290 CALL cgelqf( m, n, a, lda, work( itau ),
3291 $ work( iwork ), lwork-iwork+1, ierr )
3292 CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3293*
3294* Generate Q in VT
3295* (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
3296* (RWorkspace: 0)
3297*
3298 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3299 $ work( iwork ), lwork-iwork+1, ierr )
3300*
3301* Copy L to WORK(IU), zeroing out above it
3302*
3303 CALL clacpy( 'L', m, m, a, lda, work( iu ),
3304 $ ldwrku )
3305 CALL claset( 'U', m-1, m-1, czero, czero,
3306 $ work( iu+ldwrku ), ldwrku )
3307 ie = 1
3308 itauq = itau
3309 itaup = itauq + m
3310 iwork = itaup + m
3311*
3312* Bidiagonalize L in WORK(IU), copying result to
3313* WORK(IR)
3314* (CWorkspace: need 2*M*M+3*M,
3315* prefer 2*M*M+2*M+2*M*NB)
3316* (RWorkspace: need M)
3317*
3318 CALL cgebrd( m, m, work( iu ), ldwrku, s,
3319 $ rwork( ie ), work( itauq ),
3320 $ work( itaup ), work( iwork ),
3321 $ lwork-iwork+1, ierr )
3322 CALL clacpy( 'L', m, m, work( iu ), ldwrku,
3323 $ work( ir ), ldwrkr )
3324*
3325* Generate right bidiagonalizing vectors in WORK(IU)
3326* (CWorkspace: need 2*M*M+3*M-1,
3327* prefer 2*M*M+2*M+(M-1)*NB)
3328* (RWorkspace: 0)
3329*
3330 CALL cungbr( 'P', m, m, m, work( iu ), ldwrku,
3331 $ work( itaup ), work( iwork ),
3332 $ lwork-iwork+1, ierr )
3333*
3334* Generate left bidiagonalizing vectors in WORK(IR)
3335* (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
3336* (RWorkspace: 0)
3337*
3338 CALL cungbr( 'Q', m, m, m, work( ir ), ldwrkr,
3339 $ work( itauq ), work( iwork ),
3340 $ lwork-iwork+1, ierr )
3341 irwork = ie + m
3342*
3343* Perform bidiagonal QR iteration, computing left
3344* singular vectors of L in WORK(IR) and computing
3345* right singular vectors of L in WORK(IU)
3346* (CWorkspace: need 2*M*M)
3347* (RWorkspace: need BDSPAC)
3348*
3349 CALL cbdsqr( 'U', m, m, m, 0, s, rwork( ie ),
3350 $ work( iu ), ldwrku, work( ir ),
3351 $ ldwrkr, cdum, 1, rwork( irwork ),
3352 $ info )
3353*
3354* Multiply right singular vectors of L in WORK(IU) by
3355* Q in VT, storing result in A
3356* (CWorkspace: need M*M)
3357* (RWorkspace: 0)
3358*
3359 CALL cgemm( 'N', 'N', m, n, m, cone, work( iu ),
3360 $ ldwrku, vt, ldvt, czero, a, lda )
3361*
3362* Copy right singular vectors of A from A to VT
3363*
3364 CALL clacpy( 'F', m, n, a, lda, vt, ldvt )
3365*
3366* Copy left singular vectors of A from WORK(IR) to A
3367*
3368 CALL clacpy( 'F', m, m, work( ir ), ldwrkr, a,
3369 $ lda )
3370*
3371 ELSE
3372*
3373* Insufficient workspace for a fast algorithm
3374*
3375 itau = 1
3376 iwork = itau + m
3377*
3378* Compute A=L*Q, copying result to VT
3379* (CWorkspace: need 2*M, prefer M+M*NB)
3380* (RWorkspace: 0)
3381*
3382 CALL cgelqf( m, n, a, lda, work( itau ),
3383 $ work( iwork ), lwork-iwork+1, ierr )
3384 CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3385*
3386* Generate Q in VT
3387* (CWorkspace: need M+N, prefer M+N*NB)
3388* (RWorkspace: 0)
3389*
3390 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3391 $ work( iwork ), lwork-iwork+1, ierr )
3392 ie = 1
3393 itauq = itau
3394 itaup = itauq + m
3395 iwork = itaup + m
3396*
3397* Zero out above L in A
3398*
3399 CALL claset( 'U', m-1, m-1, czero, czero,
3400 $ a( 1, 2 ), lda )
3401*
3402* Bidiagonalize L in A
3403* (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3404* (RWorkspace: need M)
3405*
3406 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
3407 $ work( itauq ), work( itaup ),
3408 $ work( iwork ), lwork-iwork+1, ierr )
3409*
3410* Multiply right bidiagonalizing vectors in A by Q
3411* in VT
3412* (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3413* (RWorkspace: 0)
3414*
3415 CALL cunmbr( 'P', 'L', 'C', m, n, m, a, lda,
3416 $ work( itaup ), vt, ldvt,
3417 $ work( iwork ), lwork-iwork+1, ierr )
3418*
3419* Generate left bidiagonalizing vectors in A
3420* (CWorkspace: need 3*M, prefer 2*M+M*NB)
3421* (RWorkspace: 0)
3422*
3423 CALL cungbr( 'Q', m, m, m, a, lda,
3424 $ work( itauq ),
3425 $ work( iwork ), lwork-iwork+1, ierr )
3426 irwork = ie + m
3427*
3428* Perform bidiagonal QR iteration, computing left
3429* singular vectors of A in A and computing right
3430* singular vectors of A in VT
3431* (CWorkspace: 0)
3432* (RWorkspace: need BDSPAC)
3433*
3434 CALL cbdsqr( 'U', m, n, m, 0, s, rwork( ie ),
3435 $ vt,
3436 $ ldvt, a, lda, cdum, 1,
3437 $ rwork( irwork ), info )
3438*
3439 END IF
3440*
3441 ELSE IF( wntuas ) THEN
3442*
3443* Path 9t(N much larger than M, JOBU='S' or 'A',
3444* JOBVT='A')
3445* N right singular vectors to be computed in VT and
3446* M left singular vectors to be computed in U
3447*
3448 IF( lwork.GE.m*m+max( n+m, 3*m ) ) THEN
3449*
3450* Sufficient workspace for a fast algorithm
3451*
3452 iu = 1
3453 IF( lwork.GE.wrkbl+lda*m ) THEN
3454*
3455* WORK(IU) is LDA by M
3456*
3457 ldwrku = lda
3458 ELSE
3459*
3460* WORK(IU) is M by M
3461*
3462 ldwrku = m
3463 END IF
3464 itau = iu + ldwrku*m
3465 iwork = itau + m
3466*
3467* Compute A=L*Q, copying result to VT
3468* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3469* (RWorkspace: 0)
3470*
3471 CALL cgelqf( m, n, a, lda, work( itau ),
3472 $ work( iwork ), lwork-iwork+1, ierr )
3473 CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3474*
3475* Generate Q in VT
3476* (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
3477* (RWorkspace: 0)
3478*
3479 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3480 $ work( iwork ), lwork-iwork+1, ierr )
3481*
3482* Copy L to WORK(IU), zeroing out above it
3483*
3484 CALL clacpy( 'L', m, m, a, lda, work( iu ),
3485 $ ldwrku )
3486 CALL claset( 'U', m-1, m-1, czero, czero,
3487 $ work( iu+ldwrku ), ldwrku )
3488 ie = 1
3489 itauq = itau
3490 itaup = itauq + m
3491 iwork = itaup + m
3492*
3493* Bidiagonalize L in WORK(IU), copying result to U
3494* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3495* (RWorkspace: need M)
3496*
3497 CALL cgebrd( m, m, work( iu ), ldwrku, s,
3498 $ rwork( ie ), work( itauq ),
3499 $ work( itaup ), work( iwork ),
3500 $ lwork-iwork+1, ierr )
3501 CALL clacpy( 'L', m, m, work( iu ), ldwrku, u,
3502 $ ldu )
3503*
3504* Generate right bidiagonalizing vectors in WORK(IU)
3505* (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
3506* (RWorkspace: 0)
3507*
3508 CALL cungbr( 'P', m, m, m, work( iu ), ldwrku,
3509 $ work( itaup ), work( iwork ),
3510 $ lwork-iwork+1, ierr )
3511*
3512* Generate left bidiagonalizing vectors in U
3513* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
3514* (RWorkspace: 0)
3515*
3516 CALL cungbr( 'Q', m, m, m, u, ldu,
3517 $ work( itauq ),
3518 $ work( iwork ), lwork-iwork+1, ierr )
3519 irwork = ie + m
3520*
3521* Perform bidiagonal QR iteration, computing left
3522* singular vectors of L in U and computing right
3523* singular vectors of L in WORK(IU)
3524* (CWorkspace: need M*M)
3525* (RWorkspace: need BDSPAC)
3526*
3527 CALL cbdsqr( 'U', m, m, m, 0, s, rwork( ie ),
3528 $ work( iu ), ldwrku, u, ldu, cdum, 1,
3529 $ rwork( irwork ), info )
3530*
3531* Multiply right singular vectors of L in WORK(IU) by
3532* Q in VT, storing result in A
3533* (CWorkspace: need M*M)
3534* (RWorkspace: 0)
3535*
3536 CALL cgemm( 'N', 'N', m, n, m, cone, work( iu ),
3537 $ ldwrku, vt, ldvt, czero, a, lda )
3538*
3539* Copy right singular vectors of A from A to VT
3540*
3541 CALL clacpy( 'F', m, n, a, lda, vt, ldvt )
3542*
3543 ELSE
3544*
3545* Insufficient workspace for a fast algorithm
3546*
3547 itau = 1
3548 iwork = itau + m
3549*
3550* Compute A=L*Q, copying result to VT
3551* (CWorkspace: need 2*M, prefer M+M*NB)
3552* (RWorkspace: 0)
3553*
3554 CALL cgelqf( m, n, a, lda, work( itau ),
3555 $ work( iwork ), lwork-iwork+1, ierr )
3556 CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3557*
3558* Generate Q in VT
3559* (CWorkspace: need M+N, prefer M+N*NB)
3560* (RWorkspace: 0)
3561*
3562 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3563 $ work( iwork ), lwork-iwork+1, ierr )
3564*
3565* Copy L to U, zeroing out above it
3566*
3567 CALL clacpy( 'L', m, m, a, lda, u, ldu )
3568 CALL claset( 'U', m-1, m-1, czero, czero,
3569 $ u( 1, 2 ), ldu )
3570 ie = 1
3571 itauq = itau
3572 itaup = itauq + m
3573 iwork = itaup + m
3574*
3575* Bidiagonalize L in U
3576* (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3577* (RWorkspace: need M)
3578*
3579 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
3580 $ work( itauq ), work( itaup ),
3581 $ work( iwork ), lwork-iwork+1, ierr )
3582*
3583* Multiply right bidiagonalizing vectors in U by Q
3584* in VT
3585* (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3586* (RWorkspace: 0)
3587*
3588 CALL cunmbr( 'P', 'L', 'C', m, n, m, u, ldu,
3589 $ work( itaup ), vt, ldvt,
3590 $ work( iwork ), lwork-iwork+1, ierr )
3591*
3592* Generate left bidiagonalizing vectors in U
3593* (CWorkspace: need 3*M, prefer 2*M+M*NB)
3594* (RWorkspace: 0)
3595*
3596 CALL cungbr( 'Q', m, m, m, u, ldu,
3597 $ work( itauq ),
3598 $ work( iwork ), lwork-iwork+1, ierr )
3599 irwork = ie + m
3600*
3601* Perform bidiagonal QR iteration, computing left
3602* singular vectors of A in U and computing right
3603* singular vectors of A in VT
3604* (CWorkspace: 0)
3605* (RWorkspace: need BDSPAC)
3606*
3607 CALL cbdsqr( 'U', m, n, m, 0, s, rwork( ie ),
3608 $ vt,
3609 $ ldvt, u, ldu, cdum, 1,
3610 $ rwork( irwork ), info )
3611*
3612 END IF
3613*
3614 END IF
3615*
3616 END IF
3617*
3618 ELSE
3619*
3620* N .LT. MNTHR
3621*
3622* Path 10t(N greater than M, but not much larger)
3623* Reduce to bidiagonal form without LQ decomposition
3624*
3625 ie = 1
3626 itauq = 1
3627 itaup = itauq + m
3628 iwork = itaup + m
3629*
3630* Bidiagonalize A
3631* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
3632* (RWorkspace: M)
3633*
3634 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
3635 $ work( itaup ), work( iwork ), lwork-iwork+1,
3636 $ ierr )
3637 IF( wntuas ) THEN
3638*
3639* If left singular vectors desired in U, copy result to U
3640* and generate left bidiagonalizing vectors in U
3641* (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
3642* (RWorkspace: 0)
3643*
3644 CALL clacpy( 'L', m, m, a, lda, u, ldu )
3645 CALL cungbr( 'Q', m, m, n, u, ldu, work( itauq ),
3646 $ work( iwork ), lwork-iwork+1, ierr )
3647 END IF
3648 IF( wntvas ) THEN
3649*
3650* If right singular vectors desired in VT, copy result to
3651* VT and generate right bidiagonalizing vectors in VT
3652* (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB)
3653* (RWorkspace: 0)
3654*
3655 CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3656 IF( wntva )
3657 $ nrvt = n
3658 IF( wntvs )
3659 $ nrvt = m
3660 CALL cungbr( 'P', nrvt, n, m, vt, ldvt, work( itaup ),
3661 $ work( iwork ), lwork-iwork+1, ierr )
3662 END IF
3663 IF( wntuo ) THEN
3664*
3665* If left singular vectors desired in A, generate left
3666* bidiagonalizing vectors in A
3667* (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
3668* (RWorkspace: 0)
3669*
3670 CALL cungbr( 'Q', m, m, n, a, lda, work( itauq ),
3671 $ work( iwork ), lwork-iwork+1, ierr )
3672 END IF
3673 IF( wntvo ) THEN
3674*
3675* If right singular vectors desired in A, generate right
3676* bidiagonalizing vectors in A
3677* (CWorkspace: need 3*M, prefer 2*M+M*NB)
3678* (RWorkspace: 0)
3679*
3680 CALL cungbr( 'P', m, n, m, a, lda, work( itaup ),
3681 $ work( iwork ), lwork-iwork+1, ierr )
3682 END IF
3683 irwork = ie + m
3684 IF( wntuas .OR. wntuo )
3685 $ nru = m
3686 IF( wntun )
3687 $ nru = 0
3688 IF( wntvas .OR. wntvo )
3689 $ ncvt = n
3690 IF( wntvn )
3691 $ ncvt = 0
3692 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) ) THEN
3693*
3694* Perform bidiagonal QR iteration, if desired, computing
3695* left singular vectors in U and computing right singular
3696* vectors in VT
3697* (CWorkspace: 0)
3698* (RWorkspace: need BDSPAC)
3699*
3700 CALL cbdsqr( 'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3701 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
3702 $ info )
3703 ELSE IF( ( .NOT.wntuo ) .AND. wntvo ) THEN
3704*
3705* Perform bidiagonal QR iteration, if desired, computing
3706* left singular vectors in U and computing right singular
3707* vectors in A
3708* (CWorkspace: 0)
3709* (RWorkspace: need BDSPAC)
3710*
3711 CALL cbdsqr( 'L', m, ncvt, nru, 0, s, rwork( ie ), a,
3712 $ lda, u, ldu, cdum, 1, rwork( irwork ),
3713 $ info )
3714 ELSE
3715*
3716* Perform bidiagonal QR iteration, if desired, computing
3717* left singular vectors in A and computing right singular
3718* vectors in VT
3719* (CWorkspace: 0)
3720* (RWorkspace: need BDSPAC)
3721*
3722 CALL cbdsqr( 'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3723 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
3724 $ info )
3725 END IF
3726*
3727 END IF
3728*
3729 END IF
3730*
3731* Undo scaling if necessary
3732*
3733 IF( iscl.EQ.1 ) THEN
3734 IF( anrm.GT.bignum )
3735 $ CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3736 $ ierr )
3737 IF( info.NE.0 .AND. anrm.GT.bignum )
3738 $ CALL slascl( 'G', 0, 0, bignum, anrm, minmn-1, 1,
3739 $ rwork( ie ), minmn, ierr )
3740 IF( anrm.LT.smlnum )
3741 $ CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3742 $ ierr )
3743 IF( info.NE.0 .AND. anrm.LT.smlnum )
3744 $ CALL slascl( 'G', 0, 0, smlnum, anrm, minmn-1, 1,
3745 $ rwork( ie ), minmn, ierr )
3746 END IF
3747*
3748* Return optimal workspace in WORK(1)
3749*
3750 work( 1 ) = sroundup_lwork(maxwrk)
3751*
3752 RETURN
3753*
3754* End of CGESVD
3755*
3756 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cbdsqr
subroutine cbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, rwork, info)
CBDSQR
Definition cbdsqr.f:233
cgebrd
subroutine cgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
CGEBRD
Definition cgebrd.f:205
cgelqf
subroutine cgelqf(m, n, a, lda, tau, work, lwork, info)
CGELQF
Definition cgelqf.f:142
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
cgeqrf
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
Definition cgeqrf.f:144
cgesvd
subroutine cgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, info)
CGESVD computes the singular value decomposition (SVD) for GE matrices
Definition cgesvd.f:213
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:101
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:142
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:142
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:104
cungbr
subroutine cungbr(vect, m, n, k, a, lda, tau, work, lwork, info)
CUNGBR
Definition cungbr.f:156
cunglq
subroutine cunglq(m, n, k, a, lda, tau, work, lwork, info)
CUNGLQ
Definition cunglq.f:125
cungqr
subroutine cungqr(m, n, k, a, lda, tau, work, lwork, info)
CUNGQR
Definition cungqr.f:126
cunmbr
subroutine cunmbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMBR
Definition cunmbr.f:195