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