Annotation of rpl/lapack/lapack/dgesdd.f, revision 1.16
1.8 bertrand 1: *> \brief \b DGESDD
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DGESDD + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesdd.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesdd.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesdd.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
1.16 ! bertrand 21: * SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
! 22: * WORK, LWORK, IWORK, INFO )
1.8 bertrand 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
1.16 ! bertrand 157: *> The leading dimension of the array VT. LDVT >= 1;
! 158: *> if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
1.8 bertrand 159: *> if JOBZ = 'S', LDVT >= min(M,N).
160: *> \endverbatim
161: *>
162: *> \param[out] WORK
163: *> \verbatim
164: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
165: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
166: *> \endverbatim
167: *>
168: *> \param[in] LWORK
169: *> \verbatim
170: *> LWORK is INTEGER
171: *> The dimension of the array WORK. LWORK >= 1.
1.16 ! bertrand 172: *> If LWORK = -1, a workspace query is assumed. The optimal
! 173: *> size for the WORK array is calculated and stored in WORK(1),
! 174: *> and no other work except argument checking is performed.
! 175: *>
! 176: *> Let mx = max(M,N) and mn = min(M,N).
! 177: *> If JOBZ = 'N', LWORK >= 3*mn + max( mx, 7*mn ).
! 178: *> If JOBZ = 'O', LWORK >= 3*mn + max( mx, 5*mn*mn + 4*mn ).
! 179: *> If JOBZ = 'S', LWORK >= 4*mn*mn + 7*mn.
! 180: *> If JOBZ = 'A', LWORK >= 4*mn*mn + 6*mn + mx.
! 181: *> These are not tight minimums in all cases; see comments inside code.
! 182: *> For good performance, LWORK should generally be larger;
! 183: *> a query is recommended.
1.8 bertrand 184: *> \endverbatim
185: *>
186: *> \param[out] IWORK
187: *> \verbatim
188: *> IWORK is INTEGER array, dimension (8*min(M,N))
189: *> \endverbatim
190: *>
191: *> \param[out] INFO
192: *> \verbatim
193: *> INFO is INTEGER
194: *> = 0: successful exit.
195: *> < 0: if INFO = -i, the i-th argument had an illegal value.
196: *> > 0: DBDSDC did not converge, updating process failed.
197: *> \endverbatim
198: *
199: * Authors:
200: * ========
201: *
202: *> \author Univ. of Tennessee
203: *> \author Univ. of California Berkeley
204: *> \author Univ. of Colorado Denver
205: *> \author NAG Ltd.
206: *
1.16 ! bertrand 207: *> \date June 2016
1.8 bertrand 208: *
209: *> \ingroup doubleGEsing
210: *
211: *> \par Contributors:
212: * ==================
213: *>
214: *> Ming Gu and Huan Ren, Computer Science Division, University of
215: *> California at Berkeley, USA
216: *>
1.16 ! bertrand 217: *> @precisions fortran d -> s
1.8 bertrand 218: * =====================================================================
1.16 ! bertrand 219: SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
! 220: $ WORK, LWORK, IWORK, INFO )
! 221: implicit none
1.1 bertrand 222: *
1.16 ! bertrand 223: * -- LAPACK driver routine (version 3.6.1) --
1.1 bertrand 224: * -- LAPACK is a software package provided by Univ. of Tennessee, --
225: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.16 ! bertrand 226: * June 2016
1.1 bertrand 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
1.16 ! bertrand 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
1.1 bertrand 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
1.16 ! bertrand 273: EXTERNAL DLAMCH, DLANGE, LSAME
1.1 bertrand 274: * ..
275: * .. Intrinsic Functions ..
276: INTRINSIC INT, MAX, MIN, SQRT
277: * ..
278: * .. Executable Statements ..
279: *
280: * Test the input arguments
281: *
1.16 ! bertrand 282: INFO = 0
! 283: MINMN = MIN( M, N )
! 284: WNTQA = LSAME( JOBZ, 'A' )
! 285: WNTQS = LSAME( JOBZ, 'S' )
1.1 bertrand 286: WNTQAS = WNTQA .OR. WNTQS
1.16 ! bertrand 287: WNTQO = LSAME( JOBZ, 'O' )
! 288: WNTQN = LSAME( JOBZ, 'N' )
1.1 bertrand 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
1.16 ! bertrand 309: * Note: Comments in the code beginning "Workspace:" describe the
! 310: * minimal amount of workspace allocated at that point in the code,
1.1 bertrand 311: * as well as the preferred amount for good performance.
312: * NB refers to the optimal block size for the immediately
1.16 ! bertrand 313: * following subroutine, as returned by ILAENV.
1.1 bertrand 314: *
315: IF( INFO.EQ.0 ) THEN
316: MINWRK = 1
317: MAXWRK = 1
1.16 ! bertrand 318: BDSPAC = 0
! 319: MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
1.1 bertrand 320: IF( M.GE.N .AND. MINMN.GT.0 ) THEN
321: *
322: * Compute space needed for DBDSDC
323: *
324: IF( WNTQN ) THEN
1.16 ! bertrand 325: * dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
! 326: * keep 7*N for backwards compatability.
1.1 bertrand 327: BDSPAC = 7*N
328: ELSE
329: BDSPAC = 3*N*N + 4*N
330: END IF
1.16 ! bertrand 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: *
1.1 bertrand 370: IF( M.GE.MNTHR ) THEN
371: IF( WNTQN ) THEN
372: *
1.16 ! bertrand 373: * Path 1 (M >> N, JOBZ='N')
1.1 bertrand 374: *
1.16 ! bertrand 375: WRKBL = N + LWORK_DGEQRF_MN
! 376: WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN )
! 377: MAXWRK = MAX( WRKBL, BDSPAC + N )
1.1 bertrand 378: MINWRK = BDSPAC + N
379: ELSE IF( WNTQO ) THEN
380: *
1.16 ! bertrand 381: * Path 2 (M >> N, JOBZ='O')
1.1 bertrand 382: *
1.16 ! bertrand 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 )
1.1 bertrand 389: MAXWRK = WRKBL + 2*N*N
390: MINWRK = BDSPAC + 2*N*N + 3*N
391: ELSE IF( WNTQS ) THEN
392: *
1.16 ! bertrand 393: * Path 3 (M >> N, JOBZ='S')
1.1 bertrand 394: *
1.16 ! bertrand 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 )
1.1 bertrand 401: MAXWRK = WRKBL + N*N
402: MINWRK = BDSPAC + N*N + 3*N
403: ELSE IF( WNTQA ) THEN
404: *
1.16 ! bertrand 405: * Path 4 (M >> N, JOBZ='A')
1.1 bertrand 406: *
1.16 ! bertrand 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 )
1.1 bertrand 413: MAXWRK = WRKBL + N*N
1.16 ! bertrand 414: MINWRK = N*N + MAX( 3*N + BDSPAC, N + M )
1.1 bertrand 415: END IF
416: ELSE
417: *
1.16 ! bertrand 418: * Path 5 (M >= N, but not much larger)
1.1 bertrand 419: *
1.16 ! bertrand 420: WRKBL = 3*N + LWORK_DGEBRD_MN
1.1 bertrand 421: IF( WNTQN ) THEN
1.16 ! bertrand 422: * Path 5n (M >= N, jobz='N')
! 423: MAXWRK = MAX( WRKBL, 3*N + BDSPAC )
1.1 bertrand 424: MINWRK = 3*N + MAX( M, BDSPAC )
425: ELSE IF( WNTQO ) THEN
1.16 ! bertrand 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 )
1.1 bertrand 430: MAXWRK = WRKBL + M*N
1.16 ! bertrand 431: MINWRK = 3*N + MAX( M, N*N + BDSPAC )
1.1 bertrand 432: ELSE IF( WNTQS ) THEN
1.16 ! bertrand 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 )
1.1 bertrand 437: MINWRK = 3*N + MAX( M, BDSPAC )
438: ELSE IF( WNTQA ) THEN
1.16 ! bertrand 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 )
1.1 bertrand 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
1.16 ! bertrand 451: * dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
! 452: * keep 7*N for backwards compatability.
1.1 bertrand 453: BDSPAC = 7*M
454: ELSE
455: BDSPAC = 3*M*M + 4*M
456: END IF
1.16 ! bertrand 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: *
1.1 bertrand 495: IF( N.GE.MNTHR ) THEN
496: IF( WNTQN ) THEN
497: *
1.16 ! bertrand 498: * Path 1t (N >> M, JOBZ='N')
1.1 bertrand 499: *
1.16 ! bertrand 500: WRKBL = M + LWORK_DGELQF_MN
! 501: WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM )
! 502: MAXWRK = MAX( WRKBL, BDSPAC + M )
1.1 bertrand 503: MINWRK = BDSPAC + M
504: ELSE IF( WNTQO ) THEN
505: *
1.16 ! bertrand 506: * Path 2t (N >> M, JOBZ='O')
1.1 bertrand 507: *
1.16 ! bertrand 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 )
1.1 bertrand 514: MAXWRK = WRKBL + 2*M*M
515: MINWRK = BDSPAC + 2*M*M + 3*M
516: ELSE IF( WNTQS ) THEN
517: *
1.16 ! bertrand 518: * Path 3t (N >> M, JOBZ='S')
1.1 bertrand 519: *
1.16 ! bertrand 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 )
1.1 bertrand 526: MAXWRK = WRKBL + M*M
527: MINWRK = BDSPAC + M*M + 3*M
528: ELSE IF( WNTQA ) THEN
529: *
1.16 ! bertrand 530: * Path 4t (N >> M, JOBZ='A')
1.1 bertrand 531: *
1.16 ! bertrand 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 )
1.1 bertrand 538: MAXWRK = WRKBL + M*M
1.16 ! bertrand 539: MINWRK = M*M + MAX( 3*M + BDSPAC, M + N )
1.1 bertrand 540: END IF
541: ELSE
542: *
1.16 ! bertrand 543: * Path 5t (N > M, but not much larger)
1.1 bertrand 544: *
1.16 ! bertrand 545: WRKBL = 3*M + LWORK_DGEBRD_MN
1.1 bertrand 546: IF( WNTQN ) THEN
1.16 ! bertrand 547: * Path 5tn (N > M, jobz='N')
! 548: MAXWRK = MAX( WRKBL, 3*M + BDSPAC )
1.1 bertrand 549: MINWRK = 3*M + MAX( N, BDSPAC )
550: ELSE IF( WNTQO ) THEN
1.16 ! bertrand 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 )
1.1 bertrand 555: MAXWRK = WRKBL + M*N
1.16 ! bertrand 556: MINWRK = 3*M + MAX( N, M*M + BDSPAC )
1.1 bertrand 557: ELSE IF( WNTQS ) THEN
1.16 ! bertrand 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 )
1.1 bertrand 562: MINWRK = 3*M + MAX( N, BDSPAC )
563: ELSE IF( WNTQA ) THEN
1.16 ! bertrand 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 )
1.1 bertrand 568: MINWRK = 3*M + MAX( N, BDSPAC )
569: END IF
570: END IF
571: END IF
1.16 ! bertrand 572:
1.1 bertrand 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: *
1.16 ! bertrand 622: * Path 1 (M >> N, JOBZ='N')
1.1 bertrand 623: * No singular vectors to be computed
624: *
625: ITAU = 1
626: NWORK = ITAU + N
627: *
628: * Compute A=Q*R
1.16 ! bertrand 629: * Workspace: need N [tau] + N [work]
! 630: * Workspace: prefer N [tau] + N*NB [work]
1.1 bertrand 631: *
632: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1.16 ! bertrand 633: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 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
1.16 ! bertrand 644: * Workspace: need 3*N [e, tauq, taup] + N [work]
! 645: * Workspace: prefer 3*N [e, tauq, taup] + 2*N*NB [work]
1.1 bertrand 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
1.16 ! bertrand 653: * Workspace: need N [e] + BDSPAC
1.1 bertrand 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: *
1.16 ! bertrand 660: * Path 2 (M >> N, JOBZ = 'O')
1.1 bertrand 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: *
1.16 ! bertrand 668: IF( LWORK .GE. LDA*N + N*N + 3*N + BDSPAC ) THEN
1.1 bertrand 669: LDWRKR = LDA
670: ELSE
1.16 ! bertrand 671: LDWRKR = ( LWORK - N*N - 3*N - BDSPAC ) / N
1.1 bertrand 672: END IF
673: ITAU = IR + LDWRKR*N
674: NWORK = ITAU + N
675: *
676: * Compute A=Q*R
1.16 ! bertrand 677: * Workspace: need N*N [R] + N [tau] + N [work]
! 678: * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
1.1 bertrand 679: *
680: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1.16 ! bertrand 681: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 682: *
683: * Copy R to WORK(IR), zeroing out below it
684: *
685: CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
1.16 ! bertrand 686: CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
1.1 bertrand 687: $ LDWRKR )
688: *
689: * Generate Q in A
1.16 ! bertrand 690: * Workspace: need N*N [R] + N [tau] + N [work]
! 691: * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
1.1 bertrand 692: *
693: CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
1.16 ! bertrand 694: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 695: IE = ITAU
696: ITAUQ = IE + N
697: ITAUP = ITAUQ + N
698: NWORK = ITAUP + N
699: *
1.16 ! bertrand 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]
1.1 bertrand 703: *
704: CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
705: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1.16 ! bertrand 706: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 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
1.16 ! bertrand 716: * Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 726: *
727: CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
728: $ WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
1.16 ! bertrand 729: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 730: CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
731: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1.16 ! bertrand 732: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 738: *
739: DO 10 I = 1, M, LDWRKR
1.16 ! bertrand 740: CHUNK = MIN( M - I + 1, LDWRKR )
1.1 bertrand 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: *
1.16 ! bertrand 750: * Path 3 (M >> N, JOBZ='S')
1.1 bertrand 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
1.16 ! bertrand 763: * Workspace: need N*N [R] + N [tau] + N [work]
! 764: * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
1.1 bertrand 765: *
766: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1.16 ! bertrand 767: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 768: *
769: * Copy R to WORK(IR), zeroing out below it
770: *
771: CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
1.16 ! bertrand 772: CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1),
1.1 bertrand 773: $ LDWRKR )
774: *
775: * Generate Q in A
1.16 ! bertrand 776: * Workspace: need N*N [R] + N [tau] + N [work]
! 777: * Workspace: prefer N*N [R] + N [tau] + N*NB [work]
1.1 bertrand 778: *
779: CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
1.16 ! bertrand 780: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 781: IE = ITAU
782: ITAUQ = IE + N
783: ITAUP = ITAUQ + N
784: NWORK = ITAUP + N
785: *
786: * Bidiagonalize R in WORK(IR)
1.16 ! bertrand 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]
1.1 bertrand 789: *
790: CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
791: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1.16 ! bertrand 792: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 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
1.16 ! bertrand 797: * Workspace: need N*N [R] + 3*N [e, tauq, taup] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 807: *
808: CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
809: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1.16 ! bertrand 810: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 811: *
812: CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
813: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1.16 ! bertrand 814: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 815: *
816: * Multiply Q in A by left singular vectors of R in
817: * WORK(IR), storing result in U
1.16 ! bertrand 818: * Workspace: need N*N [R]
1.1 bertrand 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: *
1.16 ! bertrand 826: * Path 4 (M >> N, JOBZ='A')
1.1 bertrand 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
1.16 ! bertrand 839: * Workspace: need N*N [U] + N [tau] + N [work]
! 840: * Workspace: prefer N*N [U] + N [tau] + N*NB [work]
1.1 bertrand 841: *
842: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1.16 ! bertrand 843: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 844: CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
845: *
846: * Generate Q in U
1.16 ! bertrand 847: * Workspace: need N*N [U] + N [tau] + M [work]
! 848: * Workspace: prefer N*N [U] + N [tau] + M*NB [work]
1.1 bertrand 849: CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
1.16 ! bertrand 850: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 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
1.16 ! bertrand 871: * Workspace: need N*N [U] + 3*N [e, tauq, taup] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 881: *
882: CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
883: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
1.16 ! bertrand 884: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 885: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
886: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1.16 ! bertrand 887: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 888: *
889: * Multiply Q in U by left singular vectors of R in
890: * WORK(IU), storing result in A
1.16 ! bertrand 891: * Workspace: need N*N [U]
1.1 bertrand 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: *
1.16 ! bertrand 906: * Path 5 (M >= N, but not much larger)
1.1 bertrand 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
1.16 ! bertrand 915: * Workspace: need 3*N [e, tauq, taup] + M [work]
! 916: * Workspace: prefer 3*N [e, tauq, taup] + (M+N)*NB [work]
1.1 bertrand 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: *
1.16 ! bertrand 923: * Path 5n (M >= N, JOBZ='N')
1.1 bertrand 924: * Perform bidiagonal SVD, only computing singular values
1.16 ! bertrand 925: * Workspace: need 3*N [e, tauq, taup] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 930: * Path 5o (M >= N, JOBZ='O')
1.1 bertrand 931: IU = NWORK
1.16 ! bertrand 932: IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN
1.1 bertrand 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 )
1.16 ! bertrand 940: * IR is unused; silence compile warnings
! 941: IR = -1
1.1 bertrand 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
1.16 ! bertrand 952: LDWRKR = ( LWORK - N*N - 3*N ) / N
1.1 bertrand 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
1.16 ! bertrand 959: * Workspace: need 3*N [e, tauq, taup] + N*N [U] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 968: *
969: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
970: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1.16 ! bertrand 971: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 972: *
1.16 ! bertrand 973: IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN
1.1 bertrand 974: *
1.16 ! bertrand 975: * Path 5o-fast
1.1 bertrand 976: * Overwrite WORK(IU) by left singular vectors of A
1.16 ! bertrand 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]
1.1 bertrand 979: *
980: CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
981: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
1.16 ! bertrand 982: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 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: *
1.16 ! bertrand 989: * Path 5o-slow
1.1 bertrand 990: * Generate Q in A
1.16 ! bertrand 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]
1.1 bertrand 993: *
994: CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
1.16 ! bertrand 995: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 1002: *
1003: DO 20 I = 1, M, LDWRKR
1.16 ! bertrand 1004: CHUNK = MIN( M - I + 1, LDWRKR )
1.1 bertrand 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: *
1.16 ! bertrand 1015: * Path 5s (M >= N, JOBZ='S')
1.1 bertrand 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
1.16 ! bertrand 1019: * Workspace: need 3*N [e, tauq, taup] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 1028: * Workspace: need 3*N [e, tauq, taup] + N [work]
! 1029: * Workspace: prefer 3*N [e, tauq, taup] + N*NB [work]
1.1 bertrand 1030: *
1031: CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1032: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1.16 ! bertrand 1033: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1034: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
1035: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1.16 ! bertrand 1036: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1037: ELSE IF( WNTQA ) THEN
1038: *
1.16 ! bertrand 1039: * Path 5a (M >= N, JOBZ='A')
1.1 bertrand 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
1.16 ! bertrand 1043: * Workspace: need 3*N [e, tauq, taup] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 1053: CALL DLASET( 'F', M - N, M - N, ZERO, ONE, U(N+1,N+1),
1.1 bertrand 1054: $ LDU )
1055: END IF
1056: *
1057: * Overwrite U by left singular vectors of A and VT
1058: * by right singular vectors of A
1.16 ! bertrand 1059: * Workspace: need 3*N [e, tauq, taup] + M [work]
! 1060: * Workspace: prefer 3*N [e, tauq, taup] + M*NB [work]
1.1 bertrand 1061: *
1062: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1063: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1.16 ! bertrand 1064: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1065: CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1066: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1.16 ! bertrand 1067: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 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: *
1.16 ! bertrand 1082: * Path 1t (N >> M, JOBZ='N')
1.1 bertrand 1083: * No singular vectors to be computed
1084: *
1085: ITAU = 1
1086: NWORK = ITAU + M
1087: *
1088: * Compute A=L*Q
1.16 ! bertrand 1089: * Workspace: need M [tau] + M [work]
! 1090: * Workspace: prefer M [tau] + M*NB [work]
1.1 bertrand 1091: *
1092: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1.16 ! bertrand 1093: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 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
1.16 ! bertrand 1104: * Workspace: need 3*M [e, tauq, taup] + M [work]
! 1105: * Workspace: prefer 3*M [e, tauq, taup] + 2*M*NB [work]
1.1 bertrand 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
1.16 ! bertrand 1113: * Workspace: need M [e] + BDSPAC
1.1 bertrand 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: *
1.16 ! bertrand 1120: * Path 2t (N >> M, JOBZ='O')
1.1 bertrand 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: *
1.16 ! bertrand 1126: * WORK(IVT) is M by M
! 1127: * WORK(IL) is M by M; it is later resized to M by chunk for gemm
1.1 bertrand 1128: *
1129: IL = IVT + M*M
1.16 ! bertrand 1130: IF( LWORK .GE. M*N + M*M + 3*M + BDSPAC ) THEN
1.1 bertrand 1131: LDWRKL = M
1132: CHUNK = N
1133: ELSE
1134: LDWRKL = M
1.16 ! bertrand 1135: CHUNK = ( LWORK - M*M ) / M
1.1 bertrand 1136: END IF
1137: ITAU = IL + LDWRKL*M
1138: NWORK = ITAU + M
1139: *
1140: * Compute A=L*Q
1.16 ! bertrand 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]
1.1 bertrand 1143: *
1144: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1.16 ! bertrand 1145: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1146: *
1147: * Copy L to WORK(IL), zeroing about above it
1148: *
1149: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1.16 ! bertrand 1150: CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
! 1151: $ WORK( IL + LDWRKL ), LDWRKL )
1.1 bertrand 1152: *
1153: * Generate Q in A
1.16 ! bertrand 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]
1.1 bertrand 1156: *
1157: CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1.16 ! bertrand 1158: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 1159: IE = ITAU
1160: ITAUQ = IE + M
1161: ITAUP = ITAUQ + M
1162: NWORK = ITAUP + M
1163: *
1164: * Bidiagonalize L in WORK(IL)
1.16 ! bertrand 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]
1.1 bertrand 1167: *
1168: CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1169: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1.16 ! bertrand 1170: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 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)
1.16 ! bertrand 1175: * Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 1185: *
1186: CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1187: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1.16 ! bertrand 1188: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1189: CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1190: $ WORK( ITAUP ), WORK( IVT ), M,
1.16 ! bertrand 1191: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 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
1.16 ! bertrand 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.
1.1 bertrand 1198: *
1199: DO 30 I = 1, N, CHUNK
1.16 ! bertrand 1200: BLK = MIN( N - I + 1, CHUNK )
1.1 bertrand 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: *
1.16 ! bertrand 1209: * Path 3t (N >> M, JOBZ='S')
1.1 bertrand 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
1.16 ! bertrand 1222: * Workspace: need M*M [L] + M [tau] + M [work]
! 1223: * Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1.1 bertrand 1224: *
1225: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1.16 ! bertrand 1226: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1227: *
1228: * Copy L to WORK(IL), zeroing out above it
1229: *
1230: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1.16 ! bertrand 1231: CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO,
! 1232: $ WORK( IL + LDWRKL ), LDWRKL )
1.1 bertrand 1233: *
1234: * Generate Q in A
1.16 ! bertrand 1235: * Workspace: need M*M [L] + M [tau] + M [work]
! 1236: * Workspace: prefer M*M [L] + M [tau] + M*NB [work]
1.1 bertrand 1237: *
1238: CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1.16 ! bertrand 1239: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 1240: IE = ITAU
1241: ITAUQ = IE + M
1242: ITAUP = ITAUQ + M
1243: NWORK = ITAUP + M
1244: *
1.16 ! bertrand 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]
1.1 bertrand 1248: *
1249: CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1250: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1.16 ! bertrand 1251: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 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
1.16 ! bertrand 1256: * Workspace: need M*M [L] + 3*M [e, tauq, taup] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 1266: *
1267: CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1268: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1.16 ! bertrand 1269: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1270: CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1271: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1.16 ! bertrand 1272: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1273: *
1274: * Multiply right singular vectors of L in WORK(IL) by
1275: * Q in A, storing result in VT
1.16 ! bertrand 1276: * Workspace: need M*M [L]
1.1 bertrand 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: *
1.16 ! bertrand 1284: * Path 4t (N >> M, JOBZ='A')
1.1 bertrand 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
1.16 ! bertrand 1297: * Workspace: need M*M [VT] + M [tau] + M [work]
! 1298: * Workspace: prefer M*M [VT] + M [tau] + M*NB [work]
1.1 bertrand 1299: *
1300: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1.16 ! bertrand 1301: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1302: CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
1303: *
1304: * Generate Q in VT
1.16 ! bertrand 1305: * Workspace: need M*M [VT] + M [tau] + N [work]
! 1306: * Workspace: prefer M*M [VT] + M [tau] + N*NB [work]
1.1 bertrand 1307: *
1308: CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1.16 ! bertrand 1309: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 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)
1.16 ! bertrand 1330: * Workspace: need M*M [VT] + 3*M [e, tauq, taup] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 1340: *
1341: CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1342: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1.16 ! bertrand 1343: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1344: CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
1345: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1.16 ! bertrand 1346: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 1347: *
1348: * Multiply right singular vectors of L in WORK(IVT) by
1349: * Q in VT, storing result in A
1.16 ! bertrand 1350: * Workspace: need M*M [VT]
1.1 bertrand 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: *
1.16 ! bertrand 1365: * Path 5t (N > M, but not much larger)
1.1 bertrand 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
1.16 ! bertrand 1374: * Workspace: need 3*M [e, tauq, taup] + N [work]
! 1375: * Workspace: prefer 3*M [e, tauq, taup] + (M+N)*NB [work]
1.1 bertrand 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: *
1.16 ! bertrand 1382: * Path 5tn (N > M, JOBZ='N')
1.1 bertrand 1383: * Perform bidiagonal SVD, only computing singular values
1.16 ! bertrand 1384: * Workspace: need 3*M [e, tauq, taup] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 1389: * Path 5to (N > M, JOBZ='O')
1.1 bertrand 1390: LDWKVT = M
1391: IVT = NWORK
1.16 ! bertrand 1392: IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN
1.1 bertrand 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
1.16 ! bertrand 1399: * IL is unused; silence compile warnings
! 1400: IL = -1
1.1 bertrand 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: *
1.16 ! bertrand 1410: CHUNK = ( LWORK - M*M - 3*M ) / M
1.1 bertrand 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)
1.16 ! bertrand 1416: * Workspace: need 3*M [e, tauq, taup] + M*M [VT] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 1425: *
1426: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1427: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1.16 ! bertrand 1428: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1429: *
1.16 ! bertrand 1430: IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN
1.1 bertrand 1431: *
1.16 ! bertrand 1432: * Path 5to-fast
1.1 bertrand 1433: * Overwrite WORK(IVT) by left singular vectors of A
1.16 ! bertrand 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]
1.1 bertrand 1436: *
1437: CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1438: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1.16 ! bertrand 1439: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 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: *
1.16 ! bertrand 1446: * Path 5to-slow
1.1 bertrand 1447: * Generate P**T in A
1.16 ! bertrand 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]
1.1 bertrand 1450: *
1451: CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1.16 ! bertrand 1452: $ WORK( NWORK ), LWORK - NWORK + 1, IERR )
1.1 bertrand 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
1.16 ! bertrand 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]
1.1 bertrand 1459: *
1460: DO 40 I = 1, N, CHUNK
1.16 ! bertrand 1461: BLK = MIN( N - I + 1, CHUNK )
1.1 bertrand 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: *
1.16 ! bertrand 1471: * Path 5ts (N > M, JOBZ='S')
1.1 bertrand 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
1.16 ! bertrand 1475: * Workspace: need 3*M [e, tauq, taup] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 1484: * Workspace: need 3*M [e, tauq, taup] + M [work]
! 1485: * Workspace: prefer 3*M [e, tauq, taup] + M*NB [work]
1.1 bertrand 1486: *
1487: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1488: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1.16 ! bertrand 1489: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1490: CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1491: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1.16 ! bertrand 1492: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1493: ELSE IF( WNTQA ) THEN
1494: *
1.16 ! bertrand 1495: * Path 5ta (N > M, JOBZ='A')
1.1 bertrand 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
1.16 ! bertrand 1499: * Workspace: need 3*M [e, tauq, taup] + BDSPAC
1.1 bertrand 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
1.16 ! bertrand 1509: CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT(M+1,M+1),
1.1 bertrand 1510: $ LDVT )
1511: END IF
1512: *
1513: * Overwrite U by left singular vectors of A and VT
1514: * by right singular vectors of A
1.16 ! bertrand 1515: * Workspace: need 3*M [e, tauq, taup] + N [work]
! 1516: * Workspace: prefer 3*M [e, tauq, taup] + N*NB [work]
1.1 bertrand 1517: *
1518: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1519: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1.16 ! bertrand 1520: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 1521: CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1522: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1.16 ! bertrand 1523: $ LWORK - NWORK + 1, IERR )
1.1 bertrand 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