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