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