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