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