1: *> \brief \b DGESDD
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DGESDD + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesdd.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesdd.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesdd.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
22: * LWORK, 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 A( LDA, * ), S( * ), U( LDU, * ),
31: * $ VT( LDVT, * ), WORK( * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> DGESDD computes the singular value decomposition (SVD) of a real
41: *> M-by-N matrix A, optionally computing the left and right singular
42: *> vectors. If singular vectors are desired, it uses a
43: *> divide-and-conquer algorithm.
44: *>
45: *> The SVD is written
46: *>
47: *> A = U * SIGMA * transpose(V)
48: *>
49: *> where SIGMA is an M-by-N matrix which is zero except for its
50: *> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
51: *> V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
52: *> are the singular values of A; they are real and non-negative, and
53: *> are returned in descending order. The first min(m,n) columns of
54: *> U and V are the left and right singular vectors of A.
55: *>
56: *> Note that the routine returns VT = V**T, not V.
57: *>
58: *> The divide and conquer algorithm makes very mild assumptions about
59: *> floating point arithmetic. It will work on machines with a guard
60: *> digit in add/subtract, or on those binary machines without guard
61: *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
62: *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
63: *> without guard digits, but we know of none.
64: *> \endverbatim
65: *
66: * Arguments:
67: * ==========
68: *
69: *> \param[in] JOBZ
70: *> \verbatim
71: *> JOBZ is CHARACTER*1
72: *> Specifies options for computing all or part of the matrix U:
73: *> = 'A': all M columns of U and all N rows of V**T are
74: *> returned in the arrays U and VT;
75: *> = 'S': the first min(M,N) columns of U and the first
76: *> min(M,N) rows of V**T are returned in the arrays U
77: *> and VT;
78: *> = 'O': If M >= N, the first N columns of U are overwritten
79: *> on the array A and all rows of V**T are returned in
80: *> the array VT;
81: *> otherwise, all columns of U are returned in the
82: *> array U and the first M rows of V**T are overwritten
83: *> in the array A;
84: *> = 'N': no columns of U or rows of V**T are computed.
85: *> \endverbatim
86: *>
87: *> \param[in] M
88: *> \verbatim
89: *> M is INTEGER
90: *> The number of rows of the input matrix A. M >= 0.
91: *> \endverbatim
92: *>
93: *> \param[in] N
94: *> \verbatim
95: *> N is INTEGER
96: *> The number of columns of the input matrix A. N >= 0.
97: *> \endverbatim
98: *>
99: *> \param[in,out] A
100: *> \verbatim
101: *> A is DOUBLE PRECISION array, dimension (LDA,N)
102: *> On entry, the M-by-N matrix A.
103: *> On exit,
104: *> if JOBZ = 'O', A is overwritten with the first N columns
105: *> of U (the left singular vectors, stored
106: *> columnwise) if M >= N;
107: *> A is overwritten with the first M rows
108: *> of V**T (the right singular vectors, stored
109: *> rowwise) otherwise.
110: *> if JOBZ .ne. 'O', the contents of A are destroyed.
111: *> \endverbatim
112: *>
113: *> \param[in] LDA
114: *> \verbatim
115: *> LDA is INTEGER
116: *> The leading dimension of the array A. LDA >= max(1,M).
117: *> \endverbatim
118: *>
119: *> \param[out] S
120: *> \verbatim
121: *> S is DOUBLE PRECISION array, dimension (min(M,N))
122: *> The singular values of A, sorted so that S(i) >= S(i+1).
123: *> \endverbatim
124: *>
125: *> \param[out] U
126: *> \verbatim
127: *> U is DOUBLE PRECISION array, dimension (LDU,UCOL)
128: *> UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
129: *> UCOL = min(M,N) if JOBZ = 'S'.
130: *> If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
131: *> orthogonal matrix U;
132: *> if JOBZ = 'S', U contains the first min(M,N) columns of U
133: *> (the left singular vectors, stored columnwise);
134: *> if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
135: *> \endverbatim
136: *>
137: *> \param[in] LDU
138: *> \verbatim
139: *> LDU is INTEGER
140: *> The leading dimension of the array U. LDU >= 1; if
141: *> JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
142: *> \endverbatim
143: *>
144: *> \param[out] VT
145: *> \verbatim
146: *> VT is DOUBLE PRECISION array, dimension (LDVT,N)
147: *> If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
148: *> N-by-N orthogonal matrix V**T;
149: *> if JOBZ = 'S', VT contains the first min(M,N) rows of
150: *> V**T (the right singular vectors, stored rowwise);
151: *> if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
152: *> \endverbatim
153: *>
154: *> \param[in] LDVT
155: *> \verbatim
156: *> LDVT is INTEGER
157: *> The leading dimension of the array VT. LDVT >= 1; if
158: *> JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
159: *> if JOBZ = 'S', LDVT >= min(M,N).
160: *> \endverbatim
161: *>
162: *> \param[out] WORK
163: *> \verbatim
164: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
165: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
166: *> \endverbatim
167: *>
168: *> \param[in] LWORK
169: *> \verbatim
170: *> LWORK is INTEGER
171: *> The dimension of the array WORK. LWORK >= 1.
172: *> If JOBZ = 'N',
173: *> LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)).
174: *> If JOBZ = 'O',
175: *> LWORK >= 3*min(M,N) +
176: *> max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
177: *> If JOBZ = 'S' or 'A'
178: *> LWORK >= 3*min(M,N) +
179: *> max(max(M,N),4*min(M,N)*min(M,N)+3*min(M,N)+max(M,N)).
180: *> For good performance, LWORK should generally be larger.
181: *> If LWORK = -1 but other input arguments are legal, WORK(1)
182: *> returns the optimal LWORK.
183: *> \endverbatim
184: *>
185: *> \param[out] IWORK
186: *> \verbatim
187: *> IWORK is INTEGER array, dimension (8*min(M,N))
188: *> \endverbatim
189: *>
190: *> \param[out] INFO
191: *> \verbatim
192: *> INFO is INTEGER
193: *> = 0: successful exit.
194: *> < 0: if INFO = -i, the i-th argument had an illegal value.
195: *> > 0: DBDSDC did not converge, updating process failed.
196: *> \endverbatim
197: *
198: * Authors:
199: * ========
200: *
201: *> \author Univ. of Tennessee
202: *> \author Univ. of California Berkeley
203: *> \author Univ. of Colorado Denver
204: *> \author NAG Ltd.
205: *
206: *> \date September 2012
207: *
208: *> \ingroup doubleGEsing
209: *
210: *> \par Contributors:
211: * ==================
212: *>
213: *> Ming Gu and Huan Ren, Computer Science Division, University of
214: *> California at Berkeley, USA
215: *>
216: * =====================================================================
217: SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
218: $ LWORK, IWORK, INFO )
219: *
220: * -- LAPACK driver routine (version 3.4.2) --
221: * -- LAPACK is a software package provided by Univ. of Tennessee, --
222: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
223: * September 2012
224: *
225: * .. Scalar Arguments ..
226: CHARACTER JOBZ
227: INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
228: * ..
229: * .. Array Arguments ..
230: INTEGER IWORK( * )
231: DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
232: $ VT( LDVT, * ), WORK( * )
233: * ..
234: *
235: * =====================================================================
236: *
237: * .. Parameters ..
238: DOUBLE PRECISION ZERO, ONE
239: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
240: * ..
241: * .. Local Scalars ..
242: LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
243: INTEGER BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
244: $ IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
245: $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
246: $ MNTHR, NWORK, WRKBL
247: DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
248: * ..
249: * .. Local Arrays ..
250: INTEGER IDUM( 1 )
251: DOUBLE PRECISION DUM( 1 )
252: * ..
253: * .. External Subroutines ..
254: EXTERNAL DBDSDC, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY,
255: $ DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR,
256: $ XERBLA
257: * ..
258: * .. External Functions ..
259: LOGICAL LSAME
260: INTEGER ILAENV
261: DOUBLE PRECISION DLAMCH, DLANGE
262: EXTERNAL DLAMCH, DLANGE, ILAENV, LSAME
263: * ..
264: * .. Intrinsic Functions ..
265: INTRINSIC INT, MAX, MIN, SQRT
266: * ..
267: * .. Executable Statements ..
268: *
269: * Test the input arguments
270: *
271: INFO = 0
272: MINMN = MIN( M, N )
273: WNTQA = LSAME( JOBZ, 'A' )
274: WNTQS = LSAME( JOBZ, 'S' )
275: WNTQAS = WNTQA .OR. WNTQS
276: WNTQO = LSAME( JOBZ, 'O' )
277: WNTQN = LSAME( JOBZ, 'N' )
278: LQUERY = ( LWORK.EQ.-1 )
279: *
280: IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
281: INFO = -1
282: ELSE IF( M.LT.0 ) THEN
283: INFO = -2
284: ELSE IF( N.LT.0 ) THEN
285: INFO = -3
286: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
287: INFO = -5
288: ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
289: $ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
290: INFO = -8
291: ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
292: $ ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
293: $ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
294: INFO = -10
295: END IF
296: *
297: * Compute workspace
298: * (Note: Comments in the code beginning "Workspace:" describe the
299: * minimal amount of workspace needed at that point in the code,
300: * as well as the preferred amount for good performance.
301: * NB refers to the optimal block size for the immediately
302: * following subroutine, as returned by ILAENV.)
303: *
304: IF( INFO.EQ.0 ) THEN
305: MINWRK = 1
306: MAXWRK = 1
307: IF( M.GE.N .AND. MINMN.GT.0 ) THEN
308: *
309: * Compute space needed for DBDSDC
310: *
311: MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
312: IF( WNTQN ) THEN
313: BDSPAC = 7*N
314: ELSE
315: BDSPAC = 3*N*N + 4*N
316: END IF
317: IF( M.GE.MNTHR ) THEN
318: IF( WNTQN ) THEN
319: *
320: * Path 1 (M much larger than N, JOBZ='N')
321: *
322: WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1,
323: $ -1 )
324: WRKBL = MAX( WRKBL, 3*N+2*N*
325: $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
326: MAXWRK = MAX( WRKBL, BDSPAC+N )
327: MINWRK = BDSPAC + N
328: ELSE IF( WNTQO ) THEN
329: *
330: * Path 2 (M much larger than N, JOBZ='O')
331: *
332: WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
333: WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
334: $ N, N, -1 ) )
335: WRKBL = MAX( WRKBL, 3*N+2*N*
336: $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
337: WRKBL = MAX( WRKBL, 3*N+N*
338: $ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
339: WRKBL = MAX( WRKBL, 3*N+N*
340: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
341: WRKBL = MAX( WRKBL, BDSPAC+3*N )
342: MAXWRK = WRKBL + 2*N*N
343: MINWRK = BDSPAC + 2*N*N + 3*N
344: ELSE IF( WNTQS ) THEN
345: *
346: * Path 3 (M much larger than N, JOBZ='S')
347: *
348: WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
349: WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
350: $ N, N, -1 ) )
351: WRKBL = MAX( WRKBL, 3*N+2*N*
352: $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
353: WRKBL = MAX( WRKBL, 3*N+N*
354: $ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
355: WRKBL = MAX( WRKBL, 3*N+N*
356: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
357: WRKBL = MAX( WRKBL, BDSPAC+3*N )
358: MAXWRK = WRKBL + N*N
359: MINWRK = BDSPAC + N*N + 3*N
360: ELSE IF( WNTQA ) THEN
361: *
362: * Path 4 (M much larger than N, JOBZ='A')
363: *
364: WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
365: WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,
366: $ M, N, -1 ) )
367: WRKBL = MAX( WRKBL, 3*N+2*N*
368: $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
369: WRKBL = MAX( WRKBL, 3*N+N*
370: $ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
371: WRKBL = MAX( WRKBL, 3*N+N*
372: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
373: WRKBL = MAX( WRKBL, BDSPAC+3*N )
374: MAXWRK = WRKBL + N*N
375: MINWRK = BDSPAC + N*N + 2*N + M
376: END IF
377: ELSE
378: *
379: * Path 5 (M at least N, but not much larger)
380: *
381: WRKBL = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
382: $ -1 )
383: IF( WNTQN ) THEN
384: MAXWRK = MAX( WRKBL, BDSPAC+3*N )
385: MINWRK = 3*N + MAX( M, BDSPAC )
386: ELSE IF( WNTQO ) THEN
387: WRKBL = MAX( WRKBL, 3*N+N*
388: $ ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
389: WRKBL = MAX( WRKBL, 3*N+N*
390: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
391: WRKBL = MAX( WRKBL, BDSPAC+3*N )
392: MAXWRK = WRKBL + M*N
393: MINWRK = 3*N + MAX( M, N*N+BDSPAC )
394: ELSE IF( WNTQS ) THEN
395: WRKBL = MAX( WRKBL, 3*N+N*
396: $ ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
397: WRKBL = MAX( WRKBL, 3*N+N*
398: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
399: MAXWRK = MAX( WRKBL, BDSPAC+3*N )
400: MINWRK = 3*N + MAX( M, BDSPAC )
401: ELSE IF( WNTQA ) THEN
402: WRKBL = MAX( WRKBL, 3*N+M*
403: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
404: WRKBL = MAX( WRKBL, 3*N+N*
405: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
406: MAXWRK = MAX( MAXWRK, BDSPAC+3*N )
407: MINWRK = 3*N + MAX( M, BDSPAC )
408: END IF
409: END IF
410: ELSE IF( MINMN.GT.0 ) THEN
411: *
412: * Compute space needed for DBDSDC
413: *
414: MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
415: IF( WNTQN ) THEN
416: BDSPAC = 7*M
417: ELSE
418: BDSPAC = 3*M*M + 4*M
419: END IF
420: IF( N.GE.MNTHR ) THEN
421: IF( WNTQN ) THEN
422: *
423: * Path 1t (N much larger than M, JOBZ='N')
424: *
425: WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1,
426: $ -1 )
427: WRKBL = MAX( WRKBL, 3*M+2*M*
428: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
429: MAXWRK = MAX( WRKBL, BDSPAC+M )
430: MINWRK = BDSPAC + M
431: ELSE IF( WNTQO ) THEN
432: *
433: * Path 2t (N much larger than M, JOBZ='O')
434: *
435: WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
436: WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
437: $ N, M, -1 ) )
438: WRKBL = MAX( WRKBL, 3*M+2*M*
439: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
440: WRKBL = MAX( WRKBL, 3*M+M*
441: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
442: WRKBL = MAX( WRKBL, 3*M+M*
443: $ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
444: WRKBL = MAX( WRKBL, BDSPAC+3*M )
445: MAXWRK = WRKBL + 2*M*M
446: MINWRK = BDSPAC + 2*M*M + 3*M
447: ELSE IF( WNTQS ) THEN
448: *
449: * Path 3t (N much larger than M, JOBZ='S')
450: *
451: WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
452: WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
453: $ N, M, -1 ) )
454: WRKBL = MAX( WRKBL, 3*M+2*M*
455: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
456: WRKBL = MAX( WRKBL, 3*M+M*
457: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
458: WRKBL = MAX( WRKBL, 3*M+M*
459: $ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
460: WRKBL = MAX( WRKBL, BDSPAC+3*M )
461: MAXWRK = WRKBL + M*M
462: MINWRK = BDSPAC + M*M + 3*M
463: ELSE IF( WNTQA ) THEN
464: *
465: * Path 4t (N much larger than M, JOBZ='A')
466: *
467: WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
468: WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,
469: $ N, M, -1 ) )
470: WRKBL = MAX( WRKBL, 3*M+2*M*
471: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
472: WRKBL = MAX( WRKBL, 3*M+M*
473: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
474: WRKBL = MAX( WRKBL, 3*M+M*
475: $ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
476: WRKBL = MAX( WRKBL, BDSPAC+3*M )
477: MAXWRK = WRKBL + M*M
478: MINWRK = BDSPAC + M*M + 3*M
479: END IF
480: ELSE
481: *
482: * Path 5t (N greater than M, but not much larger)
483: *
484: WRKBL = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
485: $ -1 )
486: IF( WNTQN ) THEN
487: MAXWRK = MAX( WRKBL, BDSPAC+3*M )
488: MINWRK = 3*M + MAX( N, BDSPAC )
489: ELSE IF( WNTQO ) THEN
490: WRKBL = MAX( WRKBL, 3*M+M*
491: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
492: WRKBL = MAX( WRKBL, 3*M+M*
493: $ ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
494: WRKBL = MAX( WRKBL, BDSPAC+3*M )
495: MAXWRK = WRKBL + M*N
496: MINWRK = 3*M + MAX( N, M*M+BDSPAC )
497: ELSE IF( WNTQS ) THEN
498: WRKBL = MAX( WRKBL, 3*M+M*
499: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
500: WRKBL = MAX( WRKBL, 3*M+M*
501: $ ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
502: MAXWRK = MAX( WRKBL, BDSPAC+3*M )
503: MINWRK = 3*M + MAX( N, BDSPAC )
504: ELSE IF( WNTQA ) THEN
505: WRKBL = MAX( WRKBL, 3*M+M*
506: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
507: WRKBL = MAX( WRKBL, 3*M+M*
508: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, M, -1 ) )
509: MAXWRK = MAX( WRKBL, BDSPAC+3*M )
510: MINWRK = 3*M + MAX( N, BDSPAC )
511: END IF
512: END IF
513: END IF
514: MAXWRK = MAX( MAXWRK, MINWRK )
515: WORK( 1 ) = MAXWRK
516: *
517: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
518: INFO = -12
519: END IF
520: END IF
521: *
522: IF( INFO.NE.0 ) THEN
523: CALL XERBLA( 'DGESDD', -INFO )
524: RETURN
525: ELSE IF( LQUERY ) THEN
526: RETURN
527: END IF
528: *
529: * Quick return if possible
530: *
531: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
532: RETURN
533: END IF
534: *
535: * Get machine constants
536: *
537: EPS = DLAMCH( 'P' )
538: SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
539: BIGNUM = ONE / SMLNUM
540: *
541: * Scale A if max element outside range [SMLNUM,BIGNUM]
542: *
543: ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
544: ISCL = 0
545: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
546: ISCL = 1
547: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
548: ELSE IF( ANRM.GT.BIGNUM ) THEN
549: ISCL = 1
550: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
551: END IF
552: *
553: IF( M.GE.N ) THEN
554: *
555: * A has at least as many rows as columns. If A has sufficiently
556: * more rows than columns, first reduce using the QR
557: * decomposition (if sufficient workspace available)
558: *
559: IF( M.GE.MNTHR ) THEN
560: *
561: IF( WNTQN ) THEN
562: *
563: * Path 1 (M much larger than N, JOBZ='N')
564: * No singular vectors to be computed
565: *
566: ITAU = 1
567: NWORK = ITAU + N
568: *
569: * Compute A=Q*R
570: * (Workspace: need 2*N, prefer N+N*NB)
571: *
572: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
573: $ LWORK-NWORK+1, IERR )
574: *
575: * Zero out below R
576: *
577: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
578: IE = 1
579: ITAUQ = IE + N
580: ITAUP = ITAUQ + N
581: NWORK = ITAUP + N
582: *
583: * Bidiagonalize R in A
584: * (Workspace: need 4*N, prefer 3*N+2*N*NB)
585: *
586: CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
587: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
588: $ IERR )
589: NWORK = IE + N
590: *
591: * Perform bidiagonal SVD, computing singular values only
592: * (Workspace: need N+BDSPAC)
593: *
594: CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
595: $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
596: *
597: ELSE IF( WNTQO ) THEN
598: *
599: * Path 2 (M much larger than N, JOBZ = 'O')
600: * N left singular vectors to be overwritten on A and
601: * N right singular vectors to be computed in VT
602: *
603: IR = 1
604: *
605: * WORK(IR) is LDWRKR by N
606: *
607: IF( LWORK.GE.LDA*N+N*N+3*N+BDSPAC ) THEN
608: LDWRKR = LDA
609: ELSE
610: LDWRKR = ( LWORK-N*N-3*N-BDSPAC ) / N
611: END IF
612: ITAU = IR + LDWRKR*N
613: NWORK = ITAU + N
614: *
615: * Compute A=Q*R
616: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
617: *
618: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
619: $ LWORK-NWORK+1, IERR )
620: *
621: * Copy R to WORK(IR), zeroing out below it
622: *
623: CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
624: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
625: $ LDWRKR )
626: *
627: * Generate Q in A
628: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
629: *
630: CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
631: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
632: IE = ITAU
633: ITAUQ = IE + N
634: ITAUP = ITAUQ + N
635: NWORK = ITAUP + N
636: *
637: * Bidiagonalize R in VT, copying result to WORK(IR)
638: * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
639: *
640: CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
641: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
642: $ LWORK-NWORK+1, IERR )
643: *
644: * WORK(IU) is N by N
645: *
646: IU = NWORK
647: NWORK = IU + N*N
648: *
649: * Perform bidiagonal SVD, computing left singular vectors
650: * of bidiagonal matrix in WORK(IU) and computing right
651: * singular vectors of bidiagonal matrix in VT
652: * (Workspace: need N+N*N+BDSPAC)
653: *
654: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
655: $ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
656: $ INFO )
657: *
658: * Overwrite WORK(IU) by left singular vectors of R
659: * and VT by right singular vectors of R
660: * (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
661: *
662: CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
663: $ WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
664: $ LWORK-NWORK+1, IERR )
665: CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
666: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
667: $ LWORK-NWORK+1, IERR )
668: *
669: * Multiply Q in A by left singular vectors of R in
670: * WORK(IU), storing result in WORK(IR) and copying to A
671: * (Workspace: need 2*N*N, prefer N*N+M*N)
672: *
673: DO 10 I = 1, M, LDWRKR
674: CHUNK = MIN( M-I+1, LDWRKR )
675: CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
676: $ LDA, WORK( IU ), N, ZERO, WORK( IR ),
677: $ LDWRKR )
678: CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
679: $ A( I, 1 ), LDA )
680: 10 CONTINUE
681: *
682: ELSE IF( WNTQS ) THEN
683: *
684: * Path 3 (M much larger than N, JOBZ='S')
685: * N left singular vectors to be computed in U and
686: * N right singular vectors to be computed in VT
687: *
688: IR = 1
689: *
690: * WORK(IR) is N by N
691: *
692: LDWRKR = N
693: ITAU = IR + LDWRKR*N
694: NWORK = ITAU + N
695: *
696: * Compute A=Q*R
697: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
698: *
699: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
700: $ LWORK-NWORK+1, IERR )
701: *
702: * Copy R to WORK(IR), zeroing out below it
703: *
704: CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
705: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
706: $ LDWRKR )
707: *
708: * Generate Q in A
709: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
710: *
711: CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
712: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
713: IE = ITAU
714: ITAUQ = IE + N
715: ITAUP = ITAUQ + N
716: NWORK = ITAUP + N
717: *
718: * Bidiagonalize R in WORK(IR)
719: * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
720: *
721: CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
722: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
723: $ LWORK-NWORK+1, IERR )
724: *
725: * Perform bidiagonal SVD, computing left singular vectors
726: * of bidiagoal matrix in U and computing right singular
727: * vectors of bidiagonal matrix in VT
728: * (Workspace: need N+BDSPAC)
729: *
730: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
731: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
732: $ INFO )
733: *
734: * Overwrite U by left singular vectors of R and VT
735: * by right singular vectors of R
736: * (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
737: *
738: CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
739: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
740: $ LWORK-NWORK+1, IERR )
741: *
742: CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
743: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
744: $ LWORK-NWORK+1, IERR )
745: *
746: * Multiply Q in A by left singular vectors of R in
747: * WORK(IR), storing result in U
748: * (Workspace: need N*N)
749: *
750: CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
751: CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),
752: $ LDWRKR, ZERO, U, LDU )
753: *
754: ELSE IF( WNTQA ) THEN
755: *
756: * Path 4 (M much larger than N, JOBZ='A')
757: * M left singular vectors to be computed in U and
758: * N right singular vectors to be computed in VT
759: *
760: IU = 1
761: *
762: * WORK(IU) is N by N
763: *
764: LDWRKU = N
765: ITAU = IU + LDWRKU*N
766: NWORK = ITAU + N
767: *
768: * Compute A=Q*R, copying result to U
769: * (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
770: *
771: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
772: $ LWORK-NWORK+1, IERR )
773: CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
774: *
775: * Generate Q in U
776: * (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
777: CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
778: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
779: *
780: * Produce R in A, zeroing out other entries
781: *
782: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
783: IE = ITAU
784: ITAUQ = IE + N
785: ITAUP = ITAUQ + N
786: NWORK = ITAUP + N
787: *
788: * Bidiagonalize R in A
789: * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
790: *
791: CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
792: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
793: $ IERR )
794: *
795: * Perform bidiagonal SVD, computing left singular vectors
796: * of bidiagonal matrix in WORK(IU) and computing right
797: * singular vectors of bidiagonal matrix in VT
798: * (Workspace: need N+N*N+BDSPAC)
799: *
800: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
801: $ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
802: $ INFO )
803: *
804: * Overwrite WORK(IU) by left singular vectors of R and VT
805: * by right singular vectors of R
806: * (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
807: *
808: CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
809: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
810: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
811: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
812: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
813: $ LWORK-NWORK+1, IERR )
814: *
815: * Multiply Q in U by left singular vectors of R in
816: * WORK(IU), storing result in A
817: * (Workspace: need N*N)
818: *
819: CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),
820: $ LDWRKU, ZERO, A, LDA )
821: *
822: * Copy left singular vectors of A from A to U
823: *
824: CALL DLACPY( 'F', M, N, A, LDA, U, LDU )
825: *
826: END IF
827: *
828: ELSE
829: *
830: * M .LT. MNTHR
831: *
832: * Path 5 (M at least N, but not much larger)
833: * Reduce to bidiagonal form without QR decomposition
834: *
835: IE = 1
836: ITAUQ = IE + N
837: ITAUP = ITAUQ + N
838: NWORK = ITAUP + N
839: *
840: * Bidiagonalize A
841: * (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
842: *
843: CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
844: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
845: $ IERR )
846: IF( WNTQN ) THEN
847: *
848: * Perform bidiagonal SVD, only computing singular values
849: * (Workspace: need N+BDSPAC)
850: *
851: CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
852: $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
853: ELSE IF( WNTQO ) THEN
854: IU = NWORK
855: IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
856: *
857: * WORK( IU ) is M by N
858: *
859: LDWRKU = M
860: NWORK = IU + LDWRKU*N
861: CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),
862: $ LDWRKU )
863: ELSE
864: *
865: * WORK( IU ) is N by N
866: *
867: LDWRKU = N
868: NWORK = IU + LDWRKU*N
869: *
870: * WORK(IR) is LDWRKR by N
871: *
872: IR = NWORK
873: LDWRKR = ( LWORK-N*N-3*N ) / N
874: END IF
875: NWORK = IU + LDWRKU*N
876: *
877: * Perform bidiagonal SVD, computing left singular vectors
878: * of bidiagonal matrix in WORK(IU) and computing right
879: * singular vectors of bidiagonal matrix in VT
880: * (Workspace: need N+N*N+BDSPAC)
881: *
882: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
883: $ LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
884: $ IWORK, INFO )
885: *
886: * Overwrite VT by right singular vectors of A
887: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
888: *
889: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
890: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
891: $ LWORK-NWORK+1, IERR )
892: *
893: IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
894: *
895: * Overwrite WORK(IU) by left singular vectors of A
896: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
897: *
898: CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
899: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
900: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
901: *
902: * Copy left singular vectors of A from WORK(IU) to A
903: *
904: CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
905: ELSE
906: *
907: * Generate Q in A
908: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
909: *
910: CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
911: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
912: *
913: * Multiply Q in A by left singular vectors of
914: * bidiagonal matrix in WORK(IU), storing result in
915: * WORK(IR) and copying to A
916: * (Workspace: need 2*N*N, prefer N*N+M*N)
917: *
918: DO 20 I = 1, M, LDWRKR
919: CHUNK = MIN( M-I+1, LDWRKR )
920: CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
921: $ LDA, WORK( IU ), LDWRKU, ZERO,
922: $ WORK( IR ), LDWRKR )
923: CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
924: $ A( I, 1 ), LDA )
925: 20 CONTINUE
926: END IF
927: *
928: ELSE IF( WNTQS ) THEN
929: *
930: * Perform bidiagonal SVD, computing left singular vectors
931: * of bidiagonal matrix in U and computing right singular
932: * vectors of bidiagonal matrix in VT
933: * (Workspace: need N+BDSPAC)
934: *
935: CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU )
936: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
937: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
938: $ INFO )
939: *
940: * Overwrite U by left singular vectors of A and VT
941: * by right singular vectors of A
942: * (Workspace: need 3*N, prefer 2*N+N*NB)
943: *
944: CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
945: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
946: $ LWORK-NWORK+1, IERR )
947: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
948: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
949: $ LWORK-NWORK+1, IERR )
950: ELSE IF( WNTQA ) THEN
951: *
952: * Perform bidiagonal SVD, computing left singular vectors
953: * of bidiagonal matrix in U and computing right singular
954: * vectors of bidiagonal matrix in VT
955: * (Workspace: need N+BDSPAC)
956: *
957: CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU )
958: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
959: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
960: $ INFO )
961: *
962: * Set the right corner of U to identity matrix
963: *
964: IF( M.GT.N ) THEN
965: CALL DLASET( 'F', M-N, M-N, ZERO, ONE, U( N+1, N+1 ),
966: $ LDU )
967: END IF
968: *
969: * Overwrite U by left singular vectors of A and VT
970: * by right singular vectors of A
971: * (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB)
972: *
973: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
974: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
975: $ LWORK-NWORK+1, IERR )
976: CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
977: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
978: $ LWORK-NWORK+1, IERR )
979: END IF
980: *
981: END IF
982: *
983: ELSE
984: *
985: * A has more columns than rows. If A has sufficiently more
986: * columns than rows, first reduce using the LQ decomposition (if
987: * sufficient workspace available)
988: *
989: IF( N.GE.MNTHR ) THEN
990: *
991: IF( WNTQN ) THEN
992: *
993: * Path 1t (N much larger than M, JOBZ='N')
994: * No singular vectors to be computed
995: *
996: ITAU = 1
997: NWORK = ITAU + M
998: *
999: * Compute A=L*Q
1000: * (Workspace: need 2*M, prefer M+M*NB)
1001: *
1002: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1003: $ LWORK-NWORK+1, IERR )
1004: *
1005: * Zero out above L
1006: *
1007: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1008: IE = 1
1009: ITAUQ = IE + M
1010: ITAUP = ITAUQ + M
1011: NWORK = ITAUP + M
1012: *
1013: * Bidiagonalize L in A
1014: * (Workspace: need 4*M, prefer 3*M+2*M*NB)
1015: *
1016: CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1017: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1018: $ IERR )
1019: NWORK = IE + M
1020: *
1021: * Perform bidiagonal SVD, computing singular values only
1022: * (Workspace: need M+BDSPAC)
1023: *
1024: CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1025: $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1026: *
1027: ELSE IF( WNTQO ) THEN
1028: *
1029: * Path 2t (N much larger than M, JOBZ='O')
1030: * M right singular vectors to be overwritten on A and
1031: * M left singular vectors to be computed in U
1032: *
1033: IVT = 1
1034: *
1035: * IVT is M by M
1036: *
1037: IL = IVT + M*M
1038: IF( LWORK.GE.M*N+M*M+3*M+BDSPAC ) THEN
1039: *
1040: * WORK(IL) is M by N
1041: *
1042: LDWRKL = M
1043: CHUNK = N
1044: ELSE
1045: LDWRKL = M
1046: CHUNK = ( LWORK-M*M ) / M
1047: END IF
1048: ITAU = IL + LDWRKL*M
1049: NWORK = ITAU + M
1050: *
1051: * Compute A=L*Q
1052: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1053: *
1054: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1055: $ LWORK-NWORK+1, IERR )
1056: *
1057: * Copy L to WORK(IL), zeroing about above it
1058: *
1059: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1060: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
1061: $ WORK( IL+LDWRKL ), LDWRKL )
1062: *
1063: * Generate Q in A
1064: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1065: *
1066: CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1067: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1068: IE = ITAU
1069: ITAUQ = IE + M
1070: ITAUP = ITAUQ + M
1071: NWORK = ITAUP + M
1072: *
1073: * Bidiagonalize L in WORK(IL)
1074: * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1075: *
1076: CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1077: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1078: $ LWORK-NWORK+1, IERR )
1079: *
1080: * Perform bidiagonal SVD, computing left singular vectors
1081: * of bidiagonal matrix in U, and computing right singular
1082: * vectors of bidiagonal matrix in WORK(IVT)
1083: * (Workspace: need M+M*M+BDSPAC)
1084: *
1085: CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1086: $ WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
1087: $ IWORK, INFO )
1088: *
1089: * Overwrite U by left singular vectors of L and WORK(IVT)
1090: * by right singular vectors of L
1091: * (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
1092: *
1093: CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1094: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1095: $ LWORK-NWORK+1, IERR )
1096: CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1097: $ WORK( ITAUP ), WORK( IVT ), M,
1098: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1099: *
1100: * Multiply right singular vectors of L in WORK(IVT) by Q
1101: * in A, storing result in WORK(IL) and copying to A
1102: * (Workspace: need 2*M*M, prefer M*M+M*N)
1103: *
1104: DO 30 I = 1, N, CHUNK
1105: BLK = MIN( N-I+1, CHUNK )
1106: CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,
1107: $ A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
1108: CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
1109: $ A( 1, I ), LDA )
1110: 30 CONTINUE
1111: *
1112: ELSE IF( WNTQS ) THEN
1113: *
1114: * Path 3t (N much larger than M, JOBZ='S')
1115: * M right singular vectors to be computed in VT and
1116: * M left singular vectors to be computed in U
1117: *
1118: IL = 1
1119: *
1120: * WORK(IL) is M by M
1121: *
1122: LDWRKL = M
1123: ITAU = IL + LDWRKL*M
1124: NWORK = ITAU + M
1125: *
1126: * Compute A=L*Q
1127: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1128: *
1129: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1130: $ LWORK-NWORK+1, IERR )
1131: *
1132: * Copy L to WORK(IL), zeroing out above it
1133: *
1134: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1135: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
1136: $ WORK( IL+LDWRKL ), LDWRKL )
1137: *
1138: * Generate Q in A
1139: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1140: *
1141: CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1142: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1143: IE = ITAU
1144: ITAUQ = IE + M
1145: ITAUP = ITAUQ + M
1146: NWORK = ITAUP + M
1147: *
1148: * Bidiagonalize L in WORK(IU), copying result to U
1149: * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1150: *
1151: CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1152: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1153: $ LWORK-NWORK+1, IERR )
1154: *
1155: * Perform bidiagonal SVD, computing left singular vectors
1156: * of bidiagonal matrix in U and computing right singular
1157: * vectors of bidiagonal matrix in VT
1158: * (Workspace: need M+BDSPAC)
1159: *
1160: CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
1161: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1162: $ INFO )
1163: *
1164: * Overwrite U by left singular vectors of L and VT
1165: * by right singular vectors of L
1166: * (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1167: *
1168: CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1169: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1170: $ LWORK-NWORK+1, IERR )
1171: CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1172: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1173: $ LWORK-NWORK+1, IERR )
1174: *
1175: * Multiply right singular vectors of L in WORK(IL) by
1176: * Q in A, storing result in VT
1177: * (Workspace: need M*M)
1178: *
1179: CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1180: CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,
1181: $ A, LDA, ZERO, VT, LDVT )
1182: *
1183: ELSE IF( WNTQA ) THEN
1184: *
1185: * Path 4t (N much larger than M, JOBZ='A')
1186: * N right singular vectors to be computed in VT and
1187: * M left singular vectors to be computed in U
1188: *
1189: IVT = 1
1190: *
1191: * WORK(IVT) is M by M
1192: *
1193: LDWKVT = M
1194: ITAU = IVT + LDWKVT*M
1195: NWORK = ITAU + M
1196: *
1197: * Compute A=L*Q, copying result to VT
1198: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1199: *
1200: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1201: $ LWORK-NWORK+1, IERR )
1202: CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
1203: *
1204: * Generate Q in VT
1205: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1206: *
1207: CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1208: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1209: *
1210: * Produce L in A, zeroing out other entries
1211: *
1212: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1213: IE = ITAU
1214: ITAUQ = IE + M
1215: ITAUP = ITAUQ + M
1216: NWORK = ITAUP + M
1217: *
1218: * Bidiagonalize L in A
1219: * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1220: *
1221: CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1222: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1223: $ IERR )
1224: *
1225: * Perform bidiagonal SVD, computing left singular vectors
1226: * of bidiagonal matrix in U and computing right singular
1227: * vectors of bidiagonal matrix in WORK(IVT)
1228: * (Workspace: need M+M*M+BDSPAC)
1229: *
1230: CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1231: $ WORK( IVT ), LDWKVT, DUM, IDUM,
1232: $ WORK( NWORK ), IWORK, INFO )
1233: *
1234: * Overwrite U by left singular vectors of L and WORK(IVT)
1235: * by right singular vectors of L
1236: * (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1237: *
1238: CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1239: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1240: $ LWORK-NWORK+1, IERR )
1241: CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
1242: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1243: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1244: *
1245: * Multiply right singular vectors of L in WORK(IVT) by
1246: * Q in VT, storing result in A
1247: * (Workspace: need M*M)
1248: *
1249: CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,
1250: $ VT, LDVT, ZERO, A, LDA )
1251: *
1252: * Copy right singular vectors of A from A to VT
1253: *
1254: CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT )
1255: *
1256: END IF
1257: *
1258: ELSE
1259: *
1260: * N .LT. MNTHR
1261: *
1262: * Path 5t (N greater than M, but not much larger)
1263: * Reduce to bidiagonal form without LQ decomposition
1264: *
1265: IE = 1
1266: ITAUQ = IE + M
1267: ITAUP = ITAUQ + M
1268: NWORK = ITAUP + M
1269: *
1270: * Bidiagonalize A
1271: * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
1272: *
1273: CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1274: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1275: $ IERR )
1276: IF( WNTQN ) THEN
1277: *
1278: * Perform bidiagonal SVD, only computing singular values
1279: * (Workspace: need M+BDSPAC)
1280: *
1281: CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1282: $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1283: ELSE IF( WNTQO ) THEN
1284: LDWKVT = M
1285: IVT = NWORK
1286: IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
1287: *
1288: * WORK( IVT ) is M by N
1289: *
1290: CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),
1291: $ LDWKVT )
1292: NWORK = IVT + LDWKVT*N
1293: ELSE
1294: *
1295: * WORK( IVT ) is M by M
1296: *
1297: NWORK = IVT + LDWKVT*M
1298: IL = NWORK
1299: *
1300: * WORK(IL) is M by CHUNK
1301: *
1302: CHUNK = ( LWORK-M*M-3*M ) / M
1303: END IF
1304: *
1305: * Perform bidiagonal SVD, computing left singular vectors
1306: * of bidiagonal matrix in U and computing right singular
1307: * vectors of bidiagonal matrix in WORK(IVT)
1308: * (Workspace: need M*M+BDSPAC)
1309: *
1310: CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,
1311: $ WORK( IVT ), LDWKVT, DUM, IDUM,
1312: $ WORK( NWORK ), IWORK, INFO )
1313: *
1314: * Overwrite U by left singular vectors of A
1315: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1316: *
1317: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1318: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1319: $ LWORK-NWORK+1, IERR )
1320: *
1321: IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
1322: *
1323: * Overwrite WORK(IVT) by left singular vectors of A
1324: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1325: *
1326: CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1327: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1328: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1329: *
1330: * Copy right singular vectors of A from WORK(IVT) to A
1331: *
1332: CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
1333: ELSE
1334: *
1335: * Generate P**T in A
1336: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1337: *
1338: CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1339: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
1340: *
1341: * Multiply Q in A by right singular vectors of
1342: * bidiagonal matrix in WORK(IVT), storing result in
1343: * WORK(IL) and copying to A
1344: * (Workspace: need 2*M*M, prefer M*M+M*N)
1345: *
1346: DO 40 I = 1, N, CHUNK
1347: BLK = MIN( N-I+1, CHUNK )
1348: CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),
1349: $ LDWKVT, A( 1, I ), LDA, ZERO,
1350: $ WORK( IL ), M )
1351: CALL DLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ),
1352: $ LDA )
1353: 40 CONTINUE
1354: END IF
1355: ELSE IF( WNTQS ) THEN
1356: *
1357: * Perform bidiagonal SVD, computing left singular vectors
1358: * of bidiagonal matrix in U and computing right singular
1359: * vectors of bidiagonal matrix in VT
1360: * (Workspace: need M+BDSPAC)
1361: *
1362: CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
1363: CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1364: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1365: $ INFO )
1366: *
1367: * Overwrite U by left singular vectors of A and VT
1368: * by right singular vectors of A
1369: * (Workspace: need 3*M, prefer 2*M+M*NB)
1370: *
1371: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1372: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1373: $ LWORK-NWORK+1, IERR )
1374: CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1375: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1376: $ LWORK-NWORK+1, IERR )
1377: ELSE IF( WNTQA ) THEN
1378: *
1379: * Perform bidiagonal SVD, computing left singular vectors
1380: * of bidiagonal matrix in U and computing right singular
1381: * vectors of bidiagonal matrix in VT
1382: * (Workspace: need M+BDSPAC)
1383: *
1384: CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
1385: CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1386: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1387: $ INFO )
1388: *
1389: * Set the right corner of VT to identity matrix
1390: *
1391: IF( N.GT.M ) THEN
1392: CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT( M+1, M+1 ),
1393: $ LDVT )
1394: END IF
1395: *
1396: * Overwrite U by left singular vectors of A and VT
1397: * by right singular vectors of A
1398: * (Workspace: need 2*M+N, prefer 2*M+N*NB)
1399: *
1400: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1401: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1402: $ LWORK-NWORK+1, IERR )
1403: CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1404: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1405: $ LWORK-NWORK+1, IERR )
1406: END IF
1407: *
1408: END IF
1409: *
1410: END IF
1411: *
1412: * Undo scaling if necessary
1413: *
1414: IF( ISCL.EQ.1 ) THEN
1415: IF( ANRM.GT.BIGNUM )
1416: $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
1417: $ IERR )
1418: IF( ANRM.LT.SMLNUM )
1419: $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
1420: $ IERR )
1421: END IF
1422: *
1423: * Return optimal workspace in WORK(1)
1424: *
1425: WORK( 1 ) = MAXWRK
1426: *
1427: RETURN
1428: *
1429: * End of DGESDD
1430: *
1431: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgesdd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack