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