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