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