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