Annotation of rpl/lapack/lapack/dgesvdx.f, revision 1.8
1.1 bertrand 1: *> \brief <b> DGESVDX computes the singular value decomposition (SVD) for GE matrices</b>
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.4 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.1 bertrand 7: *
8: *> \htmlonly
1.4 bertrand 9: *> Download DGESVDX + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvdx.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvdx.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvdx.f">
1.1 bertrand 15: *> [TXT]</a>
1.4 bertrand 16: *> \endhtmlonly
1.1 bertrand 17: *
18: * Definition:
19: * ===========
20: *
1.4 bertrand 21: * SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
22: * $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
1.1 bertrand 23: * $ LWORK, IWORK, INFO )
1.4 bertrand 24: *
1.1 bertrand 25: *
26: * .. Scalar Arguments ..
27: * CHARACTER JOBU, JOBVT, RANGE
28: * INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
29: * DOUBLE PRECISION VL, VU
30: * ..
31: * .. Array Arguments ..
32: * INTEGER IWORK( * )
33: * DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
34: * $ VT( LDVT, * ), WORK( * )
35: * ..
1.4 bertrand 36: *
1.1 bertrand 37: *
38: *> \par Purpose:
39: * =============
40: *>
41: *> \verbatim
42: *>
43: *> DGESVDX computes the singular value decomposition (SVD) of a real
44: *> M-by-N matrix A, optionally computing the left and/or right singular
45: *> vectors. The SVD is written
1.4 bertrand 46: *>
1.1 bertrand 47: *> A = U * SIGMA * transpose(V)
1.4 bertrand 48: *>
1.1 bertrand 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.
1.4 bertrand 55: *>
56: *> DGESVDX uses an eigenvalue problem for obtaining the SVD, which
57: *> allows for the computation of a subset of singular values and
1.1 bertrand 58: *> vectors. See DBDSVDX for details.
1.4 bertrand 59: *>
1.1 bertrand 60: *> Note that the routine returns V**T, not V.
61: *> \endverbatim
1.4 bertrand 62: *
1.1 bertrand 63: * Arguments:
64: * ==========
65: *
66: *> \param[in] JOBU
67: *> \verbatim
68: *> JOBU is CHARACTER*1
69: *> Specifies options for computing all or part of the matrix U:
70: *> = 'V': the first min(m,n) columns of U (the left singular
1.4 bertrand 71: *> vectors) or as specified by RANGE are returned in
1.1 bertrand 72: *> the array U;
73: *> = 'N': no columns of U (no left singular vectors) are
74: *> computed.
75: *> \endverbatim
76: *>
77: *> \param[in] JOBVT
78: *> \verbatim
79: *> JOBVT is CHARACTER*1
80: *> Specifies options for computing all or part of the matrix
81: *> V**T:
82: *> = 'V': the first min(m,n) rows of V**T (the right singular
1.4 bertrand 83: *> vectors) or as specified by RANGE are returned in
1.1 bertrand 84: *> the array VT;
85: *> = 'N': no rows of V**T (no right singular vectors) are
86: *> computed.
87: *> \endverbatim
88: *>
89: *> \param[in] RANGE
90: *> \verbatim
91: *> RANGE is CHARACTER*1
92: *> = 'A': all singular values will be found.
93: *> = 'V': all singular values in the half-open interval (VL,VU]
94: *> will be found.
1.4 bertrand 95: *> = 'I': the IL-th through IU-th singular values will be found.
1.1 bertrand 96: *> \endverbatim
97: *>
98: *> \param[in] M
99: *> \verbatim
100: *> M is INTEGER
101: *> The number of rows of the input matrix A. M >= 0.
102: *> \endverbatim
103: *>
104: *> \param[in] N
105: *> \verbatim
106: *> N is INTEGER
107: *> The number of columns of the input matrix A. N >= 0.
108: *> \endverbatim
109: *>
110: *> \param[in,out] A
111: *> \verbatim
112: *> A is DOUBLE PRECISION array, dimension (LDA,N)
113: *> On entry, the M-by-N matrix A.
114: *> On exit, the contents of A are destroyed.
115: *> \endverbatim
116: *>
117: *> \param[in] LDA
118: *> \verbatim
119: *> LDA is INTEGER
120: *> The leading dimension of the array A. LDA >= max(1,M).
121: *> \endverbatim
122: *>
123: *> \param[in] VL
124: *> \verbatim
125: *> VL is DOUBLE PRECISION
1.2 bertrand 126: *> If RANGE='V', the lower bound of the interval to
127: *> be searched for singular values. VU > VL.
128: *> Not referenced if RANGE = 'A' or 'I'.
1.1 bertrand 129: *> \endverbatim
130: *>
131: *> \param[in] VU
132: *> \verbatim
133: *> VU is DOUBLE PRECISION
1.2 bertrand 134: *> If RANGE='V', the upper bound of the interval to
1.1 bertrand 135: *> be searched for singular values. VU > VL.
136: *> Not referenced if RANGE = 'A' or 'I'.
137: *> \endverbatim
138: *>
139: *> \param[in] IL
140: *> \verbatim
141: *> IL is INTEGER
1.2 bertrand 142: *> If RANGE='I', the index of the
143: *> smallest singular value to be returned.
144: *> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
145: *> Not referenced if RANGE = 'A' or 'V'.
1.1 bertrand 146: *> \endverbatim
147: *>
148: *> \param[in] IU
149: *> \verbatim
150: *> IU is INTEGER
1.2 bertrand 151: *> If RANGE='I', the index of the
152: *> largest singular value to be returned.
1.1 bertrand 153: *> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
154: *> Not referenced if RANGE = 'A' or 'V'.
155: *> \endverbatim
156: *>
157: *> \param[out] NS
158: *> \verbatim
159: *> NS is INTEGER
1.4 bertrand 160: *> The total number of singular values found,
1.1 bertrand 161: *> 0 <= NS <= min(M,N).
162: *> If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1.
163: *> \endverbatim
164: *>
165: *> \param[out] S
166: *> \verbatim
167: *> S is DOUBLE PRECISION array, dimension (min(M,N))
168: *> The singular values of A, sorted so that S(i) >= S(i+1).
169: *> \endverbatim
170: *>
171: *> \param[out] U
172: *> \verbatim
173: *> U is DOUBLE PRECISION array, dimension (LDU,UCOL)
1.4 bertrand 174: *> If JOBU = 'V', U contains columns of U (the left singular
175: *> vectors, stored columnwise) as specified by RANGE; if
1.1 bertrand 176: *> JOBU = 'N', U is not referenced.
1.4 bertrand 177: *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V',
1.2 bertrand 178: *> the exact value of NS is not known in advance and an upper
1.1 bertrand 179: *> bound must be used.
180: *> \endverbatim
181: *>
182: *> \param[in] LDU
183: *> \verbatim
184: *> LDU is INTEGER
185: *> The leading dimension of the array U. LDU >= 1; if
186: *> JOBU = 'V', LDU >= M.
187: *> \endverbatim
188: *>
189: *> \param[out] VT
190: *> \verbatim
191: *> VT is DOUBLE PRECISION array, dimension (LDVT,N)
1.4 bertrand 192: *> If JOBVT = 'V', VT contains the rows of V**T (the right singular
193: *> vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N',
1.1 bertrand 194: *> VT is not referenced.
1.4 bertrand 195: *> Note: The user must ensure that LDVT >= NS; if RANGE = 'V',
196: *> the exact value of NS is not known in advance and an upper
1.1 bertrand 197: *> bound must be used.
198: *> \endverbatim
199: *>
200: *> \param[in] LDVT
201: *> \verbatim
202: *> LDVT is INTEGER
203: *> The leading dimension of the array VT. LDVT >= 1; if
204: *> JOBVT = 'V', LDVT >= NS (see above).
205: *> \endverbatim
206: *>
207: *> \param[out] WORK
208: *> \verbatim
209: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
210: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
211: *> \endverbatim
212: *>
213: *> \param[in] LWORK
214: *> \verbatim
215: *> LWORK is INTEGER
216: *> The dimension of the array WORK.
1.4 bertrand 217: *> LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see
1.1 bertrand 218: *> comments inside the code):
1.4 bertrand 219: *> - PATH 1 (M much larger than N)
1.1 bertrand 220: *> - PATH 1t (N much larger than M)
221: *> LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths.
222: *> For good performance, LWORK should generally be larger.
223: *>
224: *> If LWORK = -1, then a workspace query is assumed; the routine
225: *> only calculates the optimal size of the WORK array, returns
226: *> this value as the first entry of the WORK array, and no error
227: *> message related to LWORK is issued by XERBLA.
228: *> \endverbatim
229: *>
230: *> \param[out] IWORK
231: *> \verbatim
232: *> IWORK is INTEGER array, dimension (12*MIN(M,N))
1.4 bertrand 233: *> If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0,
234: *> then IWORK contains the indices of the eigenvectors that failed
1.1 bertrand 235: *> to converge in DBDSVDX/DSTEVX.
236: *> \endverbatim
237: *>
238: *> \param[out] INFO
239: *> \verbatim
240: *> INFO is INTEGER
241: *> = 0: successful exit
242: *> < 0: if INFO = -i, the i-th argument had an illegal value
243: *> > 0: if INFO = i, then i eigenvectors failed to converge
244: *> in DBDSVDX/DSTEVX.
245: *> if INFO = N*2 + 1, an internal error occurred in
246: *> DBDSVDX
247: *> \endverbatim
248: *
249: * Authors:
250: * ========
251: *
1.4 bertrand 252: *> \author Univ. of Tennessee
253: *> \author Univ. of California Berkeley
254: *> \author Univ. of Colorado Denver
255: *> \author NAG Ltd.
1.1 bertrand 256: *
257: *> \ingroup doubleGEsing
258: *
259: * =====================================================================
1.4 bertrand 260: SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
261: $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
1.1 bertrand 262: $ LWORK, IWORK, INFO )
263: *
1.8 ! bertrand 264: * -- LAPACK driver routine --
1.1 bertrand 265: * -- LAPACK is a software package provided by Univ. of Tennessee, --
266: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
267: *
268: * .. Scalar Arguments ..
269: CHARACTER JOBU, JOBVT, RANGE
270: INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
271: DOUBLE PRECISION VL, VU
272: * ..
273: * .. Array Arguments ..
274: INTEGER IWORK( * )
275: DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
276: $ VT( LDVT, * ), WORK( * )
277: * ..
278: *
279: * =====================================================================
280: *
281: * .. Parameters ..
282: DOUBLE PRECISION ZERO, ONE
283: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
284: * ..
285: * .. Local Scalars ..
286: CHARACTER JOBZ, RNGTGK
287: LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
288: INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
1.4 bertrand 289: $ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK,
1.1 bertrand 290: $ J, MAXWRK, MINMN, MINWRK, MNTHR
291: DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
292: * ..
293: * .. Local Arrays ..
294: DOUBLE PRECISION DUM( 1 )
295: * ..
296: * .. External Subroutines ..
297: EXTERNAL DBDSVDX, DGEBRD, DGELQF, DGEQRF, DLACPY,
298: $ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR,
1.6 bertrand 299: $ DCOPY, XERBLA
1.1 bertrand 300: * ..
301: * .. External Functions ..
302: LOGICAL LSAME
303: INTEGER ILAENV
1.4 bertrand 304: DOUBLE PRECISION DLAMCH, DLANGE
305: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
1.1 bertrand 306: * ..
307: * .. Intrinsic Functions ..
308: INTRINSIC MAX, MIN, SQRT
309: * ..
310: * .. Executable Statements ..
311: *
312: * Test the input arguments.
313: *
314: NS = 0
315: INFO = 0
316: ABSTOL = 2*DLAMCH('S')
317: LQUERY = ( LWORK.EQ.-1 )
318: MINMN = MIN( M, N )
319:
320: WANTU = LSAME( JOBU, 'V' )
321: WANTVT = LSAME( JOBVT, 'V' )
322: IF( WANTU .OR. WANTVT ) THEN
323: JOBZ = 'V'
324: ELSE
325: JOBZ = 'N'
326: END IF
327: ALLS = LSAME( RANGE, 'A' )
328: VALS = LSAME( RANGE, 'V' )
329: INDS = LSAME( RANGE, 'I' )
330: *
331: INFO = 0
332: IF( .NOT.LSAME( JOBU, 'V' ) .AND.
333: $ .NOT.LSAME( JOBU, 'N' ) ) THEN
334: INFO = -1
335: ELSE IF( .NOT.LSAME( JOBVT, 'V' ) .AND.
336: $ .NOT.LSAME( JOBVT, 'N' ) ) THEN
337: INFO = -2
338: ELSE IF( .NOT.( ALLS .OR. VALS .OR. INDS ) ) THEN
339: INFO = -3
340: ELSE IF( M.LT.0 ) THEN
341: INFO = -4
342: ELSE IF( N.LT.0 ) THEN
343: INFO = -5
344: ELSE IF( M.GT.LDA ) THEN
345: INFO = -7
346: ELSE IF( MINMN.GT.0 ) THEN
347: IF( VALS ) THEN
348: IF( VL.LT.ZERO ) THEN
349: INFO = -8
350: ELSE IF( VU.LE.VL ) THEN
351: INFO = -9
352: END IF
353: ELSE IF( INDS ) THEN
354: IF( IL.LT.1 .OR. IL.GT.MAX( 1, MINMN ) ) THEN
355: INFO = -10
356: ELSE IF( IU.LT.MIN( MINMN, IL ) .OR. IU.GT.MINMN ) THEN
357: INFO = -11
358: END IF
359: END IF
360: IF( INFO.EQ.0 ) THEN
361: IF( WANTU .AND. LDU.LT.M ) THEN
362: INFO = -15
1.2 bertrand 363: ELSE IF( WANTVT ) THEN
364: IF( INDS ) THEN
365: IF( LDVT.LT.IU-IL+1 ) THEN
366: INFO = -17
367: END IF
368: ELSE IF( LDVT.LT.MINMN ) THEN
369: INFO = -17
370: END IF
1.1 bertrand 371: END IF
372: END IF
373: END IF
374: *
375: * Compute workspace
376: * (Note: Comments in the code beginning "Workspace:" describe the
377: * minimal amount of workspace needed at that point in the code,
378: * as well as the preferred amount for good performance.
379: * NB refers to the optimal block size for the immediately
380: * following subroutine, as returned by ILAENV.)
381: *
382: IF( INFO.EQ.0 ) THEN
383: MINWRK = 1
384: MAXWRK = 1
385: IF( MINMN.GT.0 ) THEN
386: IF( M.GE.N ) THEN
387: MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
388: IF( M.GE.MNTHR ) THEN
389: *
390: * Path 1 (M much larger than N)
391: *
1.4 bertrand 392: MAXWRK = N +
1.1 bertrand 393: $ N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
1.2 bertrand 394: MAXWRK = MAX( MAXWRK, N*(N+5) + 2*N*
1.1 bertrand 395: $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
1.2 bertrand 396: IF (WANTU) THEN
397: MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
398: $ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) )
399: END IF
400: IF (WANTVT) THEN
401: MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
402: $ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) )
403: END IF
404: MINWRK = N*(N*3+20)
1.1 bertrand 405: ELSE
406: *
407: * Path 2 (M at least N, but not much larger)
408: *
1.2 bertrand 409: MAXWRK = 4*N + ( M+N )*
1.1 bertrand 410: $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
1.2 bertrand 411: IF (WANTU) THEN
412: MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
413: $ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) )
414: END IF
415: IF (WANTVT) THEN
416: MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
417: $ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) )
418: END IF
419: MINWRK = MAX(N*(N*2+19),4*N+M)
1.1 bertrand 420: END IF
421: ELSE
422: MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
423: IF( N.GE.MNTHR ) THEN
424: *
425: * Path 1t (N much larger than M)
426: *
1.4 bertrand 427: MAXWRK = M +
1.1 bertrand 428: $ M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
1.2 bertrand 429: MAXWRK = MAX( MAXWRK, M*(M+5) + 2*M*
1.1 bertrand 430: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
1.2 bertrand 431: IF (WANTU) THEN
432: MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
433: $ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) )
434: END IF
435: IF (WANTVT) THEN
436: MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
437: $ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) )
438: END IF
439: MINWRK = M*(M*3+20)
1.1 bertrand 440: ELSE
441: *
1.2 bertrand 442: * Path 2t (N at least M, but not much larger)
1.1 bertrand 443: *
1.2 bertrand 444: MAXWRK = 4*M + ( M+N )*
1.1 bertrand 445: $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
1.2 bertrand 446: IF (WANTU) THEN
447: MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
448: $ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) )
449: END IF
450: IF (WANTVT) THEN
451: MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
452: $ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) )
453: END IF
454: MINWRK = MAX(M*(M*2+19),4*M+N)
1.1 bertrand 455: END IF
456: END IF
457: END IF
458: MAXWRK = MAX( MAXWRK, MINWRK )
459: WORK( 1 ) = DBLE( MAXWRK )
460: *
461: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
462: INFO = -19
463: END IF
464: END IF
465: *
466: IF( INFO.NE.0 ) THEN
467: CALL XERBLA( 'DGESVDX', -INFO )
468: RETURN
469: ELSE IF( LQUERY ) THEN
470: RETURN
471: END IF
472: *
473: * Quick return if possible
474: *
475: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
476: RETURN
477: END IF
478: *
479: * Set singular values indices accord to RANGE.
480: *
481: IF( ALLS ) THEN
482: RNGTGK = 'I'
483: ILTGK = 1
484: IUTGK = MIN( M, N )
485: ELSE IF( INDS ) THEN
486: RNGTGK = 'I'
487: ILTGK = IL
488: IUTGK = IU
1.4 bertrand 489: ELSE
1.1 bertrand 490: RNGTGK = 'V'
491: ILTGK = 0
492: IUTGK = 0
493: END IF
494: *
495: * Get machine constants
496: *
497: EPS = DLAMCH( 'P' )
498: SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
499: BIGNUM = ONE / SMLNUM
500: *
501: * Scale A if max element outside range [SMLNUM,BIGNUM]
502: *
503: ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
504: ISCL = 0
505: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
506: ISCL = 1
507: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
508: ELSE IF( ANRM.GT.BIGNUM ) THEN
509: ISCL = 1
510: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
511: END IF
512: *
513: IF( M.GE.N ) THEN
514: *
515: * A has at least as many rows as columns. If A has sufficiently
516: * more rows than columns, first reduce A using the QR
517: * decomposition.
518: *
519: IF( M.GE.MNTHR ) THEN
520: *
521: * Path 1 (M much larger than N):
522: * A = Q * R = Q * ( QB * B * PB**T )
523: * = Q * ( QB * ( UB * S * VB**T ) * PB**T )
524: * U = Q * QB * UB; V**T = VB**T * PB**T
525: *
526: * Compute A=Q*R
527: * (Workspace: need 2*N, prefer N+N*NB)
528: *
529: ITAU = 1
530: ITEMP = ITAU + N
531: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),
532: $ LWORK-ITEMP+1, INFO )
1.4 bertrand 533: *
1.1 bertrand 534: * Copy R into WORK and bidiagonalize it:
535: * (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB)
536: *
537: IQRF = ITEMP
538: ID = IQRF + N*N
539: IE = ID + N
540: ITAUQ = IE + N
541: ITAUP = ITAUQ + N
1.4 bertrand 542: ITEMP = ITAUP + N
1.1 bertrand 543: CALL DLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N )
544: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IQRF+1 ), N )
1.4 bertrand 545: CALL DGEBRD( N, N, WORK( IQRF ), N, WORK( ID ), WORK( IE ),
1.1 bertrand 546: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
547: $ LWORK-ITEMP+1, INFO )
548: *
549: * Solve eigenvalue problem TGK*Z=Z*S.
1.4 bertrand 550: * (Workspace: need 14*N + 2*N*(N+1))
551: *
1.1 bertrand 552: ITGKZ = ITEMP
553: ITEMP = ITGKZ + N*(N*2+1)
1.4 bertrand 554: CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
1.1 bertrand 555: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
556: $ N*2, WORK( ITEMP ), IWORK, INFO)
557: *
558: * If needed, compute left singular vectors.
559: *
560: IF( WANTU ) THEN
561: J = ITGKZ
562: DO I = 1, NS
563: CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
564: J = J + N*2
565: END DO
1.2 bertrand 566: CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU )
1.1 bertrand 567: *
568: * Call DORMBR to compute QB*UB.
569: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
570: *
1.4 bertrand 571: CALL DORMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N,
572: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
1.1 bertrand 573: $ LWORK-ITEMP+1, INFO )
574: *
575: * Call DORMQR to compute Q*(QB*UB).
576: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
577: *
1.4 bertrand 578: CALL DORMQR( 'L', 'N', M, NS, N, A, LDA,
1.1 bertrand 579: $ WORK( ITAU ), U, LDU, WORK( ITEMP ),
580: $ LWORK-ITEMP+1, INFO )
1.4 bertrand 581: END IF
582: *
1.1 bertrand 583: * If needed, compute right singular vectors.
584: *
585: IF( WANTVT) THEN
586: J = ITGKZ + N
587: DO I = 1, NS
588: CALL DCOPY( N, WORK( J ), 1, VT( I,1 ), LDVT )
589: J = J + N*2
590: END DO
591: *
592: * Call DORMBR to compute VB**T * PB**T
593: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
594: *
1.4 bertrand 595: CALL DORMBR( 'P', 'R', 'T', NS, N, N, WORK( IQRF ), N,
1.1 bertrand 596: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
597: $ LWORK-ITEMP+1, INFO )
598: END IF
599: ELSE
600: *
601: * Path 2 (M at least N, but not much larger)
602: * Reduce A to bidiagonal form without QR decomposition
603: * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
604: * U = QB * UB; V**T = VB**T * PB**T
605: *
606: * Bidiagonalize A
607: * (Workspace: need 4*N+M, prefer 4*N+(M+N)*NB)
608: *
609: ID = 1
610: IE = ID + N
611: ITAUQ = IE + N
612: ITAUP = ITAUQ + N
1.4 bertrand 613: ITEMP = ITAUP + N
614: CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
1.1 bertrand 615: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
616: $ LWORK-ITEMP+1, INFO )
617: *
618: * Solve eigenvalue problem TGK*Z=Z*S.
1.4 bertrand 619: * (Workspace: need 14*N + 2*N*(N+1))
620: *
1.1 bertrand 621: ITGKZ = ITEMP
622: ITEMP = ITGKZ + N*(N*2+1)
1.4 bertrand 623: CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
1.1 bertrand 624: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
625: $ N*2, WORK( ITEMP ), IWORK, INFO)
626: *
627: * If needed, compute left singular vectors.
628: *
629: IF( WANTU ) THEN
630: J = ITGKZ
631: DO I = 1, NS
632: CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
633: J = J + N*2
634: END DO
1.2 bertrand 635: CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU )
1.1 bertrand 636: *
637: * Call DORMBR to compute QB*UB.
638: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
1.4 bertrand 639: *
640: CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
641: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
1.1 bertrand 642: $ LWORK-ITEMP+1, IERR )
1.4 bertrand 643: END IF
644: *
1.1 bertrand 645: * If needed, compute right singular vectors.
646: *
647: IF( WANTVT) THEN
648: J = ITGKZ + N
649: DO I = 1, NS
650: CALL DCOPY( N, WORK( J ), 1, VT( I,1 ), LDVT )
651: J = J + N*2
652: END DO
653: *
654: * Call DORMBR to compute VB**T * PB**T
655: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
656: *
1.4 bertrand 657: CALL DORMBR( 'P', 'R', 'T', NS, N, N, A, LDA,
1.1 bertrand 658: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
659: $ LWORK-ITEMP+1, IERR )
660: END IF
1.4 bertrand 661: END IF
1.1 bertrand 662: ELSE
663: *
664: * A has more columns than rows. If A has sufficiently more
665: * columns than rows, first reduce A using the LQ decomposition.
666: *
667: IF( N.GE.MNTHR ) THEN
668: *
669: * Path 1t (N much larger than M):
1.4 bertrand 670: * A = L * Q = ( QB * B * PB**T ) * Q
1.1 bertrand 671: * = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
672: * U = QB * UB ; V**T = VB**T * PB**T * Q
673: *
674: * Compute A=L*Q
675: * (Workspace: need 2*M, prefer M+M*NB)
676: *
677: ITAU = 1
678: ITEMP = ITAU + M
679: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),
680: $ LWORK-ITEMP+1, INFO )
681:
682: * Copy L into WORK and bidiagonalize it:
683: * (Workspace in WORK( ITEMP ): need M*M+5*N, prefer M*M+4*M+2*M*NB)
684: *
685: ILQF = ITEMP
686: ID = ILQF + M*M
687: IE = ID + M
688: ITAUQ = IE + M
689: ITAUP = ITAUQ + M
690: ITEMP = ITAUP + M
691: CALL DLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M )
692: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( ILQF+M ), M )
1.4 bertrand 693: CALL DGEBRD( M, M, WORK( ILQF ), M, WORK( ID ), WORK( IE ),
1.1 bertrand 694: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
695: $ LWORK-ITEMP+1, INFO )
696: *
697: * Solve eigenvalue problem TGK*Z=Z*S.
1.4 bertrand 698: * (Workspace: need 2*M*M+14*M)
1.1 bertrand 699: *
700: ITGKZ = ITEMP
701: ITEMP = ITGKZ + M*(M*2+1)
1.4 bertrand 702: CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
1.1 bertrand 703: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
704: $ M*2, WORK( ITEMP ), IWORK, INFO)
705: *
706: * If needed, compute left singular vectors.
707: *
708: IF( WANTU ) THEN
709: J = ITGKZ
710: DO I = 1, NS
711: CALL DCOPY( M, WORK( J ), 1, U( 1,I ), 1 )
712: J = J + M*2
713: END DO
714: *
715: * Call DORMBR to compute QB*UB.
716: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
717: *
1.4 bertrand 718: CALL DORMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M,
719: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
1.1 bertrand 720: $ LWORK-ITEMP+1, INFO )
1.4 bertrand 721: END IF
722: *
1.1 bertrand 723: * If needed, compute right singular vectors.
724: *
725: IF( WANTVT) THEN
726: J = ITGKZ + M
727: DO I = 1, NS
728: CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
729: J = J + M*2
730: END DO
1.2 bertrand 731: CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT)
1.1 bertrand 732: *
733: * Call DORMBR to compute (VB**T)*(PB**T)
734: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
735: *
1.4 bertrand 736: CALL DORMBR( 'P', 'R', 'T', NS, M, M, WORK( ILQF ), M,
1.1 bertrand 737: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
738: $ LWORK-ITEMP+1, INFO )
739: *
740: * Call DORMLQ to compute ((VB**T)*(PB**T))*Q.
741: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
742: *
1.4 bertrand 743: CALL DORMLQ( 'R', 'N', NS, N, M, A, LDA,
1.1 bertrand 744: $ WORK( ITAU ), VT, LDVT, WORK( ITEMP ),
745: $ LWORK-ITEMP+1, INFO )
1.4 bertrand 746: END IF
1.1 bertrand 747: ELSE
748: *
749: * Path 2t (N greater than M, but not much larger)
750: * Reduce to bidiagonal form without LQ decomposition
751: * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
1.4 bertrand 752: * U = QB * UB; V**T = VB**T * PB**T
1.1 bertrand 753: *
754: * Bidiagonalize A
755: * (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB)
756: *
757: ID = 1
758: IE = ID + M
759: ITAUQ = IE + M
760: ITAUP = ITAUQ + M
761: ITEMP = ITAUP + M
1.4 bertrand 762: CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
1.1 bertrand 763: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
764: $ LWORK-ITEMP+1, INFO )
765: *
766: * Solve eigenvalue problem TGK*Z=Z*S.
1.4 bertrand 767: * (Workspace: need 2*M*M+14*M)
1.1 bertrand 768: *
769: ITGKZ = ITEMP
770: ITEMP = ITGKZ + M*(M*2+1)
1.4 bertrand 771: CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
1.1 bertrand 772: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
773: $ M*2, WORK( ITEMP ), IWORK, INFO)
1.4 bertrand 774: *
1.1 bertrand 775: * If needed, compute left singular vectors.
776: *
777: IF( WANTU ) THEN
778: J = ITGKZ
779: DO I = 1, NS
780: CALL DCOPY( M, WORK( J ), 1, U( 1,I ), 1 )
781: J = J + M*2
782: END DO
783: *
784: * Call DORMBR to compute QB*UB.
785: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
786: *
1.4 bertrand 787: CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
788: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
1.1 bertrand 789: $ LWORK-ITEMP+1, INFO )
1.4 bertrand 790: END IF
791: *
1.1 bertrand 792: * If needed, compute right singular vectors.
793: *
794: IF( WANTVT) THEN
795: J = ITGKZ + M
796: DO I = 1, NS
797: CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
798: J = J + M*2
799: END DO
1.2 bertrand 800: CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT)
1.1 bertrand 801: *
802: * Call DORMBR to compute VB**T * PB**T
803: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
804: *
1.4 bertrand 805: CALL DORMBR( 'P', 'R', 'T', NS, N, M, A, LDA,
1.1 bertrand 806: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
807: $ LWORK-ITEMP+1, INFO )
1.4 bertrand 808: END IF
1.1 bertrand 809: END IF
810: END IF
811: *
812: * Undo scaling if necessary
813: *
814: IF( ISCL.EQ.1 ) THEN
815: IF( ANRM.GT.BIGNUM )
816: $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1,
817: $ S, MINMN, INFO )
818: IF( ANRM.LT.SMLNUM )
819: $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1,
820: $ S, MINMN, INFO )
821: END IF
822: *
823: * Return optimal workspace in WORK(1)
824: *
825: WORK( 1 ) = DBLE( MAXWRK )
826: *
827: RETURN
828: *
829: * End of DGESVDX
830: *
831: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgesvdx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack