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,
22: * WORK, 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;
158: *> if 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 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.
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: *
202: *> \author Univ. of Tennessee
203: *> \author Univ. of California Berkeley
204: *> \author Univ. of Colorado Denver
205: *> \author NAG Ltd.
206: *
207: *> \date June 2016
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: * =====================================================================
218: SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
219: $ WORK, LWORK, IWORK, INFO )
220: implicit none
221: *
222: * -- LAPACK driver routine (version 3.7.0) --
223: * -- LAPACK is a software package provided by Univ. of Tennessee, --
224: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
225: * June 2016
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
249: INTEGER LWORK_DGEBRD_MN, LWORK_DGEBRD_MM,
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
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
272: EXTERNAL DLAMCH, DLANGE, LSAME
273: * ..
274: * .. Intrinsic Functions ..
275: INTRINSIC INT, MAX, MIN, SQRT
276: * ..
277: * .. Executable Statements ..
278: *
279: * Test the input arguments
280: *
281: INFO = 0
282: MINMN = MIN( M, N )
283: WNTQA = LSAME( JOBZ, 'A' )
284: WNTQS = LSAME( JOBZ, 'S' )
285: WNTQAS = WNTQA .OR. WNTQS
286: WNTQO = LSAME( JOBZ, 'O' )
287: WNTQN = LSAME( JOBZ, 'N' )
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
308: * Note: Comments in the code beginning "Workspace:" describe the
309: * minimal amount of workspace allocated at that point in the code,
310: * as well as the preferred amount for good performance.
311: * NB refers to the optimal block size for the immediately
312: * following subroutine, as returned by ILAENV.
313: *
314: IF( INFO.EQ.0 ) THEN
315: MINWRK = 1
316: MAXWRK = 1
317: BDSPAC = 0
318: MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
319: IF( M.GE.N .AND. MINMN.GT.0 ) THEN
320: *
321: * Compute space needed for DBDSDC
322: *
323: IF( WNTQN ) THEN
324: * dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
325: * keep 7*N for backwards compatability.
326: BDSPAC = 7*N
327: ELSE
328: BDSPAC = 3*N*N + 4*N
329: END IF
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: *
369: IF( M.GE.MNTHR ) THEN
370: IF( WNTQN ) THEN
371: *
372: * Path 1 (M >> N, JOBZ='N')
373: *
374: WRKBL = N + LWORK_DGEQRF_MN
375: WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
376: MAXWRK = MAX( WRKBL, BDSPAC + N )
377: MINWRK = BDSPAC + N
378: ELSE IF( WNTQO ) THEN
379: *
380: * Path 2 (M >> N, JOBZ='O')
381: *
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 )
388: MAXWRK = WRKBL + 2*N*N
389: MINWRK = BDSPAC + 2*N*N + 3*N
390: ELSE IF( WNTQS ) THEN
391: *
392: * Path 3 (M >> N, JOBZ='S')
393: *
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 )
400: MAXWRK = WRKBL + N*N
401: MINWRK = BDSPAC + N*N + 3*N
402: ELSE IF( WNTQA ) THEN
403: *
404: * Path 4 (M >> N, JOBZ='A')
405: *
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 )
412: MAXWRK = WRKBL + N*N
413: MINWRK = N*N + MAX( 3*N + BDSPAC, N + M )
414: END IF
415: ELSE
416: *
417: * Path 5 (M >= N, but not much larger)
418: *
419: WRKBL = 3*N + LWORK_DGEBRD_MN
420: IF( WNTQN ) THEN
421: * Path 5n (M >= N, jobz='N')
422: MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
423: MINWRK = 3*N + MAX( M, BDSPAC )
424: ELSE IF( WNTQO ) THEN
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 )
429: MAXWRK = WRKBL + M*N
430: MINWRK = 3*N + MAX( M, N*N + BDSPAC )
431: ELSE IF( WNTQS ) THEN
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 )
436: MINWRK = 3*N + MAX( M, BDSPAC )
437: ELSE IF( WNTQA ) THEN
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 )
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
450: * dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
451: * keep 7*N for backwards compatability.
452: BDSPAC = 7*M
453: ELSE
454: BDSPAC = 3*M*M + 4*M
455: END IF
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: *
494: IF( N.GE.MNTHR ) THEN
495: IF( WNTQN ) THEN
496: *
497: * Path 1t (N >> M, JOBZ='N')
498: *
499: WRKBL = M + LWORK_DGELQF_MN
500: WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
501: MAXWRK = MAX( WRKBL, BDSPAC + M )
502: MINWRK = BDSPAC + M
503: ELSE IF( WNTQO ) THEN
504: *
505: * Path 2t (N >> M, JOBZ='O')
506: *
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 )
513: MAXWRK = WRKBL + 2*M*M
514: MINWRK = BDSPAC + 2*M*M + 3*M
515: ELSE IF( WNTQS ) THEN
516: *
517: * Path 3t (N >> M, JOBZ='S')
518: *
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 )
525: MAXWRK = WRKBL + M*M
526: MINWRK = BDSPAC + M*M + 3*M
527: ELSE IF( WNTQA ) THEN
528: *
529: * Path 4t (N >> M, JOBZ='A')
530: *
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 )
537: MAXWRK = WRKBL + M*M
538: MINWRK = M*M + MAX( 3*M + BDSPAC, M + N )
539: END IF
540: ELSE
541: *
542: * Path 5t (N > M, but not much larger)
543: *
544: WRKBL = 3*M + LWORK_DGEBRD_MN
545: IF( WNTQN ) THEN
546: * Path 5tn (N > M, jobz='N')
547: MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
548: MINWRK = 3*M + MAX( N, BDSPAC )
549: ELSE IF( WNTQO ) THEN
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 )
554: MAXWRK = WRKBL + M*N
555: MINWRK = 3*M + MAX( N, M*M + BDSPAC )
556: ELSE IF( WNTQS ) THEN
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 )
561: MINWRK = 3*M + MAX( N, BDSPAC )
562: ELSE IF( WNTQA ) THEN
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 )
567: MINWRK = 3*M + MAX( N, BDSPAC )
568: END IF
569: END IF
570: END IF
571:
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: *
621: * Path 1 (M >> N, JOBZ='N')
622: * No singular vectors to be computed
623: *
624: ITAU = 1
625: NWORK = ITAU + N
626: *
627: * Compute A=Q*R
628: * Workspace: need N [tau] + N [work]
629: * Workspace: prefer N [tau] + N*NB [work]
630: *
631: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
632: $ LWORK - NWORK + 1, IERR )
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
643: * Workspace: need 3*N [e, tauq, taup] + N [work]
644: * Workspace: prefer 3*N [e, tauq, taup] + 2*N*NB [work]
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
652: * Workspace: need N [e] + BDSPAC
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: *
659: * Path 2 (M >> N, JOBZ = 'O')
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: *
667: IF( LWORK .GE. LDA*N + N*N + 3*N + BDSPAC ) THEN
668: LDWRKR = LDA
669: ELSE
670: LDWRKR = ( LWORK - N*N - 3*N - BDSPAC ) / N
671: END IF
672: ITAU = IR + LDWRKR*N
673: NWORK = ITAU + N
674: *
675: * Compute A=Q*R
676: * Workspace: need N*N [R] + N [tau] + N [work]
677: * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
678: *
679: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
680: $ LWORK - NWORK + 1, IERR )
681: *
682: * Copy R to WORK(IR), zeroing out below it
683: *
684: CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
685: CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
686: $ LDWRKR )
687: *
688: * Generate Q in A
689: * Workspace: need N*N [R] + N [tau] + N [work]
690: * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
691: *
692: CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
693: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
694: IE = ITAU
695: ITAUQ = IE + N
696: ITAUP = ITAUQ + N
697: NWORK = ITAUP + N
698: *
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]
702: *
703: CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
704: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
705: $ LWORK - NWORK + 1, IERR )
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
715: * Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U] + BDSPAC
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
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]
725: *
726: CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
727: $ WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
728: $ LWORK - NWORK + 1, IERR )
729: CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
730: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
731: $ LWORK - NWORK + 1, IERR )
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
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]
737: *
738: DO 10 I = 1, M, LDWRKR
739: CHUNK = MIN( M - I + 1, LDWRKR )
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: *
749: * Path 3 (M >> N, JOBZ='S')
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
762: * Workspace: need N*N [R] + N [tau] + N [work]
763: * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
764: *
765: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
766: $ LWORK - NWORK + 1, IERR )
767: *
768: * Copy R to WORK(IR), zeroing out below it
769: *
770: CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
771: CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
772: $ LDWRKR )
773: *
774: * Generate Q in A
775: * Workspace: need N*N [R] + N [tau] + N [work]
776: * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
777: *
778: CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
779: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
780: IE = ITAU
781: ITAUQ = IE + N
782: ITAUP = ITAUQ + N
783: NWORK = ITAUP + N
784: *
785: * Bidiagonalize R in WORK(IR)
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]
788: *
789: CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
790: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
791: $ LWORK - NWORK + 1, IERR )
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
796: * Workspace: need N*N [R] + 3*N [e, tauq, taup] + BDSPAC
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
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]
806: *
807: CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
808: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
809: $ LWORK - NWORK + 1, IERR )
810: *
811: CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
812: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
813: $ LWORK - NWORK + 1, IERR )
814: *
815: * Multiply Q in A by left singular vectors of R in
816: * WORK(IR), storing result in U
817: * Workspace: need N*N [R]
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: *
825: * Path 4 (M >> N, JOBZ='A')
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
838: * Workspace: need N*N [U] + N [tau] + N [work]
839: * Workspace: prefer N*N [U] + N [tau] + N*NB [work]
840: *
841: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
842: $ LWORK - NWORK + 1, IERR )
843: CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
844: *
845: * Generate Q in U
846: * Workspace: need N*N [U] + N [tau] + M [work]
847: * Workspace: prefer N*N [U] + N [tau] + M*NB [work]
848: CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
849: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
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
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]
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
870: * Workspace: need N*N [U] + 3*N [e, tauq, taup] + BDSPAC
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
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]
880: *
881: CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
882: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
883: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
884: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
885: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
886: $ LWORK - NWORK + 1, IERR )
887: *
888: * Multiply Q in U by left singular vectors of R in
889: * WORK(IU), storing result in A
890: * Workspace: need N*N [U]
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: *
905: * Path 5 (M >= N, but not much larger)
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
914: * Workspace: need 3*N [e, tauq, taup] + M [work]
915: * Workspace: prefer 3*N [e, tauq, taup] + (M+N)*NB [work]
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: *
922: * Path 5n (M >= N, JOBZ='N')
923: * Perform bidiagonal SVD, only computing singular values
924: * Workspace: need 3*N [e, tauq, taup] + BDSPAC
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
929: * Path 5o (M >= N, JOBZ='O')
930: IU = NWORK
931: IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN
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 )
939: * IR is unused; silence compile warnings
940: IR = -1
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
951: LDWRKR = ( LWORK - N*N - 3*N ) / N
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
958: * Workspace: need 3*N [e, tauq, taup] + N*N [U] + BDSPAC
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
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]
967: *
968: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
969: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
970: $ LWORK - NWORK + 1, IERR )
971: *
972: IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN
973: *
974: * Path 5o-fast
975: * Overwrite WORK(IU) by left singular vectors of A
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]
978: *
979: CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
980: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
981: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
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: *
988: * Path 5o-slow
989: * Generate Q in A
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]
992: *
993: CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
994: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
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
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]
1001: *
1002: DO 20 I = 1, M, LDWRKR
1003: CHUNK = MIN( M - I + 1, LDWRKR )
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: *
1014: * Path 5s (M >= N, JOBZ='S')
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
1018: * Workspace: need 3*N [e, tauq, taup] + BDSPAC
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
1027: * Workspace: need 3*N [e, tauq, taup] + N [work]
1028: * Workspace: prefer 3*N [e, tauq, taup] + N*NB [work]
1029: *
1030: CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1031: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1032: $ LWORK - NWORK + 1, IERR )
1033: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
1034: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1035: $ LWORK - NWORK + 1, IERR )
1036: ELSE IF( WNTQA ) THEN
1037: *
1038: * Path 5a (M >= N, JOBZ='A')
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
1042: * Workspace: need 3*N [e, tauq, taup] + BDSPAC
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
1052: CALL DLASET( 'F', M - N, M - N, ZERO, ONE, U(N+1,N+1),
1053: $ LDU )
1054: END IF
1055: *
1056: * Overwrite U by left singular vectors of A and VT
1057: * by right singular vectors of A
1058: * Workspace: need 3*N [e, tauq, taup] + M [work]
1059: * Workspace: prefer 3*N [e, tauq, taup] + M*NB [work]
1060: *
1061: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1062: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1063: $ LWORK - NWORK + 1, IERR )
1064: CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1065: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1066: $ LWORK - NWORK + 1, IERR )
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: *
1081: * Path 1t (N >> M, JOBZ='N')
1082: * No singular vectors to be computed
1083: *
1084: ITAU = 1
1085: NWORK = ITAU + M
1086: *
1087: * Compute A=L*Q
1088: * Workspace: need M [tau] + M [work]
1089: * Workspace: prefer M [tau] + M*NB [work]
1090: *
1091: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1092: $ LWORK - NWORK + 1, IERR )
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
1103: * Workspace: need 3*M [e, tauq, taup] + M [work]
1104: * Workspace: prefer 3*M [e, tauq, taup] + 2*M*NB [work]
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
1112: * Workspace: need M [e] + BDSPAC
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: *
1119: * Path 2t (N >> M, JOBZ='O')
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: *
1125: * WORK(IVT) is M by M
1126: * WORK(IL) is M by M; it is later resized to M by chunk for gemm
1127: *
1128: IL = IVT + M*M
1129: IF( LWORK .GE. M*N + M*M + 3*M + BDSPAC ) THEN
1130: LDWRKL = M
1131: CHUNK = N
1132: ELSE
1133: LDWRKL = M
1134: CHUNK = ( LWORK - M*M ) / M
1135: END IF
1136: ITAU = IL + LDWRKL*M
1137: NWORK = ITAU + M
1138: *
1139: * Compute A=L*Q
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]
1142: *
1143: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1144: $ LWORK - NWORK + 1, IERR )
1145: *
1146: * Copy L to WORK(IL), zeroing about above it
1147: *
1148: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1149: CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
1150: $ WORK( IL + LDWRKL ), LDWRKL )
1151: *
1152: * Generate Q in A
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]
1155: *
1156: CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1157: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1158: IE = ITAU
1159: ITAUQ = IE + M
1160: ITAUP = ITAUQ + M
1161: NWORK = ITAUP + M
1162: *
1163: * Bidiagonalize L in WORK(IL)
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]
1166: *
1167: CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1168: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1169: $ LWORK - NWORK + 1, IERR )
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)
1174: * Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + BDSPAC
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
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]
1184: *
1185: CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1186: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1187: $ LWORK - NWORK + 1, IERR )
1188: CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1189: $ WORK( ITAUP ), WORK( IVT ), M,
1190: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
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
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.
1197: *
1198: DO 30 I = 1, N, CHUNK
1199: BLK = MIN( N - I + 1, CHUNK )
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: *
1208: * Path 3t (N >> M, JOBZ='S')
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
1221: * Workspace: need M*M [L] + M [tau] + M [work]
1222: * Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1223: *
1224: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1225: $ LWORK - NWORK + 1, IERR )
1226: *
1227: * Copy L to WORK(IL), zeroing out above it
1228: *
1229: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1230: CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
1231: $ WORK( IL + LDWRKL ), LDWRKL )
1232: *
1233: * Generate Q in A
1234: * Workspace: need M*M [L] + M [tau] + M [work]
1235: * Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1236: *
1237: CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1238: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1239: IE = ITAU
1240: ITAUQ = IE + M
1241: ITAUP = ITAUQ + M
1242: NWORK = ITAUP + M
1243: *
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]
1247: *
1248: CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1249: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1250: $ LWORK - NWORK + 1, IERR )
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
1255: * Workspace: need M*M [L] + 3*M [e, tauq, taup] + BDSPAC
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
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]
1265: *
1266: CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1267: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1268: $ LWORK - NWORK + 1, IERR )
1269: CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1270: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1271: $ LWORK - NWORK + 1, IERR )
1272: *
1273: * Multiply right singular vectors of L in WORK(IL) by
1274: * Q in A, storing result in VT
1275: * Workspace: need M*M [L]
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: *
1283: * Path 4t (N >> M, JOBZ='A')
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
1296: * Workspace: need M*M [VT] + M [tau] + M [work]
1297: * Workspace: prefer M*M [VT] + M [tau] + M*NB [work]
1298: *
1299: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1300: $ LWORK - NWORK + 1, IERR )
1301: CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
1302: *
1303: * Generate Q in VT
1304: * Workspace: need M*M [VT] + M [tau] + N [work]
1305: * Workspace: prefer M*M [VT] + M [tau] + N*NB [work]
1306: *
1307: CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1308: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
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
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]
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)
1329: * Workspace: need M*M [VT] + 3*M [e, tauq, taup] + BDSPAC
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
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]
1339: *
1340: CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1341: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1342: $ LWORK - NWORK + 1, IERR )
1343: CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
1344: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1345: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1346: *
1347: * Multiply right singular vectors of L in WORK(IVT) by
1348: * Q in VT, storing result in A
1349: * Workspace: need M*M [VT]
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: *
1364: * Path 5t (N > M, but not much larger)
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
1373: * Workspace: need 3*M [e, tauq, taup] + N [work]
1374: * Workspace: prefer 3*M [e, tauq, taup] + (M+N)*NB [work]
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: *
1381: * Path 5tn (N > M, JOBZ='N')
1382: * Perform bidiagonal SVD, only computing singular values
1383: * Workspace: need 3*M [e, tauq, taup] + BDSPAC
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
1388: * Path 5to (N > M, JOBZ='O')
1389: LDWKVT = M
1390: IVT = NWORK
1391: IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN
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
1398: * IL is unused; silence compile warnings
1399: IL = -1
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: *
1409: CHUNK = ( LWORK - M*M - 3*M ) / M
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)
1415: * Workspace: need 3*M [e, tauq, taup] + M*M [VT] + BDSPAC
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
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]
1424: *
1425: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1426: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1427: $ LWORK - NWORK + 1, IERR )
1428: *
1429: IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN
1430: *
1431: * Path 5to-fast
1432: * Overwrite WORK(IVT) by left singular vectors of A
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]
1435: *
1436: CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1437: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1438: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
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: *
1445: * Path 5to-slow
1446: * Generate P**T in A
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]
1449: *
1450: CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1451: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
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
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]
1458: *
1459: DO 40 I = 1, N, CHUNK
1460: BLK = MIN( N - I + 1, CHUNK )
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: *
1470: * Path 5ts (N > M, JOBZ='S')
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
1474: * Workspace: need 3*M [e, tauq, taup] + BDSPAC
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
1483: * Workspace: need 3*M [e, tauq, taup] + M [work]
1484: * Workspace: prefer 3*M [e, tauq, taup] + M*NB [work]
1485: *
1486: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1487: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1488: $ LWORK - NWORK + 1, IERR )
1489: CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1490: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1491: $ LWORK - NWORK + 1, IERR )
1492: ELSE IF( WNTQA ) THEN
1493: *
1494: * Path 5ta (N > M, JOBZ='A')
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
1498: * Workspace: need 3*M [e, tauq, taup] + BDSPAC
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
1508: CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT(M+1,M+1),
1509: $ LDVT )
1510: END IF
1511: *
1512: * Overwrite U by left singular vectors of A and VT
1513: * by right singular vectors of A
1514: * Workspace: need 3*M [e, tauq, taup] + N [work]
1515: * Workspace: prefer 3*M [e, tauq, taup] + N*NB [work]
1516: *
1517: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1518: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1519: $ LWORK - NWORK + 1, IERR )
1520: CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1521: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1522: $ LWORK - NWORK + 1, IERR )
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