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