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