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