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