1: *> \brief \b ZGESDD
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZGESDD + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesdd.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesdd.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesdd.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
22: * WORK, LWORK, RWORK, IWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER JOBZ
26: * INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
27: * ..
28: * .. Array Arguments ..
29: * INTEGER IWORK( * )
30: * DOUBLE PRECISION RWORK( * ), S( * )
31: * COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
32: * $ WORK( * )
33: * ..
34: *
35: *
36: *> \par Purpose:
37: * =============
38: *>
39: *> \verbatim
40: *>
41: *> ZGESDD computes the singular value decomposition (SVD) of a complex
42: *> M-by-N matrix A, optionally computing the left and/or right singular
43: *> vectors, by using divide-and-conquer method. The SVD is written
44: *>
45: *> A = U * SIGMA * conjugate-transpose(V)
46: *>
47: *> where SIGMA is an M-by-N matrix which is zero except for its
48: *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
49: *> V is an N-by-N unitary matrix. The diagonal elements of SIGMA
50: *> are the singular values of A; they are real and non-negative, and
51: *> are returned in descending order. The first min(m,n) columns of
52: *> U and V are the left and right singular vectors of A.
53: *>
54: *> Note that the routine returns VT = V**H, not V.
55: *>
56: *> The divide and conquer algorithm makes very mild assumptions about
57: *> floating point arithmetic. It will work on machines with a guard
58: *> digit in add/subtract, or on those binary machines without guard
59: *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
60: *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
61: *> without guard digits, but we know of none.
62: *> \endverbatim
63: *
64: * Arguments:
65: * ==========
66: *
67: *> \param[in] JOBZ
68: *> \verbatim
69: *> JOBZ is CHARACTER*1
70: *> Specifies options for computing all or part of the matrix U:
71: *> = 'A': all M columns of U and all N rows of V**H are
72: *> returned in the arrays U and VT;
73: *> = 'S': the first min(M,N) columns of U and the first
74: *> min(M,N) rows of V**H are returned in the arrays U
75: *> and VT;
76: *> = 'O': If M >= N, the first N columns of U are overwritten
77: *> in the array A and all rows of V**H are returned in
78: *> the array VT;
79: *> otherwise, all columns of U are returned in the
80: *> array U and the first M rows of V**H are overwritten
81: *> in the array A;
82: *> = 'N': no columns of U or rows of V**H are computed.
83: *> \endverbatim
84: *>
85: *> \param[in] M
86: *> \verbatim
87: *> M is INTEGER
88: *> The number of rows of the input matrix A. M >= 0.
89: *> \endverbatim
90: *>
91: *> \param[in] N
92: *> \verbatim
93: *> N is INTEGER
94: *> The number of columns of the input matrix A. N >= 0.
95: *> \endverbatim
96: *>
97: *> \param[in,out] A
98: *> \verbatim
99: *> A is COMPLEX*16 array, dimension (LDA,N)
100: *> On entry, the M-by-N matrix A.
101: *> On exit,
102: *> if JOBZ = 'O', A is overwritten with the first N columns
103: *> of U (the left singular vectors, stored
104: *> columnwise) if M >= N;
105: *> A is overwritten with the first M rows
106: *> of V**H (the right singular vectors, stored
107: *> rowwise) otherwise.
108: *> if JOBZ .ne. 'O', the contents of A are destroyed.
109: *> \endverbatim
110: *>
111: *> \param[in] LDA
112: *> \verbatim
113: *> LDA is INTEGER
114: *> The leading dimension of the array A. LDA >= max(1,M).
115: *> \endverbatim
116: *>
117: *> \param[out] S
118: *> \verbatim
119: *> S is DOUBLE PRECISION array, dimension (min(M,N))
120: *> The singular values of A, sorted so that S(i) >= S(i+1).
121: *> \endverbatim
122: *>
123: *> \param[out] U
124: *> \verbatim
125: *> U is COMPLEX*16 array, dimension (LDU,UCOL)
126: *> UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
127: *> UCOL = min(M,N) if JOBZ = 'S'.
128: *> If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
129: *> unitary matrix U;
130: *> if JOBZ = 'S', U contains the first min(M,N) columns of U
131: *> (the left singular vectors, stored columnwise);
132: *> if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
133: *> \endverbatim
134: *>
135: *> \param[in] LDU
136: *> \verbatim
137: *> LDU is INTEGER
138: *> The leading dimension of the array U. LDU >= 1;
139: *> if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
140: *> \endverbatim
141: *>
142: *> \param[out] VT
143: *> \verbatim
144: *> VT is COMPLEX*16 array, dimension (LDVT,N)
145: *> If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
146: *> N-by-N unitary matrix V**H;
147: *> if JOBZ = 'S', VT contains the first min(M,N) rows of
148: *> V**H (the right singular vectors, stored rowwise);
149: *> if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
150: *> \endverbatim
151: *>
152: *> \param[in] LDVT
153: *> \verbatim
154: *> LDVT is INTEGER
155: *> The leading dimension of the array VT. LDVT >= 1;
156: *> if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
157: *> if JOBZ = 'S', LDVT >= min(M,N).
158: *> \endverbatim
159: *>
160: *> \param[out] WORK
161: *> \verbatim
162: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
163: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
164: *> \endverbatim
165: *>
166: *> \param[in] LWORK
167: *> \verbatim
168: *> LWORK is INTEGER
169: *> The dimension of the array WORK. LWORK >= 1.
170: *> If LWORK = -1, a workspace query is assumed. The optimal
171: *> size for the WORK array is calculated and stored in WORK(1),
172: *> and no other work except argument checking is performed.
173: *>
174: *> Let mx = max(M,N) and mn = min(M,N).
175: *> If JOBZ = 'N', LWORK >= 2*mn + mx.
176: *> If JOBZ = 'O', LWORK >= 2*mn*mn + 2*mn + mx.
177: *> If JOBZ = 'S', LWORK >= mn*mn + 3*mn.
178: *> If JOBZ = 'A', LWORK >= mn*mn + 2*mn + mx.
179: *> These are not tight minimums in all cases; see comments inside code.
180: *> For good performance, LWORK should generally be larger;
181: *> a query is recommended.
182: *> \endverbatim
183: *>
184: *> \param[out] RWORK
185: *> \verbatim
186: *> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
187: *> Let mx = max(M,N) and mn = min(M,N).
188: *> If JOBZ = 'N', LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn);
189: *> else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn;
190: *> else LRWORK >= max( 5*mn*mn + 5*mn,
191: *> 2*mx*mn + 2*mn*mn + mn ).
192: *> \endverbatim
193: *>
194: *> \param[out] IWORK
195: *> \verbatim
196: *> IWORK is INTEGER array, dimension (8*min(M,N))
197: *> \endverbatim
198: *>
199: *> \param[out] INFO
200: *> \verbatim
201: *> INFO is INTEGER
202: *> = 0: successful exit.
203: *> < 0: if INFO = -i, the i-th argument had an illegal value.
204: *> > 0: The updating process of DBDSDC did not converge.
205: *> \endverbatim
206: *
207: * Authors:
208: * ========
209: *
210: *> \author Univ. of Tennessee
211: *> \author Univ. of California Berkeley
212: *> \author Univ. of Colorado Denver
213: *> \author NAG Ltd.
214: *
215: *> \date June 2016
216: *
217: *> \ingroup complex16GEsing
218: *
219: *> \par Contributors:
220: * ==================
221: *>
222: *> Ming Gu and Huan Ren, Computer Science Division, University of
223: *> California at Berkeley, USA
224: *>
225: * =====================================================================
226: SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
227: $ WORK, LWORK, RWORK, IWORK, INFO )
228: implicit none
229: *
230: * -- LAPACK driver routine (version 3.7.0) --
231: * -- LAPACK is a software package provided by Univ. of Tennessee, --
232: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
233: * June 2016
234: *
235: * .. Scalar Arguments ..
236: CHARACTER JOBZ
237: INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
238: * ..
239: * .. Array Arguments ..
240: INTEGER IWORK( * )
241: DOUBLE PRECISION RWORK( * ), S( * )
242: COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
243: $ WORK( * )
244: * ..
245: *
246: * =====================================================================
247: *
248: * .. Parameters ..
249: COMPLEX*16 CZERO, CONE
250: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
251: $ CONE = ( 1.0D+0, 0.0D+0 ) )
252: DOUBLE PRECISION ZERO, ONE
253: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
254: * ..
255: * .. Local Scalars ..
256: LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
257: INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
258: $ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
259: $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
260: $ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL
261: INTEGER LWORK_ZGEBRD_MN, LWORK_ZGEBRD_MM,
262: $ LWORK_ZGEBRD_NN, LWORK_ZGELQF_MN,
263: $ LWORK_ZGEQRF_MN,
264: $ LWORK_ZUNGBR_P_MN, LWORK_ZUNGBR_P_NN,
265: $ LWORK_ZUNGBR_Q_MN, LWORK_ZUNGBR_Q_MM,
266: $ LWORK_ZUNGLQ_MN, LWORK_ZUNGLQ_NN,
267: $ LWORK_ZUNGQR_MM, LWORK_ZUNGQR_MN,
268: $ LWORK_ZUNMBR_PRC_MM, LWORK_ZUNMBR_QLN_MM,
269: $ LWORK_ZUNMBR_PRC_MN, LWORK_ZUNMBR_QLN_MN,
270: $ LWORK_ZUNMBR_PRC_NN, LWORK_ZUNMBR_QLN_NN
271: DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
272: * ..
273: * .. Local Arrays ..
274: INTEGER IDUM( 1 )
275: DOUBLE PRECISION DUM( 1 )
276: COMPLEX*16 CDUM( 1 )
277: * ..
278: * .. External Subroutines ..
279: EXTERNAL DBDSDC, DLASCL, XERBLA, ZGEBRD, ZGELQF, ZGEMM,
280: $ ZGEQRF, ZLACP2, ZLACPY, ZLACRM, ZLARCM, ZLASCL,
281: $ ZLASET, ZUNGBR, ZUNGLQ, ZUNGQR, ZUNMBR
282: * ..
283: * .. External Functions ..
284: LOGICAL LSAME
285: DOUBLE PRECISION DLAMCH, ZLANGE
286: EXTERNAL LSAME, DLAMCH, ZLANGE
287: * ..
288: * .. Intrinsic Functions ..
289: INTRINSIC INT, MAX, MIN, SQRT
290: * ..
291: * .. Executable Statements ..
292: *
293: * Test the input arguments
294: *
295: INFO = 0
296: MINMN = MIN( M, N )
297: MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 )
298: MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 )
299: WNTQA = LSAME( JOBZ, 'A' )
300: WNTQS = LSAME( JOBZ, 'S' )
301: WNTQAS = WNTQA .OR. WNTQS
302: WNTQO = LSAME( JOBZ, 'O' )
303: WNTQN = LSAME( JOBZ, 'N' )
304: LQUERY = ( LWORK.EQ.-1 )
305: MINWRK = 1
306: MAXWRK = 1
307: *
308: IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
309: INFO = -1
310: ELSE IF( M.LT.0 ) THEN
311: INFO = -2
312: ELSE IF( N.LT.0 ) THEN
313: INFO = -3
314: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
315: INFO = -5
316: ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
317: $ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
318: INFO = -8
319: ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
320: $ ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
321: $ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
322: INFO = -10
323: END IF
324: *
325: * Compute workspace
326: * Note: Comments in the code beginning "Workspace:" describe the
327: * minimal amount of workspace allocated at that point in the code,
328: * as well as the preferred amount for good performance.
329: * CWorkspace refers to complex workspace, and RWorkspace to
330: * real workspace. NB refers to the optimal block size for the
331: * immediately following subroutine, as returned by ILAENV.)
332: *
333: IF( INFO.EQ.0 ) THEN
334: MINWRK = 1
335: MAXWRK = 1
336: IF( M.GE.N .AND. MINMN.GT.0 ) THEN
337: *
338: * There is no complex work space needed for bidiagonal SVD
339: * The real work space needed for bidiagonal SVD (dbdsdc) is
340: * BDSPAC = 3*N*N + 4*N for singular values and vectors;
341: * BDSPAC = 4*N for singular values only;
342: * not including e, RU, and RVT matrices.
343: *
344: * Compute space preferred for each routine
345: CALL ZGEBRD( M, N, CDUM(1), M, DUM(1), DUM(1), CDUM(1),
346: $ CDUM(1), CDUM(1), -1, IERR )
347: LWORK_ZGEBRD_MN = INT( CDUM(1) )
348: *
349: CALL ZGEBRD( N, N, CDUM(1), N, DUM(1), DUM(1), CDUM(1),
350: $ CDUM(1), CDUM(1), -1, IERR )
351: LWORK_ZGEBRD_NN = INT( CDUM(1) )
352: *
353: CALL ZGEQRF( M, N, CDUM(1), M, CDUM(1), CDUM(1), -1, IERR )
354: LWORK_ZGEQRF_MN = INT( CDUM(1) )
355: *
356: CALL ZUNGBR( 'P', N, N, N, CDUM(1), N, CDUM(1), CDUM(1),
357: $ -1, IERR )
358: LWORK_ZUNGBR_P_NN = INT( CDUM(1) )
359: *
360: CALL ZUNGBR( 'Q', M, M, N, CDUM(1), M, CDUM(1), CDUM(1),
361: $ -1, IERR )
362: LWORK_ZUNGBR_Q_MM = INT( CDUM(1) )
363: *
364: CALL ZUNGBR( 'Q', M, N, N, CDUM(1), M, CDUM(1), CDUM(1),
365: $ -1, IERR )
366: LWORK_ZUNGBR_Q_MN = INT( CDUM(1) )
367: *
368: CALL ZUNGQR( M, M, N, CDUM(1), M, CDUM(1), CDUM(1),
369: $ -1, IERR )
370: LWORK_ZUNGQR_MM = INT( CDUM(1) )
371: *
372: CALL ZUNGQR( M, N, N, CDUM(1), M, CDUM(1), CDUM(1),
373: $ -1, IERR )
374: LWORK_ZUNGQR_MN = INT( CDUM(1) )
375: *
376: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, CDUM(1), N, CDUM(1),
377: $ CDUM(1), N, CDUM(1), -1, IERR )
378: LWORK_ZUNMBR_PRC_NN = INT( CDUM(1) )
379: *
380: CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, CDUM(1), M, CDUM(1),
381: $ CDUM(1), M, CDUM(1), -1, IERR )
382: LWORK_ZUNMBR_QLN_MM = INT( CDUM(1) )
383: *
384: CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, CDUM(1), M, CDUM(1),
385: $ CDUM(1), M, CDUM(1), -1, IERR )
386: LWORK_ZUNMBR_QLN_MN = INT( CDUM(1) )
387: *
388: CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, CDUM(1), N, CDUM(1),
389: $ CDUM(1), N, CDUM(1), -1, IERR )
390: LWORK_ZUNMBR_QLN_NN = INT( CDUM(1) )
391: *
392: IF( M.GE.MNTHR1 ) THEN
393: IF( WNTQN ) THEN
394: *
395: * Path 1 (M >> N, JOBZ='N')
396: *
397: MAXWRK = N + LWORK_ZGEQRF_MN
398: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZGEBRD_NN )
399: MINWRK = 3*N
400: ELSE IF( WNTQO ) THEN
401: *
402: * Path 2 (M >> N, JOBZ='O')
403: *
404: WRKBL = N + LWORK_ZGEQRF_MN
405: WRKBL = MAX( WRKBL, N + LWORK_ZUNGQR_MN )
406: WRKBL = MAX( WRKBL, 2*N + LWORK_ZGEBRD_NN )
407: WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_QLN_NN )
408: WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_PRC_NN )
409: MAXWRK = M*N + N*N + WRKBL
410: MINWRK = 2*N*N + 3*N
411: ELSE IF( WNTQS ) THEN
412: *
413: * Path 3 (M >> N, JOBZ='S')
414: *
415: WRKBL = N + LWORK_ZGEQRF_MN
416: WRKBL = MAX( WRKBL, N + LWORK_ZUNGQR_MN )
417: WRKBL = MAX( WRKBL, 2*N + LWORK_ZGEBRD_NN )
418: WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_QLN_NN )
419: WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_PRC_NN )
420: MAXWRK = N*N + WRKBL
421: MINWRK = N*N + 3*N
422: ELSE IF( WNTQA ) THEN
423: *
424: * Path 4 (M >> N, JOBZ='A')
425: *
426: WRKBL = N + LWORK_ZGEQRF_MN
427: WRKBL = MAX( WRKBL, N + LWORK_ZUNGQR_MM )
428: WRKBL = MAX( WRKBL, 2*N + LWORK_ZGEBRD_NN )
429: WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_QLN_NN )
430: WRKBL = MAX( WRKBL, 2*N + LWORK_ZUNMBR_PRC_NN )
431: MAXWRK = N*N + WRKBL
432: MINWRK = N*N + MAX( 3*N, N + M )
433: END IF
434: ELSE IF( M.GE.MNTHR2 ) THEN
435: *
436: * Path 5 (M >> N, but not as much as MNTHR1)
437: *
438: MAXWRK = 2*N + LWORK_ZGEBRD_MN
439: MINWRK = 2*N + M
440: IF( WNTQO ) THEN
441: * Path 5o (M >> N, JOBZ='O')
442: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_P_NN )
443: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_Q_MN )
444: MAXWRK = MAXWRK + M*N
445: MINWRK = MINWRK + N*N
446: ELSE IF( WNTQS ) THEN
447: * Path 5s (M >> N, JOBZ='S')
448: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_P_NN )
449: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_Q_MN )
450: ELSE IF( WNTQA ) THEN
451: * Path 5a (M >> N, JOBZ='A')
452: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_P_NN )
453: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR_Q_MM )
454: END IF
455: ELSE
456: *
457: * Path 6 (M >= N, but not much larger)
458: *
459: MAXWRK = 2*N + LWORK_ZGEBRD_MN
460: MINWRK = 2*N + M
461: IF( WNTQO ) THEN
462: * Path 6o (M >= N, JOBZ='O')
463: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_PRC_NN )
464: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_QLN_MN )
465: MAXWRK = MAXWRK + M*N
466: MINWRK = MINWRK + N*N
467: ELSE IF( WNTQS ) THEN
468: * Path 6s (M >= N, JOBZ='S')
469: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_QLN_MN )
470: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_PRC_NN )
471: ELSE IF( WNTQA ) THEN
472: * Path 6a (M >= N, JOBZ='A')
473: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_QLN_MM )
474: MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR_PRC_NN )
475: END IF
476: END IF
477: ELSE IF( MINMN.GT.0 ) THEN
478: *
479: * There is no complex work space needed for bidiagonal SVD
480: * The real work space needed for bidiagonal SVD (dbdsdc) is
481: * BDSPAC = 3*M*M + 4*M for singular values and vectors;
482: * BDSPAC = 4*M for singular values only;
483: * not including e, RU, and RVT matrices.
484: *
485: * Compute space preferred for each routine
486: CALL ZGEBRD( M, N, CDUM(1), M, DUM(1), DUM(1), CDUM(1),
487: $ CDUM(1), CDUM(1), -1, IERR )
488: LWORK_ZGEBRD_MN = INT( CDUM(1) )
489: *
490: CALL ZGEBRD( M, M, CDUM(1), M, DUM(1), DUM(1), CDUM(1),
491: $ CDUM(1), CDUM(1), -1, IERR )
492: LWORK_ZGEBRD_MM = INT( CDUM(1) )
493: *
494: CALL ZGELQF( M, N, CDUM(1), M, CDUM(1), CDUM(1), -1, IERR )
495: LWORK_ZGELQF_MN = INT( CDUM(1) )
496: *
497: CALL ZUNGBR( 'P', M, N, M, CDUM(1), M, CDUM(1), CDUM(1),
498: $ -1, IERR )
499: LWORK_ZUNGBR_P_MN = INT( CDUM(1) )
500: *
501: CALL ZUNGBR( 'P', N, N, M, CDUM(1), N, CDUM(1), CDUM(1),
502: $ -1, IERR )
503: LWORK_ZUNGBR_P_NN = INT( CDUM(1) )
504: *
505: CALL ZUNGBR( 'Q', M, M, N, CDUM(1), M, CDUM(1), CDUM(1),
506: $ -1, IERR )
507: LWORK_ZUNGBR_Q_MM = INT( CDUM(1) )
508: *
509: CALL ZUNGLQ( M, N, M, CDUM(1), M, CDUM(1), CDUM(1),
510: $ -1, IERR )
511: LWORK_ZUNGLQ_MN = INT( CDUM(1) )
512: *
513: CALL ZUNGLQ( N, N, M, CDUM(1), N, CDUM(1), CDUM(1),
514: $ -1, IERR )
515: LWORK_ZUNGLQ_NN = INT( CDUM(1) )
516: *
517: CALL ZUNMBR( 'P', 'R', 'C', M, M, M, CDUM(1), M, CDUM(1),
518: $ CDUM(1), M, CDUM(1), -1, IERR )
519: LWORK_ZUNMBR_PRC_MM = INT( CDUM(1) )
520: *
521: CALL ZUNMBR( 'P', 'R', 'C', M, N, M, CDUM(1), M, CDUM(1),
522: $ CDUM(1), M, CDUM(1), -1, IERR )
523: LWORK_ZUNMBR_PRC_MN = INT( CDUM(1) )
524: *
525: CALL ZUNMBR( 'P', 'R', 'C', N, N, M, CDUM(1), N, CDUM(1),
526: $ CDUM(1), N, CDUM(1), -1, IERR )
527: LWORK_ZUNMBR_PRC_NN = INT( CDUM(1) )
528: *
529: CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, CDUM(1), M, CDUM(1),
530: $ CDUM(1), M, CDUM(1), -1, IERR )
531: LWORK_ZUNMBR_QLN_MM = INT( CDUM(1) )
532: *
533: IF( N.GE.MNTHR1 ) THEN
534: IF( WNTQN ) THEN
535: *
536: * Path 1t (N >> M, JOBZ='N')
537: *
538: MAXWRK = M + LWORK_ZGELQF_MN
539: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZGEBRD_MM )
540: MINWRK = 3*M
541: ELSE IF( WNTQO ) THEN
542: *
543: * Path 2t (N >> M, JOBZ='O')
544: *
545: WRKBL = M + LWORK_ZGELQF_MN
546: WRKBL = MAX( WRKBL, M + LWORK_ZUNGLQ_MN )
547: WRKBL = MAX( WRKBL, 2*M + LWORK_ZGEBRD_MM )
548: WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_QLN_MM )
549: WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_PRC_MM )
550: MAXWRK = M*N + M*M + WRKBL
551: MINWRK = 2*M*M + 3*M
552: ELSE IF( WNTQS ) THEN
553: *
554: * Path 3t (N >> M, JOBZ='S')
555: *
556: WRKBL = M + LWORK_ZGELQF_MN
557: WRKBL = MAX( WRKBL, M + LWORK_ZUNGLQ_MN )
558: WRKBL = MAX( WRKBL, 2*M + LWORK_ZGEBRD_MM )
559: WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_QLN_MM )
560: WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_PRC_MM )
561: MAXWRK = M*M + WRKBL
562: MINWRK = M*M + 3*M
563: ELSE IF( WNTQA ) THEN
564: *
565: * Path 4t (N >> M, JOBZ='A')
566: *
567: WRKBL = M + LWORK_ZGELQF_MN
568: WRKBL = MAX( WRKBL, M + LWORK_ZUNGLQ_NN )
569: WRKBL = MAX( WRKBL, 2*M + LWORK_ZGEBRD_MM )
570: WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_QLN_MM )
571: WRKBL = MAX( WRKBL, 2*M + LWORK_ZUNMBR_PRC_MM )
572: MAXWRK = M*M + WRKBL
573: MINWRK = M*M + MAX( 3*M, M + N )
574: END IF
575: ELSE IF( N.GE.MNTHR2 ) THEN
576: *
577: * Path 5t (N >> M, but not as much as MNTHR1)
578: *
579: MAXWRK = 2*M + LWORK_ZGEBRD_MN
580: MINWRK = 2*M + N
581: IF( WNTQO ) THEN
582: * Path 5to (N >> M, JOBZ='O')
583: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_Q_MM )
584: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_P_MN )
585: MAXWRK = MAXWRK + M*N
586: MINWRK = MINWRK + M*M
587: ELSE IF( WNTQS ) THEN
588: * Path 5ts (N >> M, JOBZ='S')
589: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_Q_MM )
590: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_P_MN )
591: ELSE IF( WNTQA ) THEN
592: * Path 5ta (N >> M, JOBZ='A')
593: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_Q_MM )
594: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR_P_NN )
595: END IF
596: ELSE
597: *
598: * Path 6t (N > M, but not much larger)
599: *
600: MAXWRK = 2*M + LWORK_ZGEBRD_MN
601: MINWRK = 2*M + N
602: IF( WNTQO ) THEN
603: * Path 6to (N > M, JOBZ='O')
604: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_QLN_MM )
605: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_PRC_MN )
606: MAXWRK = MAXWRK + M*N
607: MINWRK = MINWRK + M*M
608: ELSE IF( WNTQS ) THEN
609: * Path 6ts (N > M, JOBZ='S')
610: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_QLN_MM )
611: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_PRC_MN )
612: ELSE IF( WNTQA ) THEN
613: * Path 6ta (N > M, JOBZ='A')
614: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_QLN_MM )
615: MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR_PRC_NN )
616: END IF
617: END IF
618: END IF
619: MAXWRK = MAX( MAXWRK, MINWRK )
620: END IF
621: IF( INFO.EQ.0 ) THEN
622: WORK( 1 ) = MAXWRK
623: IF( LWORK.LT.MINWRK .AND. .NOT. LQUERY ) THEN
624: INFO = -12
625: END IF
626: END IF
627: *
628: IF( INFO.NE.0 ) THEN
629: CALL XERBLA( 'ZGESDD', -INFO )
630: RETURN
631: ELSE IF( LQUERY ) THEN
632: RETURN
633: END IF
634: *
635: * Quick return if possible
636: *
637: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
638: RETURN
639: END IF
640: *
641: * Get machine constants
642: *
643: EPS = DLAMCH( 'P' )
644: SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
645: BIGNUM = ONE / SMLNUM
646: *
647: * Scale A if max element outside range [SMLNUM,BIGNUM]
648: *
649: ANRM = ZLANGE( 'M', M, N, A, LDA, DUM )
650: ISCL = 0
651: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
652: ISCL = 1
653: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
654: ELSE IF( ANRM.GT.BIGNUM ) THEN
655: ISCL = 1
656: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
657: END IF
658: *
659: IF( M.GE.N ) THEN
660: *
661: * A has at least as many rows as columns. If A has sufficiently
662: * more rows than columns, first reduce using the QR
663: * decomposition (if sufficient workspace available)
664: *
665: IF( M.GE.MNTHR1 ) THEN
666: *
667: IF( WNTQN ) THEN
668: *
669: * Path 1 (M >> N, JOBZ='N')
670: * No singular vectors to be computed
671: *
672: ITAU = 1
673: NWORK = ITAU + N
674: *
675: * Compute A=Q*R
676: * CWorkspace: need N [tau] + N [work]
677: * CWorkspace: prefer N [tau] + N*NB [work]
678: * RWorkspace: need 0
679: *
680: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
681: $ LWORK-NWORK+1, IERR )
682: *
683: * Zero out below R
684: *
685: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
686: $ LDA )
687: IE = 1
688: ITAUQ = 1
689: ITAUP = ITAUQ + N
690: NWORK = ITAUP + N
691: *
692: * Bidiagonalize R in A
693: * CWorkspace: need 2*N [tauq, taup] + N [work]
694: * CWorkspace: prefer 2*N [tauq, taup] + 2*N*NB [work]
695: * RWorkspace: need N [e]
696: *
697: CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
698: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
699: $ IERR )
700: NRWORK = IE + N
701: *
702: * Perform bidiagonal SVD, compute singular values only
703: * CWorkspace: need 0
704: * RWorkspace: need N [e] + BDSPAC
705: *
706: CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM,1,DUM,1,
707: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
708: *
709: ELSE IF( WNTQO ) THEN
710: *
711: * Path 2 (M >> N, JOBZ='O')
712: * N left singular vectors to be overwritten on A and
713: * N right singular vectors to be computed in VT
714: *
715: IU = 1
716: *
717: * WORK(IU) is N by N
718: *
719: LDWRKU = N
720: IR = IU + LDWRKU*N
721: IF( LWORK .GE. M*N + N*N + 3*N ) THEN
722: *
723: * WORK(IR) is M by N
724: *
725: LDWRKR = M
726: ELSE
727: LDWRKR = ( LWORK - N*N - 3*N ) / N
728: END IF
729: ITAU = IR + LDWRKR*N
730: NWORK = ITAU + N
731: *
732: * Compute A=Q*R
733: * CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work]
734: * CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
735: * RWorkspace: need 0
736: *
737: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
738: $ LWORK-NWORK+1, IERR )
739: *
740: * Copy R to WORK( IR ), zeroing out below it
741: *
742: CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
743: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
744: $ LDWRKR )
745: *
746: * Generate Q in A
747: * CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work]
748: * CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
749: * RWorkspace: need 0
750: *
751: CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
752: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
753: IE = 1
754: ITAUQ = ITAU
755: ITAUP = ITAUQ + N
756: NWORK = ITAUP + N
757: *
758: * Bidiagonalize R in WORK(IR)
759: * CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
760: * CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
761: * RWorkspace: need N [e]
762: *
763: CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
764: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
765: $ LWORK-NWORK+1, IERR )
766: *
767: * Perform bidiagonal SVD, computing left singular vectors
768: * of R in WORK(IRU) and computing right singular vectors
769: * of R in WORK(IRVT)
770: * CWorkspace: need 0
771: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
772: *
773: IRU = IE + N
774: IRVT = IRU + N*N
775: NRWORK = IRVT + N*N
776: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
777: $ N, RWORK( IRVT ), N, DUM, IDUM,
778: $ RWORK( NRWORK ), IWORK, INFO )
779: *
780: * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
781: * Overwrite WORK(IU) by the left singular vectors of R
782: * CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
783: * CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
784: * RWorkspace: need 0
785: *
786: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
787: $ LDWRKU )
788: CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
789: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
790: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
791: *
792: * Copy real matrix RWORK(IRVT) to complex matrix VT
793: * Overwrite VT by the right singular vectors of R
794: * CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
795: * CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
796: * RWorkspace: need 0
797: *
798: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
799: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
800: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
801: $ LWORK-NWORK+1, IERR )
802: *
803: * Multiply Q in A by left singular vectors of R in
804: * WORK(IU), storing result in WORK(IR) and copying to A
805: * CWorkspace: need N*N [U] + N*N [R]
806: * CWorkspace: prefer N*N [U] + M*N [R]
807: * RWorkspace: need 0
808: *
809: DO 10 I = 1, M, LDWRKR
810: CHUNK = MIN( M-I+1, LDWRKR )
811: CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
812: $ LDA, WORK( IU ), LDWRKU, CZERO,
813: $ WORK( IR ), LDWRKR )
814: CALL ZLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
815: $ A( I, 1 ), LDA )
816: 10 CONTINUE
817: *
818: ELSE IF( WNTQS ) THEN
819: *
820: * Path 3 (M >> N, JOBZ='S')
821: * N left singular vectors to be computed in U and
822: * N right singular vectors to be computed in VT
823: *
824: IR = 1
825: *
826: * WORK(IR) is N by N
827: *
828: LDWRKR = N
829: ITAU = IR + LDWRKR*N
830: NWORK = ITAU + N
831: *
832: * Compute A=Q*R
833: * CWorkspace: need N*N [R] + N [tau] + N [work]
834: * CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
835: * RWorkspace: need 0
836: *
837: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
838: $ LWORK-NWORK+1, IERR )
839: *
840: * Copy R to WORK(IR), zeroing out below it
841: *
842: CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
843: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
844: $ LDWRKR )
845: *
846: * Generate Q in A
847: * CWorkspace: need N*N [R] + N [tau] + N [work]
848: * CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
849: * RWorkspace: need 0
850: *
851: CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
852: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
853: IE = 1
854: ITAUQ = ITAU
855: ITAUP = ITAUQ + N
856: NWORK = ITAUP + N
857: *
858: * Bidiagonalize R in WORK(IR)
859: * CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
860: * CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
861: * RWorkspace: need N [e]
862: *
863: CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
864: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
865: $ LWORK-NWORK+1, IERR )
866: *
867: * Perform bidiagonal SVD, computing left singular vectors
868: * of bidiagonal matrix in RWORK(IRU) and computing right
869: * singular vectors of bidiagonal matrix in RWORK(IRVT)
870: * CWorkspace: need 0
871: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
872: *
873: IRU = IE + N
874: IRVT = IRU + N*N
875: NRWORK = IRVT + N*N
876: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
877: $ N, RWORK( IRVT ), N, DUM, IDUM,
878: $ RWORK( NRWORK ), IWORK, INFO )
879: *
880: * Copy real matrix RWORK(IRU) to complex matrix U
881: * Overwrite U by left singular vectors of R
882: * CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
883: * CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
884: * RWorkspace: need 0
885: *
886: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
887: CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
888: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
889: $ LWORK-NWORK+1, IERR )
890: *
891: * Copy real matrix RWORK(IRVT) to complex matrix VT
892: * Overwrite VT by right singular vectors of R
893: * CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
894: * CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
895: * RWorkspace: need 0
896: *
897: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
898: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
899: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
900: $ LWORK-NWORK+1, IERR )
901: *
902: * Multiply Q in A by left singular vectors of R in
903: * WORK(IR), storing result in U
904: * CWorkspace: need N*N [R]
905: * RWorkspace: need 0
906: *
907: CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
908: CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),
909: $ LDWRKR, CZERO, U, LDU )
910: *
911: ELSE IF( WNTQA ) THEN
912: *
913: * Path 4 (M >> N, JOBZ='A')
914: * M left singular vectors to be computed in U and
915: * N right singular vectors to be computed in VT
916: *
917: IU = 1
918: *
919: * WORK(IU) is N by N
920: *
921: LDWRKU = N
922: ITAU = IU + LDWRKU*N
923: NWORK = ITAU + N
924: *
925: * Compute A=Q*R, copying result to U
926: * CWorkspace: need N*N [U] + N [tau] + N [work]
927: * CWorkspace: prefer N*N [U] + N [tau] + N*NB [work]
928: * RWorkspace: need 0
929: *
930: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
931: $ LWORK-NWORK+1, IERR )
932: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
933: *
934: * Generate Q in U
935: * CWorkspace: need N*N [U] + N [tau] + M [work]
936: * CWorkspace: prefer N*N [U] + N [tau] + M*NB [work]
937: * RWorkspace: need 0
938: *
939: CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
940: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
941: *
942: * Produce R in A, zeroing out below it
943: *
944: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
945: $ LDA )
946: IE = 1
947: ITAUQ = ITAU
948: ITAUP = ITAUQ + N
949: NWORK = ITAUP + N
950: *
951: * Bidiagonalize R in A
952: * CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
953: * CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + 2*N*NB [work]
954: * RWorkspace: need N [e]
955: *
956: CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
957: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
958: $ IERR )
959: IRU = IE + N
960: IRVT = IRU + N*N
961: NRWORK = IRVT + N*N
962: *
963: * Perform bidiagonal SVD, computing left singular vectors
964: * of bidiagonal matrix in RWORK(IRU) and computing right
965: * singular vectors of bidiagonal matrix in RWORK(IRVT)
966: * CWorkspace: need 0
967: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
968: *
969: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
970: $ N, RWORK( IRVT ), N, DUM, IDUM,
971: $ RWORK( NRWORK ), IWORK, INFO )
972: *
973: * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
974: * Overwrite WORK(IU) by left singular vectors of R
975: * CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
976: * CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
977: * RWorkspace: need 0
978: *
979: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
980: $ LDWRKU )
981: CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
982: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
983: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
984: *
985: * Copy real matrix RWORK(IRVT) to complex matrix VT
986: * Overwrite VT by right singular vectors of R
987: * CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
988: * CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
989: * RWorkspace: need 0
990: *
991: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
992: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
993: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
994: $ LWORK-NWORK+1, IERR )
995: *
996: * Multiply Q in U by left singular vectors of R in
997: * WORK(IU), storing result in A
998: * CWorkspace: need N*N [U]
999: * RWorkspace: need 0
1000: *
1001: CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),
1002: $ LDWRKU, CZERO, A, LDA )
1003: *
1004: * Copy left singular vectors of A from A to U
1005: *
1006: CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
1007: *
1008: END IF
1009: *
1010: ELSE IF( M.GE.MNTHR2 ) THEN
1011: *
1012: * MNTHR2 <= M < MNTHR1
1013: *
1014: * Path 5 (M >> N, but not as much as MNTHR1)
1015: * Reduce to bidiagonal form without QR decomposition, use
1016: * ZUNGBR and matrix multiplication to compute singular vectors
1017: *
1018: IE = 1
1019: NRWORK = IE + N
1020: ITAUQ = 1
1021: ITAUP = ITAUQ + N
1022: NWORK = ITAUP + N
1023: *
1024: * Bidiagonalize A
1025: * CWorkspace: need 2*N [tauq, taup] + M [work]
1026: * CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
1027: * RWorkspace: need N [e]
1028: *
1029: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1030: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1031: $ IERR )
1032: IF( WNTQN ) THEN
1033: *
1034: * Path 5n (M >> N, JOBZ='N')
1035: * Compute singular values only
1036: * CWorkspace: need 0
1037: * RWorkspace: need N [e] + BDSPAC
1038: *
1039: CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1,DUM,1,
1040: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1041: ELSE IF( WNTQO ) THEN
1042: IU = NWORK
1043: IRU = NRWORK
1044: IRVT = IRU + N*N
1045: NRWORK = IRVT + N*N
1046: *
1047: * Path 5o (M >> N, JOBZ='O')
1048: * Copy A to VT, generate P**H
1049: * CWorkspace: need 2*N [tauq, taup] + N [work]
1050: * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1051: * RWorkspace: need 0
1052: *
1053: CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
1054: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1055: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1056: *
1057: * Generate Q in A
1058: * CWorkspace: need 2*N [tauq, taup] + N [work]
1059: * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1060: * RWorkspace: need 0
1061: *
1062: CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
1063: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1064: *
1065: IF( LWORK .GE. M*N + 3*N ) THEN
1066: *
1067: * WORK( IU ) is M by N
1068: *
1069: LDWRKU = M
1070: ELSE
1071: *
1072: * WORK(IU) is LDWRKU by N
1073: *
1074: LDWRKU = ( LWORK - 3*N ) / N
1075: END IF
1076: NWORK = IU + LDWRKU*N
1077: *
1078: * Perform bidiagonal SVD, computing left singular vectors
1079: * of bidiagonal matrix in RWORK(IRU) and computing right
1080: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1081: * CWorkspace: need 0
1082: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1083: *
1084: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1085: $ N, RWORK( IRVT ), N, DUM, IDUM,
1086: $ RWORK( NRWORK ), IWORK, INFO )
1087: *
1088: * Multiply real matrix RWORK(IRVT) by P**H in VT,
1089: * storing the result in WORK(IU), copying to VT
1090: * CWorkspace: need 2*N [tauq, taup] + N*N [U]
1091: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1092: *
1093: CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,
1094: $ WORK( IU ), LDWRKU, RWORK( NRWORK ) )
1095: CALL ZLACPY( 'F', N, N, WORK( IU ), LDWRKU, VT, LDVT )
1096: *
1097: * Multiply Q in A by real matrix RWORK(IRU), storing the
1098: * result in WORK(IU), copying to A
1099: * CWorkspace: need 2*N [tauq, taup] + N*N [U]
1100: * CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
1101: * RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork]
1102: * RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1103: *
1104: NRWORK = IRVT
1105: DO 20 I = 1, M, LDWRKU
1106: CHUNK = MIN( M-I+1, LDWRKU )
1107: CALL ZLACRM( CHUNK, N, A( I, 1 ), LDA, RWORK( IRU ),
1108: $ N, WORK( IU ), LDWRKU, RWORK( NRWORK ) )
1109: CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
1110: $ A( I, 1 ), LDA )
1111: 20 CONTINUE
1112: *
1113: ELSE IF( WNTQS ) THEN
1114: *
1115: * Path 5s (M >> N, JOBZ='S')
1116: * Copy A to VT, generate P**H
1117: * CWorkspace: need 2*N [tauq, taup] + N [work]
1118: * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1119: * RWorkspace: need 0
1120: *
1121: CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
1122: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1123: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1124: *
1125: * Copy A to U, generate Q
1126: * CWorkspace: need 2*N [tauq, taup] + N [work]
1127: * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1128: * RWorkspace: need 0
1129: *
1130: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
1131: CALL ZUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),
1132: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1133: *
1134: * Perform bidiagonal SVD, computing left singular vectors
1135: * of bidiagonal matrix in RWORK(IRU) and computing right
1136: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1137: * CWorkspace: need 0
1138: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1139: *
1140: IRU = NRWORK
1141: IRVT = IRU + N*N
1142: NRWORK = IRVT + N*N
1143: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1144: $ N, RWORK( IRVT ), N, DUM, IDUM,
1145: $ RWORK( NRWORK ), IWORK, INFO )
1146: *
1147: * Multiply real matrix RWORK(IRVT) by P**H in VT,
1148: * storing the result in A, copying to VT
1149: * CWorkspace: need 0
1150: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1151: *
1152: CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
1153: $ RWORK( NRWORK ) )
1154: CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT )
1155: *
1156: * Multiply Q in U by real matrix RWORK(IRU), storing the
1157: * result in A, copying to U
1158: * CWorkspace: need 0
1159: * RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1160: *
1161: NRWORK = IRVT
1162: CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
1163: $ RWORK( NRWORK ) )
1164: CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
1165: ELSE
1166: *
1167: * Path 5a (M >> N, JOBZ='A')
1168: * Copy A to VT, generate P**H
1169: * CWorkspace: need 2*N [tauq, taup] + N [work]
1170: * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1171: * RWorkspace: need 0
1172: *
1173: CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
1174: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
1175: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1176: *
1177: * Copy A to U, generate Q
1178: * CWorkspace: need 2*N [tauq, taup] + M [work]
1179: * CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
1180: * RWorkspace: need 0
1181: *
1182: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
1183: CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1184: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1185: *
1186: * Perform bidiagonal SVD, computing left singular vectors
1187: * of bidiagonal matrix in RWORK(IRU) and computing right
1188: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1189: * CWorkspace: need 0
1190: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1191: *
1192: IRU = NRWORK
1193: IRVT = IRU + N*N
1194: NRWORK = IRVT + N*N
1195: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1196: $ N, RWORK( IRVT ), N, DUM, IDUM,
1197: $ RWORK( NRWORK ), IWORK, INFO )
1198: *
1199: * Multiply real matrix RWORK(IRVT) by P**H in VT,
1200: * storing the result in A, copying to VT
1201: * CWorkspace: need 0
1202: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1203: *
1204: CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
1205: $ RWORK( NRWORK ) )
1206: CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT )
1207: *
1208: * Multiply Q in U by real matrix RWORK(IRU), storing the
1209: * result in A, copying to U
1210: * CWorkspace: need 0
1211: * RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1212: *
1213: NRWORK = IRVT
1214: CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
1215: $ RWORK( NRWORK ) )
1216: CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
1217: END IF
1218: *
1219: ELSE
1220: *
1221: * M .LT. MNTHR2
1222: *
1223: * Path 6 (M >= N, but not much larger)
1224: * Reduce to bidiagonal form without QR decomposition
1225: * Use ZUNMBR to compute singular vectors
1226: *
1227: IE = 1
1228: NRWORK = IE + N
1229: ITAUQ = 1
1230: ITAUP = ITAUQ + N
1231: NWORK = ITAUP + N
1232: *
1233: * Bidiagonalize A
1234: * CWorkspace: need 2*N [tauq, taup] + M [work]
1235: * CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
1236: * RWorkspace: need N [e]
1237: *
1238: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1239: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1240: $ IERR )
1241: IF( WNTQN ) THEN
1242: *
1243: * Path 6n (M >= N, JOBZ='N')
1244: * Compute singular values only
1245: * CWorkspace: need 0
1246: * RWorkspace: need N [e] + BDSPAC
1247: *
1248: CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM,1,DUM,1,
1249: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1250: ELSE IF( WNTQO ) THEN
1251: IU = NWORK
1252: IRU = NRWORK
1253: IRVT = IRU + N*N
1254: NRWORK = IRVT + N*N
1255: IF( LWORK .GE. M*N + 3*N ) THEN
1256: *
1257: * WORK( IU ) is M by N
1258: *
1259: LDWRKU = M
1260: ELSE
1261: *
1262: * WORK( IU ) is LDWRKU by N
1263: *
1264: LDWRKU = ( LWORK - 3*N ) / N
1265: END IF
1266: NWORK = IU + LDWRKU*N
1267: *
1268: * Path 6o (M >= N, JOBZ='O')
1269: * Perform bidiagonal SVD, computing left singular vectors
1270: * of bidiagonal matrix in RWORK(IRU) and computing right
1271: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1272: * CWorkspace: need 0
1273: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1274: *
1275: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1276: $ N, RWORK( IRVT ), N, DUM, IDUM,
1277: $ RWORK( NRWORK ), IWORK, INFO )
1278: *
1279: * Copy real matrix RWORK(IRVT) to complex matrix VT
1280: * Overwrite VT by right singular vectors of A
1281: * CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work]
1282: * CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
1283: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1284: *
1285: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1286: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1287: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1288: $ LWORK-NWORK+1, IERR )
1289: *
1290: IF( LWORK .GE. M*N + 3*N ) THEN
1291: *
1292: * Path 6o-fast
1293: * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1294: * Overwrite WORK(IU) by left singular vectors of A, copying
1295: * to A
1296: * CWorkspace: need 2*N [tauq, taup] + M*N [U] + N [work]
1297: * CWorkspace: prefer 2*N [tauq, taup] + M*N [U] + N*NB [work]
1298: * RWorkspace: need N [e] + N*N [RU]
1299: *
1300: CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),
1301: $ LDWRKU )
1302: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
1303: $ LDWRKU )
1304: CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1305: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
1306: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1307: CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
1308: ELSE
1309: *
1310: * Path 6o-slow
1311: * Generate Q in A
1312: * CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work]
1313: * CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
1314: * RWorkspace: need 0
1315: *
1316: CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
1317: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1318: *
1319: * Multiply Q in A by real matrix RWORK(IRU), storing the
1320: * result in WORK(IU), copying to A
1321: * CWorkspace: need 2*N [tauq, taup] + N*N [U]
1322: * CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
1323: * RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork]
1324: * RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1325: *
1326: NRWORK = IRVT
1327: DO 30 I = 1, M, LDWRKU
1328: CHUNK = MIN( M-I+1, LDWRKU )
1329: CALL ZLACRM( CHUNK, N, A( I, 1 ), LDA,
1330: $ RWORK( IRU ), N, WORK( IU ), LDWRKU,
1331: $ RWORK( NRWORK ) )
1332: CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
1333: $ A( I, 1 ), LDA )
1334: 30 CONTINUE
1335: END IF
1336: *
1337: ELSE IF( WNTQS ) THEN
1338: *
1339: * Path 6s (M >= N, JOBZ='S')
1340: * Perform bidiagonal SVD, computing left singular vectors
1341: * of bidiagonal matrix in RWORK(IRU) and computing right
1342: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1343: * CWorkspace: need 0
1344: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1345: *
1346: IRU = NRWORK
1347: IRVT = IRU + N*N
1348: NRWORK = IRVT + N*N
1349: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1350: $ N, RWORK( IRVT ), N, DUM, IDUM,
1351: $ RWORK( NRWORK ), IWORK, INFO )
1352: *
1353: * Copy real matrix RWORK(IRU) to complex matrix U
1354: * Overwrite U by left singular vectors of A
1355: * CWorkspace: need 2*N [tauq, taup] + N [work]
1356: * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1357: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1358: *
1359: CALL ZLASET( 'F', M, N, CZERO, CZERO, U, LDU )
1360: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
1361: CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1362: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1363: $ LWORK-NWORK+1, IERR )
1364: *
1365: * Copy real matrix RWORK(IRVT) to complex matrix VT
1366: * Overwrite VT by right singular vectors of A
1367: * CWorkspace: need 2*N [tauq, taup] + N [work]
1368: * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1369: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1370: *
1371: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1372: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1373: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1374: $ LWORK-NWORK+1, IERR )
1375: ELSE
1376: *
1377: * Path 6a (M >= N, JOBZ='A')
1378: * Perform bidiagonal SVD, computing left singular vectors
1379: * of bidiagonal matrix in RWORK(IRU) and computing right
1380: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1381: * CWorkspace: need 0
1382: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1383: *
1384: IRU = NRWORK
1385: IRVT = IRU + N*N
1386: NRWORK = IRVT + N*N
1387: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1388: $ N, RWORK( IRVT ), N, DUM, IDUM,
1389: $ RWORK( NRWORK ), IWORK, INFO )
1390: *
1391: * Set the right corner of U to identity matrix
1392: *
1393: CALL ZLASET( 'F', M, M, CZERO, CZERO, U, LDU )
1394: IF( M.GT.N ) THEN
1395: CALL ZLASET( 'F', M-N, M-N, CZERO, CONE,
1396: $ U( N+1, N+1 ), LDU )
1397: END IF
1398: *
1399: * Copy real matrix RWORK(IRU) to complex matrix U
1400: * Overwrite U by left singular vectors of A
1401: * CWorkspace: need 2*N [tauq, taup] + M [work]
1402: * CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
1403: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1404: *
1405: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
1406: CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1407: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1408: $ LWORK-NWORK+1, IERR )
1409: *
1410: * Copy real matrix RWORK(IRVT) to complex matrix VT
1411: * Overwrite VT by right singular vectors of A
1412: * CWorkspace: need 2*N [tauq, taup] + N [work]
1413: * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1414: * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1415: *
1416: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1417: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1418: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1419: $ LWORK-NWORK+1, IERR )
1420: END IF
1421: *
1422: END IF
1423: *
1424: ELSE
1425: *
1426: * A has more columns than rows. If A has sufficiently more
1427: * columns than rows, first reduce using the LQ decomposition (if
1428: * sufficient workspace available)
1429: *
1430: IF( N.GE.MNTHR1 ) THEN
1431: *
1432: IF( WNTQN ) THEN
1433: *
1434: * Path 1t (N >> M, JOBZ='N')
1435: * No singular vectors to be computed
1436: *
1437: ITAU = 1
1438: NWORK = ITAU + M
1439: *
1440: * Compute A=L*Q
1441: * CWorkspace: need M [tau] + M [work]
1442: * CWorkspace: prefer M [tau] + M*NB [work]
1443: * RWorkspace: need 0
1444: *
1445: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1446: $ LWORK-NWORK+1, IERR )
1447: *
1448: * Zero out above L
1449: *
1450: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
1451: $ LDA )
1452: IE = 1
1453: ITAUQ = 1
1454: ITAUP = ITAUQ + M
1455: NWORK = ITAUP + M
1456: *
1457: * Bidiagonalize L in A
1458: * CWorkspace: need 2*M [tauq, taup] + M [work]
1459: * CWorkspace: prefer 2*M [tauq, taup] + 2*M*NB [work]
1460: * RWorkspace: need M [e]
1461: *
1462: CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1463: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1464: $ IERR )
1465: NRWORK = IE + M
1466: *
1467: * Perform bidiagonal SVD, compute singular values only
1468: * CWorkspace: need 0
1469: * RWorkspace: need M [e] + BDSPAC
1470: *
1471: CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
1472: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1473: *
1474: ELSE IF( WNTQO ) THEN
1475: *
1476: * Path 2t (N >> M, JOBZ='O')
1477: * M right singular vectors to be overwritten on A and
1478: * M left singular vectors to be computed in U
1479: *
1480: IVT = 1
1481: LDWKVT = M
1482: *
1483: * WORK(IVT) is M by M
1484: *
1485: IL = IVT + LDWKVT*M
1486: IF( LWORK .GE. M*N + M*M + 3*M ) THEN
1487: *
1488: * WORK(IL) M by N
1489: *
1490: LDWRKL = M
1491: CHUNK = N
1492: ELSE
1493: *
1494: * WORK(IL) is M by CHUNK
1495: *
1496: LDWRKL = M
1497: CHUNK = ( LWORK - M*M - 3*M ) / M
1498: END IF
1499: ITAU = IL + LDWRKL*CHUNK
1500: NWORK = ITAU + M
1501: *
1502: * Compute A=L*Q
1503: * CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1504: * CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1505: * RWorkspace: need 0
1506: *
1507: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1508: $ LWORK-NWORK+1, IERR )
1509: *
1510: * Copy L to WORK(IL), zeroing about above it
1511: *
1512: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1513: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
1514: $ WORK( IL+LDWRKL ), LDWRKL )
1515: *
1516: * Generate Q in A
1517: * CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1518: * CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1519: * RWorkspace: need 0
1520: *
1521: CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
1522: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1523: IE = 1
1524: ITAUQ = ITAU
1525: ITAUP = ITAUQ + M
1526: NWORK = ITAUP + M
1527: *
1528: * Bidiagonalize L in WORK(IL)
1529: * CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1530: * CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
1531: * RWorkspace: need M [e]
1532: *
1533: CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
1534: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1535: $ LWORK-NWORK+1, IERR )
1536: *
1537: * Perform bidiagonal SVD, computing left singular vectors
1538: * of bidiagonal matrix in RWORK(IRU) and computing right
1539: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1540: * CWorkspace: need 0
1541: * RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1542: *
1543: IRU = IE + M
1544: IRVT = IRU + M*M
1545: NRWORK = IRVT + M*M
1546: CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1547: $ M, RWORK( IRVT ), M, DUM, IDUM,
1548: $ RWORK( NRWORK ), IWORK, INFO )
1549: *
1550: * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1551: * Overwrite WORK(IU) by the left singular vectors of L
1552: * CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1553: * CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
1554: * RWorkspace: need 0
1555: *
1556: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1557: CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1558: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1559: $ LWORK-NWORK+1, IERR )
1560: *
1561: * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1562: * Overwrite WORK(IVT) by the right singular vectors of L
1563: * CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1564: * CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
1565: * RWorkspace: need 0
1566: *
1567: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
1568: $ LDWKVT )
1569: CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
1570: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1571: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1572: *
1573: * Multiply right singular vectors of L in WORK(IL) by Q
1574: * in A, storing result in WORK(IL) and copying to A
1575: * CWorkspace: need M*M [VT] + M*M [L]
1576: * CWorkspace: prefer M*M [VT] + M*N [L]
1577: * RWorkspace: need 0
1578: *
1579: DO 40 I = 1, N, CHUNK
1580: BLK = MIN( N-I+1, CHUNK )
1581: CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IVT ), M,
1582: $ A( 1, I ), LDA, CZERO, WORK( IL ),
1583: $ LDWRKL )
1584: CALL ZLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
1585: $ A( 1, I ), LDA )
1586: 40 CONTINUE
1587: *
1588: ELSE IF( WNTQS ) THEN
1589: *
1590: * Path 3t (N >> M, JOBZ='S')
1591: * M right singular vectors to be computed in VT and
1592: * M left singular vectors to be computed in U
1593: *
1594: IL = 1
1595: *
1596: * WORK(IL) is M by M
1597: *
1598: LDWRKL = M
1599: ITAU = IL + LDWRKL*M
1600: NWORK = ITAU + M
1601: *
1602: * Compute A=L*Q
1603: * CWorkspace: need M*M [L] + M [tau] + M [work]
1604: * CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
1605: * RWorkspace: need 0
1606: *
1607: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1608: $ LWORK-NWORK+1, IERR )
1609: *
1610: * Copy L to WORK(IL), zeroing out above it
1611: *
1612: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1613: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
1614: $ WORK( IL+LDWRKL ), LDWRKL )
1615: *
1616: * Generate Q in A
1617: * CWorkspace: need M*M [L] + M [tau] + M [work]
1618: * CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
1619: * RWorkspace: need 0
1620: *
1621: CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
1622: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1623: IE = 1
1624: ITAUQ = ITAU
1625: ITAUP = ITAUQ + M
1626: NWORK = ITAUP + M
1627: *
1628: * Bidiagonalize L in WORK(IL)
1629: * CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1630: * CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
1631: * RWorkspace: need M [e]
1632: *
1633: CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
1634: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1635: $ LWORK-NWORK+1, IERR )
1636: *
1637: * Perform bidiagonal SVD, computing left singular vectors
1638: * of bidiagonal matrix in RWORK(IRU) and computing right
1639: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1640: * CWorkspace: need 0
1641: * RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1642: *
1643: IRU = IE + M
1644: IRVT = IRU + M*M
1645: NRWORK = IRVT + M*M
1646: CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1647: $ M, RWORK( IRVT ), M, DUM, IDUM,
1648: $ RWORK( NRWORK ), IWORK, INFO )
1649: *
1650: * Copy real matrix RWORK(IRU) to complex matrix U
1651: * Overwrite U by left singular vectors of L
1652: * CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1653: * CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
1654: * RWorkspace: need 0
1655: *
1656: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1657: CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1658: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1659: $ LWORK-NWORK+1, IERR )
1660: *
1661: * Copy real matrix RWORK(IRVT) to complex matrix VT
1662: * Overwrite VT by left singular vectors of L
1663: * CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1664: * CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
1665: * RWorkspace: need 0
1666: *
1667: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
1668: CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
1669: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1670: $ LWORK-NWORK+1, IERR )
1671: *
1672: * Copy VT to WORK(IL), multiply right singular vectors of L
1673: * in WORK(IL) by Q in A, storing result in VT
1674: * CWorkspace: need M*M [L]
1675: * RWorkspace: need 0
1676: *
1677: CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1678: CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,
1679: $ A, LDA, CZERO, VT, LDVT )
1680: *
1681: ELSE IF( WNTQA ) THEN
1682: *
1683: * Path 4t (N >> M, JOBZ='A')
1684: * N right singular vectors to be computed in VT and
1685: * M left singular vectors to be computed in U
1686: *
1687: IVT = 1
1688: *
1689: * WORK(IVT) is M by M
1690: *
1691: LDWKVT = M
1692: ITAU = IVT + LDWKVT*M
1693: NWORK = ITAU + M
1694: *
1695: * Compute A=L*Q, copying result to VT
1696: * CWorkspace: need M*M [VT] + M [tau] + M [work]
1697: * CWorkspace: prefer M*M [VT] + M [tau] + M*NB [work]
1698: * RWorkspace: need 0
1699: *
1700: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1701: $ LWORK-NWORK+1, IERR )
1702: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
1703: *
1704: * Generate Q in VT
1705: * CWorkspace: need M*M [VT] + M [tau] + N [work]
1706: * CWorkspace: prefer M*M [VT] + M [tau] + N*NB [work]
1707: * RWorkspace: need 0
1708: *
1709: CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1710: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1711: *
1712: * Produce L in A, zeroing out above it
1713: *
1714: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
1715: $ LDA )
1716: IE = 1
1717: ITAUQ = ITAU
1718: ITAUP = ITAUQ + M
1719: NWORK = ITAUP + M
1720: *
1721: * Bidiagonalize L in A
1722: * CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1723: * CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + 2*M*NB [work]
1724: * RWorkspace: need M [e]
1725: *
1726: CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1727: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1728: $ IERR )
1729: *
1730: * Perform bidiagonal SVD, computing left singular vectors
1731: * of bidiagonal matrix in RWORK(IRU) and computing right
1732: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1733: * CWorkspace: need 0
1734: * RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1735: *
1736: IRU = IE + M
1737: IRVT = IRU + M*M
1738: NRWORK = IRVT + M*M
1739: CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1740: $ M, RWORK( IRVT ), M, DUM, IDUM,
1741: $ RWORK( NRWORK ), IWORK, INFO )
1742: *
1743: * Copy real matrix RWORK(IRU) to complex matrix U
1744: * Overwrite U by left singular vectors of L
1745: * CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1746: * CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
1747: * RWorkspace: need 0
1748: *
1749: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1750: CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1751: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1752: $ LWORK-NWORK+1, IERR )
1753: *
1754: * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1755: * Overwrite WORK(IVT) by right singular vectors of L
1756: * CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1757: * CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
1758: * RWorkspace: need 0
1759: *
1760: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
1761: $ LDWKVT )
1762: CALL ZUNMBR( 'P', 'R', 'C', M, M, M, A, LDA,
1763: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1764: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1765: *
1766: * Multiply right singular vectors of L in WORK(IVT) by
1767: * Q in VT, storing result in A
1768: * CWorkspace: need M*M [VT]
1769: * RWorkspace: need 0
1770: *
1771: CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT,
1772: $ VT, LDVT, CZERO, A, LDA )
1773: *
1774: * Copy right singular vectors of A from A to VT
1775: *
1776: CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
1777: *
1778: END IF
1779: *
1780: ELSE IF( N.GE.MNTHR2 ) THEN
1781: *
1782: * MNTHR2 <= N < MNTHR1
1783: *
1784: * Path 5t (N >> M, but not as much as MNTHR1)
1785: * Reduce to bidiagonal form without QR decomposition, use
1786: * ZUNGBR and matrix multiplication to compute singular vectors
1787: *
1788: IE = 1
1789: NRWORK = IE + M
1790: ITAUQ = 1
1791: ITAUP = ITAUQ + M
1792: NWORK = ITAUP + M
1793: *
1794: * Bidiagonalize A
1795: * CWorkspace: need 2*M [tauq, taup] + N [work]
1796: * CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
1797: * RWorkspace: need M [e]
1798: *
1799: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1800: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1801: $ IERR )
1802: *
1803: IF( WNTQN ) THEN
1804: *
1805: * Path 5tn (N >> M, JOBZ='N')
1806: * Compute singular values only
1807: * CWorkspace: need 0
1808: * RWorkspace: need M [e] + BDSPAC
1809: *
1810: CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
1811: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1812: ELSE IF( WNTQO ) THEN
1813: IRVT = NRWORK
1814: IRU = IRVT + M*M
1815: NRWORK = IRU + M*M
1816: IVT = NWORK
1817: *
1818: * Path 5to (N >> M, JOBZ='O')
1819: * Copy A to U, generate Q
1820: * CWorkspace: need 2*M [tauq, taup] + M [work]
1821: * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1822: * RWorkspace: need 0
1823: *
1824: CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
1825: CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1826: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1827: *
1828: * Generate P**H in A
1829: * CWorkspace: need 2*M [tauq, taup] + M [work]
1830: * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1831: * RWorkspace: need 0
1832: *
1833: CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1834: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1835: *
1836: LDWKVT = M
1837: IF( LWORK .GE. M*N + 3*M ) THEN
1838: *
1839: * WORK( IVT ) is M by N
1840: *
1841: NWORK = IVT + LDWKVT*N
1842: CHUNK = N
1843: ELSE
1844: *
1845: * WORK( IVT ) is M by CHUNK
1846: *
1847: CHUNK = ( LWORK - 3*M ) / M
1848: NWORK = IVT + LDWKVT*CHUNK
1849: END IF
1850: *
1851: * Perform bidiagonal SVD, computing left singular vectors
1852: * of bidiagonal matrix in RWORK(IRU) and computing right
1853: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1854: * CWorkspace: need 0
1855: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1856: *
1857: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1858: $ M, RWORK( IRVT ), M, DUM, IDUM,
1859: $ RWORK( NRWORK ), IWORK, INFO )
1860: *
1861: * Multiply Q in U by real matrix RWORK(IRVT)
1862: * storing the result in WORK(IVT), copying to U
1863: * CWorkspace: need 2*M [tauq, taup] + M*M [VT]
1864: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1865: *
1866: CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ),
1867: $ LDWKVT, RWORK( NRWORK ) )
1868: CALL ZLACPY( 'F', M, M, WORK( IVT ), LDWKVT, U, LDU )
1869: *
1870: * Multiply RWORK(IRVT) by P**H in A, storing the
1871: * result in WORK(IVT), copying to A
1872: * CWorkspace: need 2*M [tauq, taup] + M*M [VT]
1873: * CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
1874: * RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork]
1875: * RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1876: *
1877: NRWORK = IRU
1878: DO 50 I = 1, N, CHUNK
1879: BLK = MIN( N-I+1, CHUNK )
1880: CALL ZLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), LDA,
1881: $ WORK( IVT ), LDWKVT, RWORK( NRWORK ) )
1882: CALL ZLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
1883: $ A( 1, I ), LDA )
1884: 50 CONTINUE
1885: ELSE IF( WNTQS ) THEN
1886: *
1887: * Path 5ts (N >> M, JOBZ='S')
1888: * Copy A to U, generate Q
1889: * CWorkspace: need 2*M [tauq, taup] + M [work]
1890: * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1891: * RWorkspace: need 0
1892: *
1893: CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
1894: CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1895: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1896: *
1897: * Copy A to VT, generate P**H
1898: * CWorkspace: need 2*M [tauq, taup] + M [work]
1899: * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1900: * RWorkspace: need 0
1901: *
1902: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
1903: CALL ZUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ),
1904: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1905: *
1906: * Perform bidiagonal SVD, computing left singular vectors
1907: * of bidiagonal matrix in RWORK(IRU) and computing right
1908: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1909: * CWorkspace: need 0
1910: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1911: *
1912: IRVT = NRWORK
1913: IRU = IRVT + M*M
1914: NRWORK = IRU + M*M
1915: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1916: $ M, RWORK( IRVT ), M, DUM, IDUM,
1917: $ RWORK( NRWORK ), IWORK, INFO )
1918: *
1919: * Multiply Q in U by real matrix RWORK(IRU), storing the
1920: * result in A, copying to U
1921: * CWorkspace: need 0
1922: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1923: *
1924: CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
1925: $ RWORK( NRWORK ) )
1926: CALL ZLACPY( 'F', M, M, A, LDA, U, LDU )
1927: *
1928: * Multiply real matrix RWORK(IRVT) by P**H in VT,
1929: * storing the result in A, copying to VT
1930: * CWorkspace: need 0
1931: * RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1932: *
1933: NRWORK = IRU
1934: CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
1935: $ RWORK( NRWORK ) )
1936: CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
1937: ELSE
1938: *
1939: * Path 5ta (N >> M, JOBZ='A')
1940: * Copy A to U, generate Q
1941: * CWorkspace: need 2*M [tauq, taup] + M [work]
1942: * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1943: * RWorkspace: need 0
1944: *
1945: CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
1946: CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1947: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1948: *
1949: * Copy A to VT, generate P**H
1950: * CWorkspace: need 2*M [tauq, taup] + N [work]
1951: * CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
1952: * RWorkspace: need 0
1953: *
1954: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
1955: CALL ZUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ),
1956: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1957: *
1958: * Perform bidiagonal SVD, computing left singular vectors
1959: * of bidiagonal matrix in RWORK(IRU) and computing right
1960: * singular vectors of bidiagonal matrix in RWORK(IRVT)
1961: * CWorkspace: need 0
1962: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1963: *
1964: IRVT = NRWORK
1965: IRU = IRVT + M*M
1966: NRWORK = IRU + M*M
1967: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1968: $ M, RWORK( IRVT ), M, DUM, IDUM,
1969: $ RWORK( NRWORK ), IWORK, INFO )
1970: *
1971: * Multiply Q in U by real matrix RWORK(IRU), storing the
1972: * result in A, copying to U
1973: * CWorkspace: need 0
1974: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1975: *
1976: CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
1977: $ RWORK( NRWORK ) )
1978: CALL ZLACPY( 'F', M, M, A, LDA, U, LDU )
1979: *
1980: * Multiply real matrix RWORK(IRVT) by P**H in VT,
1981: * storing the result in A, copying to VT
1982: * CWorkspace: need 0
1983: * RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1984: *
1985: NRWORK = IRU
1986: CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
1987: $ RWORK( NRWORK ) )
1988: CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
1989: END IF
1990: *
1991: ELSE
1992: *
1993: * N .LT. MNTHR2
1994: *
1995: * Path 6t (N > M, but not much larger)
1996: * Reduce to bidiagonal form without LQ decomposition
1997: * Use ZUNMBR to compute singular vectors
1998: *
1999: IE = 1
2000: NRWORK = IE + M
2001: ITAUQ = 1
2002: ITAUP = ITAUQ + M
2003: NWORK = ITAUP + M
2004: *
2005: * Bidiagonalize A
2006: * CWorkspace: need 2*M [tauq, taup] + N [work]
2007: * CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
2008: * RWorkspace: need M [e]
2009: *
2010: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
2011: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
2012: $ IERR )
2013: IF( WNTQN ) THEN
2014: *
2015: * Path 6tn (N > M, JOBZ='N')
2016: * Compute singular values only
2017: * CWorkspace: need 0
2018: * RWorkspace: need M [e] + BDSPAC
2019: *
2020: CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM,1,DUM,1,
2021: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
2022: ELSE IF( WNTQO ) THEN
2023: * Path 6to (N > M, JOBZ='O')
2024: LDWKVT = M
2025: IVT = NWORK
2026: IF( LWORK .GE. M*N + 3*M ) THEN
2027: *
2028: * WORK( IVT ) is M by N
2029: *
2030: CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IVT ),
2031: $ LDWKVT )
2032: NWORK = IVT + LDWKVT*N
2033: ELSE
2034: *
2035: * WORK( IVT ) is M by CHUNK
2036: *
2037: CHUNK = ( LWORK - 3*M ) / M
2038: NWORK = IVT + LDWKVT*CHUNK
2039: END IF
2040: *
2041: * Perform bidiagonal SVD, computing left singular vectors
2042: * of bidiagonal matrix in RWORK(IRU) and computing right
2043: * singular vectors of bidiagonal matrix in RWORK(IRVT)
2044: * CWorkspace: need 0
2045: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2046: *
2047: IRVT = NRWORK
2048: IRU = IRVT + M*M
2049: NRWORK = IRU + M*M
2050: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
2051: $ M, RWORK( IRVT ), M, DUM, IDUM,
2052: $ RWORK( NRWORK ), IWORK, INFO )
2053: *
2054: * Copy real matrix RWORK(IRU) to complex matrix U
2055: * Overwrite U by left singular vectors of A
2056: * CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work]
2057: * CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
2058: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2059: *
2060: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
2061: CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
2062: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
2063: $ LWORK-NWORK+1, IERR )
2064: *
2065: IF( LWORK .GE. M*N + 3*M ) THEN
2066: *
2067: * Path 6to-fast
2068: * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
2069: * Overwrite WORK(IVT) by right singular vectors of A,
2070: * copying to A
2071: * CWorkspace: need 2*M [tauq, taup] + M*N [VT] + M [work]
2072: * CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] + M*NB [work]
2073: * RWorkspace: need M [e] + M*M [RVT]
2074: *
2075: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
2076: $ LDWKVT )
2077: CALL ZUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
2078: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
2079: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
2080: CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
2081: ELSE
2082: *
2083: * Path 6to-slow
2084: * Generate P**H in A
2085: * CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work]
2086: * CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
2087: * RWorkspace: need 0
2088: *
2089: CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
2090: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
2091: *
2092: * Multiply Q in A by real matrix RWORK(IRU), storing the
2093: * result in WORK(IU), copying to A
2094: * CWorkspace: need 2*M [tauq, taup] + M*M [VT]
2095: * CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
2096: * RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork]
2097: * RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
2098: *
2099: NRWORK = IRU
2100: DO 60 I = 1, N, CHUNK
2101: BLK = MIN( N-I+1, CHUNK )
2102: CALL ZLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ),
2103: $ LDA, WORK( IVT ), LDWKVT,
2104: $ RWORK( NRWORK ) )
2105: CALL ZLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
2106: $ A( 1, I ), LDA )
2107: 60 CONTINUE
2108: END IF
2109: ELSE IF( WNTQS ) THEN
2110: *
2111: * Path 6ts (N > M, JOBZ='S')
2112: * Perform bidiagonal SVD, computing left singular vectors
2113: * of bidiagonal matrix in RWORK(IRU) and computing right
2114: * singular vectors of bidiagonal matrix in RWORK(IRVT)
2115: * CWorkspace: need 0
2116: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2117: *
2118: IRVT = NRWORK
2119: IRU = IRVT + M*M
2120: NRWORK = IRU + M*M
2121: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
2122: $ M, RWORK( IRVT ), M, DUM, IDUM,
2123: $ RWORK( NRWORK ), IWORK, INFO )
2124: *
2125: * Copy real matrix RWORK(IRU) to complex matrix U
2126: * Overwrite U by left singular vectors of A
2127: * CWorkspace: need 2*M [tauq, taup] + M [work]
2128: * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2129: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2130: *
2131: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
2132: CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
2133: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
2134: $ LWORK-NWORK+1, IERR )
2135: *
2136: * Copy real matrix RWORK(IRVT) to complex matrix VT
2137: * Overwrite VT by right singular vectors of A
2138: * CWorkspace: need 2*M [tauq, taup] + M [work]
2139: * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2140: * RWorkspace: need M [e] + M*M [RVT]
2141: *
2142: CALL ZLASET( 'F', M, N, CZERO, CZERO, VT, LDVT )
2143: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
2144: CALL ZUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
2145: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
2146: $ LWORK-NWORK+1, IERR )
2147: ELSE
2148: *
2149: * Path 6ta (N > M, JOBZ='A')
2150: * Perform bidiagonal SVD, computing left singular vectors
2151: * of bidiagonal matrix in RWORK(IRU) and computing right
2152: * singular vectors of bidiagonal matrix in RWORK(IRVT)
2153: * CWorkspace: need 0
2154: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2155: *
2156: IRVT = NRWORK
2157: IRU = IRVT + M*M
2158: NRWORK = IRU + M*M
2159: *
2160: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
2161: $ M, RWORK( IRVT ), M, DUM, IDUM,
2162: $ RWORK( NRWORK ), IWORK, INFO )
2163: *
2164: * Copy real matrix RWORK(IRU) to complex matrix U
2165: * Overwrite U by left singular vectors of A
2166: * CWorkspace: need 2*M [tauq, taup] + M [work]
2167: * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2168: * RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2169: *
2170: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
2171: CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
2172: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
2173: $ LWORK-NWORK+1, IERR )
2174: *
2175: * Set all of VT to identity matrix
2176: *
2177: CALL ZLASET( 'F', N, N, CZERO, CONE, VT, LDVT )
2178: *
2179: * Copy real matrix RWORK(IRVT) to complex matrix VT
2180: * Overwrite VT by right singular vectors of A
2181: * CWorkspace: need 2*M [tauq, taup] + N [work]
2182: * CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
2183: * RWorkspace: need M [e] + M*M [RVT]
2184: *
2185: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
2186: CALL ZUNMBR( 'P', 'R', 'C', N, N, M, A, LDA,
2187: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
2188: $ LWORK-NWORK+1, IERR )
2189: END IF
2190: *
2191: END IF
2192: *
2193: END IF
2194: *
2195: * Undo scaling if necessary
2196: *
2197: IF( ISCL.EQ.1 ) THEN
2198: IF( ANRM.GT.BIGNUM )
2199: $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
2200: $ IERR )
2201: IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
2202: $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1,
2203: $ RWORK( IE ), MINMN, IERR )
2204: IF( ANRM.LT.SMLNUM )
2205: $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
2206: $ IERR )
2207: IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
2208: $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1,
2209: $ RWORK( IE ), MINMN, IERR )
2210: END IF
2211: *
2212: * Return optimal workspace in WORK(1)
2213: *
2214: WORK( 1 ) = MAXWRK
2215: *
2216: RETURN
2217: *
2218: * End of ZGESDD
2219: *
2220: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgesdd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack