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