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: *> 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'.
129: *> \endverbatim
130: *>
131: *> \param[in] VU
132: *> \verbatim
133: *> VU is DOUBLE PRECISION
134: *> If RANGE='V', the upper bound of the interval to
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
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'.
146: *> \endverbatim
147: *>
148: *> \param[in] IU
149: *> \verbatim
150: *> IU is INTEGER
151: *> If RANGE='I', the index of the
152: *> largest singular value to be returned.
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
160: *> The total number of singular values found,
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)
174: *> If JOBU = 'V', U contains columns of U (the left singular
175: *> vectors, stored columnwise) as specified by RANGE; if
176: *> JOBU = 'N', U is not referenced.
177: *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V',
178: *> the exact value of NS is not known in advance and an upper
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)
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',
194: *> VT is not referenced.
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
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.
217: *> LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see
218: *> comments inside the code):
219: *> - PATH 1 (M much larger than N)
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))
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
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: *
252: *> \author Univ. of Tennessee
253: *> \author Univ. of California Berkeley
254: *> \author Univ. of Colorado Denver
255: *> \author NAG Ltd.
256: *
257: *> \date June 2016
258: *
259: *> \ingroup doubleGEsing
260: *
261: * =====================================================================
262: SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
263: $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
264: $ LWORK, IWORK, INFO )
265: *
266: * -- LAPACK driver routine (version 3.7.0) --
267: * -- LAPACK is a software package provided by Univ. of Tennessee, --
268: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269: * June 2016
270: *
271: * .. Scalar Arguments ..
272: CHARACTER JOBU, JOBVT, RANGE
273: INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
274: DOUBLE PRECISION VL, VU
275: * ..
276: * .. Array Arguments ..
277: INTEGER IWORK( * )
278: DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
279: $ VT( LDVT, * ), WORK( * )
280: * ..
281: *
282: * =====================================================================
283: *
284: * .. Parameters ..
285: DOUBLE PRECISION ZERO, ONE
286: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
287: * ..
288: * .. Local Scalars ..
289: CHARACTER JOBZ, RNGTGK
290: LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
291: INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
292: $ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK,
293: $ J, MAXWRK, MINMN, MINWRK, MNTHR
294: DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
295: * ..
296: * .. Local Arrays ..
297: DOUBLE PRECISION DUM( 1 )
298: * ..
299: * .. External Subroutines ..
300: EXTERNAL DBDSVDX, DGEBRD, DGELQF, DGEQRF, DLACPY,
301: $ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR,
302: $ XERBLA
303: * ..
304: * .. External Functions ..
305: LOGICAL LSAME
306: INTEGER ILAENV
307: DOUBLE PRECISION DLAMCH, DLANGE
308: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
309: * ..
310: * .. Intrinsic Functions ..
311: INTRINSIC MAX, MIN, SQRT
312: * ..
313: * .. Executable Statements ..
314: *
315: * Test the input arguments.
316: *
317: NS = 0
318: INFO = 0
319: ABSTOL = 2*DLAMCH('S')
320: LQUERY = ( LWORK.EQ.-1 )
321: MINMN = MIN( M, N )
322:
323: WANTU = LSAME( JOBU, 'V' )
324: WANTVT = LSAME( JOBVT, 'V' )
325: IF( WANTU .OR. WANTVT ) THEN
326: JOBZ = 'V'
327: ELSE
328: JOBZ = 'N'
329: END IF
330: ALLS = LSAME( RANGE, 'A' )
331: VALS = LSAME( RANGE, 'V' )
332: INDS = LSAME( RANGE, 'I' )
333: *
334: INFO = 0
335: IF( .NOT.LSAME( JOBU, 'V' ) .AND.
336: $ .NOT.LSAME( JOBU, 'N' ) ) THEN
337: INFO = -1
338: ELSE IF( .NOT.LSAME( JOBVT, 'V' ) .AND.
339: $ .NOT.LSAME( JOBVT, 'N' ) ) THEN
340: INFO = -2
341: ELSE IF( .NOT.( ALLS .OR. VALS .OR. INDS ) ) THEN
342: INFO = -3
343: ELSE IF( M.LT.0 ) THEN
344: INFO = -4
345: ELSE IF( N.LT.0 ) THEN
346: INFO = -5
347: ELSE IF( M.GT.LDA ) THEN
348: INFO = -7
349: ELSE IF( MINMN.GT.0 ) THEN
350: IF( VALS ) THEN
351: IF( VL.LT.ZERO ) THEN
352: INFO = -8
353: ELSE IF( VU.LE.VL ) THEN
354: INFO = -9
355: END IF
356: ELSE IF( INDS ) THEN
357: IF( IL.LT.1 .OR. IL.GT.MAX( 1, MINMN ) ) THEN
358: INFO = -10
359: ELSE IF( IU.LT.MIN( MINMN, IL ) .OR. IU.GT.MINMN ) THEN
360: INFO = -11
361: END IF
362: END IF
363: IF( INFO.EQ.0 ) THEN
364: IF( WANTU .AND. LDU.LT.M ) THEN
365: INFO = -15
366: ELSE IF( WANTVT ) THEN
367: IF( INDS ) THEN
368: IF( LDVT.LT.IU-IL+1 ) THEN
369: INFO = -17
370: END IF
371: ELSE IF( LDVT.LT.MINMN ) THEN
372: INFO = -17
373: END IF
374: END IF
375: END IF
376: END IF
377: *
378: * Compute workspace
379: * (Note: Comments in the code beginning "Workspace:" describe the
380: * minimal amount of workspace needed at that point in the code,
381: * as well as the preferred amount for good performance.
382: * NB refers to the optimal block size for the immediately
383: * following subroutine, as returned by ILAENV.)
384: *
385: IF( INFO.EQ.0 ) THEN
386: MINWRK = 1
387: MAXWRK = 1
388: IF( MINMN.GT.0 ) THEN
389: IF( M.GE.N ) THEN
390: MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
391: IF( M.GE.MNTHR ) THEN
392: *
393: * Path 1 (M much larger than N)
394: *
395: MAXWRK = N +
396: $ N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
397: MAXWRK = MAX( MAXWRK, N*(N+5) + 2*N*
398: $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
399: IF (WANTU) THEN
400: MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
401: $ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) )
402: END IF
403: IF (WANTVT) THEN
404: MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
405: $ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) )
406: END IF
407: MINWRK = N*(N*3+20)
408: ELSE
409: *
410: * Path 2 (M at least N, but not much larger)
411: *
412: MAXWRK = 4*N + ( M+N )*
413: $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
414: IF (WANTU) THEN
415: MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
416: $ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) )
417: END IF
418: IF (WANTVT) THEN
419: MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
420: $ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) )
421: END IF
422: MINWRK = MAX(N*(N*2+19),4*N+M)
423: END IF
424: ELSE
425: MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
426: IF( N.GE.MNTHR ) THEN
427: *
428: * Path 1t (N much larger than M)
429: *
430: MAXWRK = M +
431: $ M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
432: MAXWRK = MAX( MAXWRK, M*(M+5) + 2*M*
433: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
434: IF (WANTU) THEN
435: MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
436: $ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) )
437: END IF
438: IF (WANTVT) THEN
439: MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
440: $ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) )
441: END IF
442: MINWRK = M*(M*3+20)
443: ELSE
444: *
445: * Path 2t (N at least M, but not much larger)
446: *
447: MAXWRK = 4*M + ( M+N )*
448: $ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
449: IF (WANTU) THEN
450: MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
451: $ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) )
452: END IF
453: IF (WANTVT) THEN
454: MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
455: $ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) )
456: END IF
457: MINWRK = MAX(M*(M*2+19),4*M+N)
458: END IF
459: END IF
460: END IF
461: MAXWRK = MAX( MAXWRK, MINWRK )
462: WORK( 1 ) = DBLE( MAXWRK )
463: *
464: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
465: INFO = -19
466: END IF
467: END IF
468: *
469: IF( INFO.NE.0 ) THEN
470: CALL XERBLA( 'DGESVDX', -INFO )
471: RETURN
472: ELSE IF( LQUERY ) THEN
473: RETURN
474: END IF
475: *
476: * Quick return if possible
477: *
478: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
479: RETURN
480: END IF
481: *
482: * Set singular values indices accord to RANGE.
483: *
484: IF( ALLS ) THEN
485: RNGTGK = 'I'
486: ILTGK = 1
487: IUTGK = MIN( M, N )
488: ELSE IF( INDS ) THEN
489: RNGTGK = 'I'
490: ILTGK = IL
491: IUTGK = IU
492: ELSE
493: RNGTGK = 'V'
494: ILTGK = 0
495: IUTGK = 0
496: END IF
497: *
498: * Get machine constants
499: *
500: EPS = DLAMCH( 'P' )
501: SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
502: BIGNUM = ONE / SMLNUM
503: *
504: * Scale A if max element outside range [SMLNUM,BIGNUM]
505: *
506: ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
507: ISCL = 0
508: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
509: ISCL = 1
510: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
511: ELSE IF( ANRM.GT.BIGNUM ) THEN
512: ISCL = 1
513: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
514: END IF
515: *
516: IF( M.GE.N ) THEN
517: *
518: * A has at least as many rows as columns. If A has sufficiently
519: * more rows than columns, first reduce A using the QR
520: * decomposition.
521: *
522: IF( M.GE.MNTHR ) THEN
523: *
524: * Path 1 (M much larger than N):
525: * A = Q * R = Q * ( QB * B * PB**T )
526: * = Q * ( QB * ( UB * S * VB**T ) * PB**T )
527: * U = Q * QB * UB; V**T = VB**T * PB**T
528: *
529: * Compute A=Q*R
530: * (Workspace: need 2*N, prefer N+N*NB)
531: *
532: ITAU = 1
533: ITEMP = ITAU + N
534: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),
535: $ LWORK-ITEMP+1, INFO )
536: *
537: * Copy R into WORK and bidiagonalize it:
538: * (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB)
539: *
540: IQRF = ITEMP
541: ID = IQRF + N*N
542: IE = ID + N
543: ITAUQ = IE + N
544: ITAUP = ITAUQ + N
545: ITEMP = ITAUP + N
546: CALL DLACPY( 'U', N, N, A, LDA, WORK( IQRF ), N )
547: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IQRF+1 ), N )
548: CALL DGEBRD( N, N, WORK( IQRF ), N, WORK( ID ), WORK( IE ),
549: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
550: $ LWORK-ITEMP+1, INFO )
551: *
552: * Solve eigenvalue problem TGK*Z=Z*S.
553: * (Workspace: need 14*N + 2*N*(N+1))
554: *
555: ITGKZ = ITEMP
556: ITEMP = ITGKZ + N*(N*2+1)
557: CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
558: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
559: $ N*2, WORK( ITEMP ), IWORK, INFO)
560: *
561: * If needed, compute left singular vectors.
562: *
563: IF( WANTU ) THEN
564: J = ITGKZ
565: DO I = 1, NS
566: CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
567: J = J + N*2
568: END DO
569: CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU )
570: *
571: * Call DORMBR to compute QB*UB.
572: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
573: *
574: CALL DORMBR( 'Q', 'L', 'N', N, NS, N, WORK( IQRF ), N,
575: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
576: $ LWORK-ITEMP+1, INFO )
577: *
578: * Call DORMQR to compute Q*(QB*UB).
579: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
580: *
581: CALL DORMQR( 'L', 'N', M, NS, N, A, LDA,
582: $ WORK( ITAU ), U, LDU, WORK( ITEMP ),
583: $ LWORK-ITEMP+1, INFO )
584: END IF
585: *
586: * If needed, compute right singular vectors.
587: *
588: IF( WANTVT) THEN
589: J = ITGKZ + N
590: DO I = 1, NS
591: CALL DCOPY( N, WORK( J ), 1, VT( I,1 ), LDVT )
592: J = J + N*2
593: END DO
594: *
595: * Call DORMBR to compute VB**T * PB**T
596: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
597: *
598: CALL DORMBR( 'P', 'R', 'T', NS, N, N, WORK( IQRF ), N,
599: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
600: $ LWORK-ITEMP+1, INFO )
601: END IF
602: ELSE
603: *
604: * Path 2 (M at least N, but not much larger)
605: * Reduce A to bidiagonal form without QR decomposition
606: * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
607: * U = QB * UB; V**T = VB**T * PB**T
608: *
609: * Bidiagonalize A
610: * (Workspace: need 4*N+M, prefer 4*N+(M+N)*NB)
611: *
612: ID = 1
613: IE = ID + N
614: ITAUQ = IE + N
615: ITAUP = ITAUQ + N
616: ITEMP = ITAUP + N
617: CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
618: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
619: $ LWORK-ITEMP+1, INFO )
620: *
621: * Solve eigenvalue problem TGK*Z=Z*S.
622: * (Workspace: need 14*N + 2*N*(N+1))
623: *
624: ITGKZ = ITEMP
625: ITEMP = ITGKZ + N*(N*2+1)
626: CALL DBDSVDX( 'U', JOBZ, RNGTGK, N, WORK( ID ), WORK( IE ),
627: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
628: $ N*2, WORK( ITEMP ), IWORK, INFO)
629: *
630: * If needed, compute left singular vectors.
631: *
632: IF( WANTU ) THEN
633: J = ITGKZ
634: DO I = 1, NS
635: CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
636: J = J + N*2
637: END DO
638: CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU )
639: *
640: * Call DORMBR to compute QB*UB.
641: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
642: *
643: CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
644: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
645: $ LWORK-ITEMP+1, IERR )
646: END IF
647: *
648: * If needed, compute right singular vectors.
649: *
650: IF( WANTVT) THEN
651: J = ITGKZ + N
652: DO I = 1, NS
653: CALL DCOPY( N, WORK( J ), 1, VT( I,1 ), LDVT )
654: J = J + N*2
655: END DO
656: *
657: * Call DORMBR to compute VB**T * PB**T
658: * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
659: *
660: CALL DORMBR( 'P', 'R', 'T', NS, N, N, A, LDA,
661: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
662: $ LWORK-ITEMP+1, IERR )
663: END IF
664: END IF
665: ELSE
666: *
667: * A has more columns than rows. If A has sufficiently more
668: * columns than rows, first reduce A using the LQ decomposition.
669: *
670: IF( N.GE.MNTHR ) THEN
671: *
672: * Path 1t (N much larger than M):
673: * A = L * Q = ( QB * B * PB**T ) * Q
674: * = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
675: * U = QB * UB ; V**T = VB**T * PB**T * Q
676: *
677: * Compute A=L*Q
678: * (Workspace: need 2*M, prefer M+M*NB)
679: *
680: ITAU = 1
681: ITEMP = ITAU + M
682: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( ITEMP ),
683: $ LWORK-ITEMP+1, INFO )
684:
685: * Copy L into WORK and bidiagonalize it:
686: * (Workspace in WORK( ITEMP ): need M*M+5*N, prefer M*M+4*M+2*M*NB)
687: *
688: ILQF = ITEMP
689: ID = ILQF + M*M
690: IE = ID + M
691: ITAUQ = IE + M
692: ITAUP = ITAUQ + M
693: ITEMP = ITAUP + M
694: CALL DLACPY( 'L', M, M, A, LDA, WORK( ILQF ), M )
695: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( ILQF+M ), M )
696: CALL DGEBRD( M, M, WORK( ILQF ), M, WORK( ID ), WORK( IE ),
697: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
698: $ LWORK-ITEMP+1, INFO )
699: *
700: * Solve eigenvalue problem TGK*Z=Z*S.
701: * (Workspace: need 2*M*M+14*M)
702: *
703: ITGKZ = ITEMP
704: ITEMP = ITGKZ + M*(M*2+1)
705: CALL DBDSVDX( 'U', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
706: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
707: $ M*2, WORK( ITEMP ), IWORK, INFO)
708: *
709: * If needed, compute left singular vectors.
710: *
711: IF( WANTU ) THEN
712: J = ITGKZ
713: DO I = 1, NS
714: CALL DCOPY( M, WORK( J ), 1, U( 1,I ), 1 )
715: J = J + M*2
716: END DO
717: *
718: * Call DORMBR to compute QB*UB.
719: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
720: *
721: CALL DORMBR( 'Q', 'L', 'N', M, NS, M, WORK( ILQF ), M,
722: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
723: $ LWORK-ITEMP+1, INFO )
724: END IF
725: *
726: * If needed, compute right singular vectors.
727: *
728: IF( WANTVT) THEN
729: J = ITGKZ + M
730: DO I = 1, NS
731: CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
732: J = J + M*2
733: END DO
734: CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT)
735: *
736: * Call DORMBR to compute (VB**T)*(PB**T)
737: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
738: *
739: CALL DORMBR( 'P', 'R', 'T', NS, M, M, WORK( ILQF ), M,
740: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
741: $ LWORK-ITEMP+1, INFO )
742: *
743: * Call DORMLQ to compute ((VB**T)*(PB**T))*Q.
744: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
745: *
746: CALL DORMLQ( 'R', 'N', NS, N, M, A, LDA,
747: $ WORK( ITAU ), VT, LDVT, WORK( ITEMP ),
748: $ LWORK-ITEMP+1, INFO )
749: END IF
750: ELSE
751: *
752: * Path 2t (N greater than M, but not much larger)
753: * Reduce to bidiagonal form without LQ decomposition
754: * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
755: * U = QB * UB; V**T = VB**T * PB**T
756: *
757: * Bidiagonalize A
758: * (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB)
759: *
760: ID = 1
761: IE = ID + M
762: ITAUQ = IE + M
763: ITAUP = ITAUQ + M
764: ITEMP = ITAUP + M
765: CALL DGEBRD( M, N, A, LDA, WORK( ID ), WORK( IE ),
766: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
767: $ LWORK-ITEMP+1, INFO )
768: *
769: * Solve eigenvalue problem TGK*Z=Z*S.
770: * (Workspace: need 2*M*M+14*M)
771: *
772: ITGKZ = ITEMP
773: ITEMP = ITGKZ + M*(M*2+1)
774: CALL DBDSVDX( 'L', JOBZ, RNGTGK, M, WORK( ID ), WORK( IE ),
775: $ VL, VU, ILTGK, IUTGK, NS, S, WORK( ITGKZ ),
776: $ M*2, WORK( ITEMP ), IWORK, INFO)
777: *
778: * If needed, compute left singular vectors.
779: *
780: IF( WANTU ) THEN
781: J = ITGKZ
782: DO I = 1, NS
783: CALL DCOPY( M, WORK( J ), 1, U( 1,I ), 1 )
784: J = J + M*2
785: END DO
786: *
787: * Call DORMBR to compute QB*UB.
788: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
789: *
790: CALL DORMBR( 'Q', 'L', 'N', M, NS, N, A, LDA,
791: $ WORK( ITAUQ ), U, LDU, WORK( ITEMP ),
792: $ LWORK-ITEMP+1, INFO )
793: END IF
794: *
795: * If needed, compute right singular vectors.
796: *
797: IF( WANTVT) THEN
798: J = ITGKZ + M
799: DO I = 1, NS
800: CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
801: J = J + M*2
802: END DO
803: CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT)
804: *
805: * Call DORMBR to compute VB**T * PB**T
806: * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
807: *
808: CALL DORMBR( 'P', 'R', 'T', NS, N, M, A, LDA,
809: $ WORK( ITAUP ), VT, LDVT, WORK( ITEMP ),
810: $ LWORK-ITEMP+1, INFO )
811: END IF
812: END IF
813: END IF
814: *
815: * Undo scaling if necessary
816: *
817: IF( ISCL.EQ.1 ) THEN
818: IF( ANRM.GT.BIGNUM )
819: $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1,
820: $ S, MINMN, INFO )
821: IF( ANRM.LT.SMLNUM )
822: $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1,
823: $ S, MINMN, INFO )
824: END IF
825: *
826: * Return optimal workspace in WORK(1)
827: *
828: WORK( 1 ) = DBLE( MAXWRK )
829: *
830: RETURN
831: *
832: * End of DGESVDX
833: *
834: 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