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