Annotation of rpl/lapack/lapack/zgesvdq.f, revision 1.2
1.1 bertrand 1: *> \brief <b> ZGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method 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 ZGESVDQ + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvdq.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvdq.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvdq.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
22: * S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
23: * CWORK, LCWORK, RWORK, LRWORK, INFO )
24: *
25: * .. Scalar Arguments ..
26: * IMPLICIT NONE
27: * CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV
28: * INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LCWORK, LRWORK,
29: * INFO
30: * ..
31: * .. Array Arguments ..
32: * COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * )
33: * DOUBLE PRECISION S( * ), RWORK( * )
34: * INTEGER IWORK( * )
35: * ..
36: *
37: *
38: *> \par Purpose:
39: * =============
40: *>
41: *> \verbatim
42: *>
1.2 ! bertrand 43: *> ZCGESVDQ computes the singular value decomposition (SVD) of a complex
1.1 bertrand 44: *> M-by-N matrix A, where M >= N. The SVD of A is written as
45: *> [++] [xx] [x0] [xx]
46: *> A = U * SIGMA * V^*, [++] = [xx] * [ox] * [xx]
47: *> [++] [xx]
48: *> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
49: *> matrix, and V is an N-by-N unitary matrix. The diagonal elements
50: *> of SIGMA are the singular values of A. The columns of U and V are the
51: *> left and the right singular vectors of A, respectively.
52: *> \endverbatim
53: *
54: * Arguments
55: * =========
56: *
57: *> \param[in] JOBA
58: *> \verbatim
59: *> JOBA is CHARACTER*1
60: *> Specifies the level of accuracy in the computed SVD
61: *> = 'A' The requested accuracy corresponds to having the backward
62: *> error bounded by || delta A ||_F <= f(m,n) * EPS * || A ||_F,
63: *> where EPS = DLAMCH('Epsilon'). This authorises ZGESVDQ to
64: *> truncate the computed triangular factor in a rank revealing
65: *> QR factorization whenever the truncated part is below the
66: *> threshold of the order of EPS * ||A||_F. This is aggressive
67: *> truncation level.
68: *> = 'M' Similarly as with 'A', but the truncation is more gentle: it
69: *> is allowed only when there is a drop on the diagonal of the
70: *> triangular factor in the QR factorization. This is medium
71: *> truncation level.
72: *> = 'H' High accuracy requested. No numerical rank determination based
73: *> on the rank revealing QR factorization is attempted.
74: *> = 'E' Same as 'H', and in addition the condition number of column
75: *> scaled A is estimated and returned in RWORK(1).
76: *> N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1)
77: *> \endverbatim
78: *>
79: *> \param[in] JOBP
80: *> \verbatim
81: *> JOBP is CHARACTER*1
82: *> = 'P' The rows of A are ordered in decreasing order with respect to
83: *> ||A(i,:)||_\infty. This enhances numerical accuracy at the cost
84: *> of extra data movement. Recommended for numerical robustness.
85: *> = 'N' No row pivoting.
86: *> \endverbatim
87: *>
88: *> \param[in] JOBR
89: *> \verbatim
90: *> JOBR is CHARACTER*1
91: *> = 'T' After the initial pivoted QR factorization, ZGESVD is applied to
92: *> the adjoint R**H of the computed triangular factor R. This involves
93: *> some extra data movement (matrix transpositions). Useful for
94: *> experiments, research and development.
95: *> = 'N' The triangular factor R is given as input to CGESVD. This may be
96: *> preferred as it involves less data movement.
97: *> \endverbatim
98: *>
99: *> \param[in] JOBU
100: *> \verbatim
101: *> JOBU is CHARACTER*1
102: *> = 'A' All M left singular vectors are computed and returned in the
103: *> matrix U. See the description of U.
104: *> = 'S' or 'U' N = min(M,N) left singular vectors are computed and returned
105: *> in the matrix U. See the description of U.
106: *> = 'R' Numerical rank NUMRANK is determined and only NUMRANK left singular
107: *> vectors are computed and returned in the matrix U.
108: *> = 'F' The N left singular vectors are returned in factored form as the
109: *> product of the Q factor from the initial QR factorization and the
110: *> N left singular vectors of (R**H , 0)**H. If row pivoting is used,
111: *> then the necessary information on the row pivoting is stored in
112: *> IWORK(N+1:N+M-1).
113: *> = 'N' The left singular vectors are not computed.
114: *> \endverbatim
115: *>
116: *> \param[in] JOBV
117: *> \verbatim
118: *> JOBV is CHARACTER*1
119: *> = 'A', 'V' All N right singular vectors are computed and returned in
120: *> the matrix V.
121: *> = 'R' Numerical rank NUMRANK is determined and only NUMRANK right singular
122: *> vectors are computed and returned in the matrix V. This option is
123: *> allowed only if JOBU = 'R' or JOBU = 'N'; otherwise it is illegal.
124: *> = 'N' The right singular vectors are not computed.
125: *> \endverbatim
126: *>
127: *> \param[in] M
128: *> \verbatim
129: *> M is INTEGER
130: *> The number of rows of the input matrix A. M >= 0.
131: *> \endverbatim
132: *>
133: *> \param[in] N
134: *> \verbatim
135: *> N is INTEGER
136: *> The number of columns of the input matrix A. M >= N >= 0.
137: *> \endverbatim
138: *>
139: *> \param[in,out] A
140: *> \verbatim
141: *> A is COMPLEX*16 array of dimensions LDA x N
142: *> On entry, the input matrix A.
143: *> On exit, if JOBU .NE. 'N' or JOBV .NE. 'N', the lower triangle of A contains
144: *> the Householder vectors as stored by ZGEQP3. If JOBU = 'F', these Householder
145: *> vectors together with CWORK(1:N) can be used to restore the Q factors from
146: *> the initial pivoted QR factorization of A. See the description of U.
147: *> \endverbatim
148: *>
149: *> \param[in] LDA
150: *> \verbatim
151: *> LDA is INTEGER.
152: *> The leading dimension of the array A. LDA >= max(1,M).
153: *> \endverbatim
154: *>
155: *> \param[out] S
156: *> \verbatim
157: *> S is DOUBLE PRECISION array of dimension N.
158: *> The singular values of A, ordered so that S(i) >= S(i+1).
159: *> \endverbatim
160: *>
161: *> \param[out] U
162: *> \verbatim
163: *> U is COMPLEX*16 array, dimension
164: *> LDU x M if JOBU = 'A'; see the description of LDU. In this case,
165: *> on exit, U contains the M left singular vectors.
166: *> LDU x N if JOBU = 'S', 'U', 'R' ; see the description of LDU. In this
167: *> case, U contains the leading N or the leading NUMRANK left singular vectors.
168: *> LDU x N if JOBU = 'F' ; see the description of LDU. In this case U
169: *> contains N x N unitary matrix that can be used to form the left
170: *> singular vectors.
171: *> If JOBU = 'N', U is not referenced.
172: *> \endverbatim
173: *>
174: *> \param[in] LDU
175: *> \verbatim
176: *> LDU is INTEGER.
177: *> The leading dimension of the array U.
178: *> If JOBU = 'A', 'S', 'U', 'R', LDU >= max(1,M).
179: *> If JOBU = 'F', LDU >= max(1,N).
180: *> Otherwise, LDU >= 1.
181: *> \endverbatim
182: *>
183: *> \param[out] V
184: *> \verbatim
185: *> V is COMPLEX*16 array, dimension
186: *> LDV x N if JOBV = 'A', 'V', 'R' or if JOBA = 'E' .
187: *> If JOBV = 'A', or 'V', V contains the N-by-N unitary matrix V**H;
188: *> If JOBV = 'R', V contains the first NUMRANK rows of V**H (the right
189: *> singular vectors, stored rowwise, of the NUMRANK largest singular values).
190: *> If JOBV = 'N' and JOBA = 'E', V is used as a workspace.
191: *> If JOBV = 'N', and JOBA.NE.'E', V is not referenced.
192: *> \endverbatim
193: *>
194: *> \param[in] LDV
195: *> \verbatim
196: *> LDV is INTEGER
197: *> The leading dimension of the array V.
198: *> If JOBV = 'A', 'V', 'R', or JOBA = 'E', LDV >= max(1,N).
199: *> Otherwise, LDV >= 1.
200: *> \endverbatim
201: *>
202: *> \param[out] NUMRANK
203: *> \verbatim
204: *> NUMRANK is INTEGER
205: *> NUMRANK is the numerical rank first determined after the rank
206: *> revealing QR factorization, following the strategy specified by the
207: *> value of JOBA. If JOBV = 'R' and JOBU = 'R', only NUMRANK
208: *> leading singular values and vectors are then requested in the call
209: *> of CGESVD. The final value of NUMRANK might be further reduced if
210: *> some singular values are computed as zeros.
211: *> \endverbatim
212: *>
213: *> \param[out] IWORK
214: *> \verbatim
215: *> IWORK is INTEGER array, dimension (max(1, LIWORK)).
216: *> On exit, IWORK(1:N) contains column pivoting permutation of the
217: *> rank revealing QR factorization.
218: *> If JOBP = 'P', IWORK(N+1:N+M-1) contains the indices of the sequence
219: *> of row swaps used in row pivoting. These can be used to restore the
220: *> left singular vectors in the case JOBU = 'F'.
1.2 ! bertrand 221: *>
1.1 bertrand 222: *> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
1.2 ! bertrand 223: *> IWORK(1) returns the minimal LIWORK.
1.1 bertrand 224: *> \endverbatim
225: *>
226: *> \param[in] LIWORK
227: *> \verbatim
228: *> LIWORK is INTEGER
229: *> The dimension of the array IWORK.
230: *> LIWORK >= N + M - 1, if JOBP = 'P';
231: *> LIWORK >= N if JOBP = 'N'.
232: *>
233: *> If LIWORK = -1, then a workspace query is assumed; the routine
234: *> only calculates and returns the optimal and minimal sizes
235: *> for the CWORK, IWORK, and RWORK arrays, and no error
236: *> message related to LCWORK is issued by XERBLA.
237: *> \endverbatim
238: *>
239: *> \param[out] CWORK
240: *> \verbatim
241: *> CWORK is COMPLEX*12 array, dimension (max(2, LCWORK)), used as a workspace.
242: *> On exit, if, on entry, LCWORK.NE.-1, CWORK(1:N) contains parameters
243: *> needed to recover the Q factor from the QR factorization computed by
244: *> ZGEQP3.
1.2 ! bertrand 245: *>
1.1 bertrand 246: *> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
247: *> CWORK(1) returns the optimal LCWORK, and
248: *> CWORK(2) returns the minimal LCWORK.
249: *> \endverbatim
250: *>
251: *> \param[in,out] LCWORK
252: *> \verbatim
253: *> LCWORK is INTEGER
254: *> The dimension of the array CWORK. It is determined as follows:
255: *> Let LWQP3 = N+1, LWCON = 2*N, and let
256: *> LWUNQ = { MAX( N, 1 ), if JOBU = 'R', 'S', or 'U'
257: *> { MAX( M, 1 ), if JOBU = 'A'
258: *> LWSVD = MAX( 3*N, 1 )
259: *> LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 3*(N/2), 1 ), LWUNLQ = MAX( N, 1 ),
260: *> LWQRF = MAX( N/2, 1 ), LWUNQ2 = MAX( N, 1 )
261: *> Then the minimal value of LCWORK is:
262: *> = MAX( N + LWQP3, LWSVD ) if only the singular values are needed;
263: *> = MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values are needed,
264: *> and a scaled condition estimate requested;
265: *>
266: *> = N + MAX( LWQP3, LWSVD, LWUNQ ) if the singular values and the left
267: *> singular vectors are requested;
268: *> = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the singular values and the left
269: *> singular vectors are requested, and also
270: *> a scaled condition estimate requested;
271: *>
272: *> = N + MAX( LWQP3, LWSVD ) if the singular values and the right
273: *> singular vectors are requested;
274: *> = N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and the right
275: *> singular vectors are requested, and also
276: *> a scaled condition etimate requested;
277: *>
278: *> = N + MAX( LWQP3, LWSVD, LWUNQ ) if the full SVD is requested with JOBV = 'R';
279: *> independent of JOBR;
280: *> = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the full SVD is requested,
281: *> JOBV = 'R' and, also a scaled condition
282: *> estimate requested; independent of JOBR;
283: *> = MAX( N + MAX( LWQP3, LWSVD, LWUNQ ),
284: *> N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ) ) if the
285: *> full SVD is requested with JOBV = 'A' or 'V', and
286: *> JOBR ='N'
287: *> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
288: *> N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ ) )
289: *> if the full SVD is requested with JOBV = 'A' or 'V', and
290: *> JOBR ='N', and also a scaled condition number estimate
291: *> requested.
292: *> = MAX( N + MAX( LWQP3, LWSVD, LWUNQ ),
293: *> N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) ) if the
294: *> full SVD is requested with JOBV = 'A', 'V', and JOBR ='T'
295: *> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
296: *> N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) )
297: *> if the full SVD is requested with JOBV = 'A', 'V' and
298: *> JOBR ='T', and also a scaled condition number estimate
299: *> requested.
300: *> Finally, LCWORK must be at least two: LCWORK = MAX( 2, LCWORK ).
301: *>
302: *> If LCWORK = -1, then a workspace query is assumed; the routine
303: *> only calculates and returns the optimal and minimal sizes
304: *> for the CWORK, IWORK, and RWORK arrays, and no error
305: *> message related to LCWORK is issued by XERBLA.
306: *> \endverbatim
307: *>
308: *> \param[out] RWORK
309: *> \verbatim
310: *> RWORK is DOUBLE PRECISION array, dimension (max(1, LRWORK)).
311: *> On exit,
312: *> 1. If JOBA = 'E', RWORK(1) contains an estimate of the condition
313: *> number of column scaled A. If A = C * D where D is diagonal and C
314: *> has unit columns in the Euclidean norm, then, assuming full column rank,
315: *> N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1).
316: *> Otherwise, RWORK(1) = -1.
317: *> 2. RWORK(2) contains the number of singular values computed as
1.2 ! bertrand 318: *> exact zeros in ZGESVD applied to the upper triangular or trapezoidal
1.1 bertrand 319: *> R (from the initial QR factorization). In case of early exit (no call to
320: *> ZGESVD, such as in the case of zero matrix) RWORK(2) = -1.
1.2 ! bertrand 321: *>
1.1 bertrand 322: *> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
323: *> RWORK(1) returns the minimal LRWORK.
324: *> \endverbatim
325: *>
326: *> \param[in] LRWORK
327: *> \verbatim
328: *> LRWORK is INTEGER.
329: *> The dimension of the array RWORK.
330: *> If JOBP ='P', then LRWORK >= MAX(2, M, 5*N);
331: *> Otherwise, LRWORK >= MAX(2, 5*N).
1.2 ! bertrand 332: *>
1.1 bertrand 333: *> If LRWORK = -1, then a workspace query is assumed; the routine
334: *> only calculates and returns the optimal and minimal sizes
335: *> for the CWORK, IWORK, and RWORK arrays, and no error
336: *> message related to LCWORK is issued by XERBLA.
337: *> \endverbatim
338: *>
339: *> \param[out] INFO
340: *> \verbatim
341: *> INFO is INTEGER
342: *> = 0: successful exit.
343: *> < 0: if INFO = -i, the i-th argument had an illegal value.
344: *> > 0: if ZBDSQR did not converge, INFO specifies how many superdiagonals
345: *> of an intermediate bidiagonal form B (computed in ZGESVD) did not
346: *> converge to zero.
347: *> \endverbatim
348: *
349: *> \par Further Details:
350: * ========================
351: *>
352: *> \verbatim
353: *>
354: *> 1. The data movement (matrix transpose) is coded using simple nested
355: *> DO-loops because BLAS and LAPACK do not provide corresponding subroutines.
356: *> Those DO-loops are easily identified in this source code - by the CONTINUE
357: *> statements labeled with 11**. In an optimized version of this code, the
358: *> nested DO loops should be replaced with calls to an optimized subroutine.
359: *> 2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause
360: *> column norm overflow. This is the minial precaution and it is left to the
361: *> SVD routine (CGESVD) to do its own preemptive scaling if potential over-
362: *> or underflows are detected. To avoid repeated scanning of the array A,
363: *> an optimal implementation would do all necessary scaling before calling
364: *> CGESVD and the scaling in CGESVD can be switched off.
365: *> 3. Other comments related to code optimization are given in comments in the
366: *> code, enlosed in [[double brackets]].
367: *> \endverbatim
368: *
369: *> \par Bugs, examples and comments
370: * ===========================
371: *
372: *> \verbatim
373: *> Please report all bugs and send interesting examples and/or comments to
374: *> drmac@math.hr. Thank you.
375: *> \endverbatim
376: *
377: *> \par References
378: * ===============
379: *
380: *> \verbatim
381: *> [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
382: *> Computing the SVD with High Accuracy. ACM Trans. Math. Softw.
383: *> 44(1): 11:1-11:30 (2017)
384: *>
385: *> SIGMA library, xGESVDQ section updated February 2016.
386: *> Developed and coded by Zlatko Drmac, Department of Mathematics
387: *> University of Zagreb, Croatia, drmac@math.hr
388: *> \endverbatim
389: *
390: *
391: *> \par Contributors:
392: * ==================
393: *>
394: *> \verbatim
395: *> Developed and coded by Zlatko Drmac, Department of Mathematics
396: *> University of Zagreb, Croatia, drmac@math.hr
397: *> \endverbatim
398: *
399: * Authors:
400: * ========
401: *
402: *> \author Univ. of Tennessee
403: *> \author Univ. of California Berkeley
404: *> \author Univ. of Colorado Denver
405: *> \author NAG Ltd.
406: *
407: *> \ingroup complex16GEsing
408: *
409: * =====================================================================
410: SUBROUTINE ZGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
411: $ S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
412: $ CWORK, LCWORK, RWORK, LRWORK, INFO )
413: * .. Scalar Arguments ..
414: IMPLICIT NONE
415: CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV
416: INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LCWORK, LRWORK,
417: $ INFO
418: * ..
419: * .. Array Arguments ..
420: COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * )
421: DOUBLE PRECISION S( * ), RWORK( * )
422: INTEGER IWORK( * )
423: *
424: * =====================================================================
425: *
426: * .. Parameters ..
427: DOUBLE PRECISION ZERO, ONE
428: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
429: COMPLEX*16 CZERO, CONE
430: PARAMETER ( CZERO = (0.0D0,0.0D0), CONE = (1.0D0,0.0D0) )
431: * ..
432: * .. Local Scalars ..
433: INTEGER IERR, NR, N1, OPTRATIO, p, q
434: INTEGER LWCON, LWQP3, LWRK_ZGELQF, LWRK_ZGESVD, LWRK_ZGESVD2,
435: $ LWRK_ZGEQP3, LWRK_ZGEQRF, LWRK_ZUNMLQ, LWRK_ZUNMQR,
436: $ LWRK_ZUNMQR2, LWLQF, LWQRF, LWSVD, LWSVD2, LWUNQ,
437: $ LWUNQ2, LWUNLQ, MINWRK, MINWRK2, OPTWRK, OPTWRK2,
438: $ IMINWRK, RMINWRK
439: LOGICAL ACCLA, ACCLM, ACCLH, ASCALED, CONDA, DNTWU, DNTWV,
440: $ LQUERY, LSVC0, LSVEC, ROWPRM, RSVEC, RTRANS, WNTUA,
441: $ WNTUF, WNTUR, WNTUS, WNTVA, WNTVR
442: DOUBLE PRECISION BIG, EPSLN, RTMP, SCONDA, SFMIN
443: COMPLEX*16 CTMP
444: * ..
445: * .. Local Arrays
446: COMPLEX*16 CDUMMY(1)
447: DOUBLE PRECISION RDUMMY(1)
448: * ..
449: * .. External Subroutines (BLAS, LAPACK)
450: EXTERNAL ZGELQF, ZGEQP3, ZGEQRF, ZGESVD, ZLACPY, ZLAPMT,
451: $ ZLASCL, ZLASET, ZLASWP, ZDSCAL, DLASET, DLASCL,
452: $ ZPOCON, ZUNMLQ, ZUNMQR, XERBLA
453: * ..
454: * .. External Functions (BLAS, LAPACK)
455: LOGICAL LSAME
456: INTEGER IDAMAX
457: DOUBLE PRECISION ZLANGE, DZNRM2, DLAMCH
458: EXTERNAL LSAME, ZLANGE, IDAMAX, DZNRM2, DLAMCH
459: * ..
460: * .. Intrinsic Functions ..
461: INTRINSIC ABS, CONJG, MAX, MIN, DBLE, SQRT
462: * ..
463: * .. Executable Statements ..
464: *
465: * Test the input arguments
466: *
467: WNTUS = LSAME( JOBU, 'S' ) .OR. LSAME( JOBU, 'U' )
468: WNTUR = LSAME( JOBU, 'R' )
469: WNTUA = LSAME( JOBU, 'A' )
470: WNTUF = LSAME( JOBU, 'F' )
471: LSVC0 = WNTUS .OR. WNTUR .OR. WNTUA
472: LSVEC = LSVC0 .OR. WNTUF
473: DNTWU = LSAME( JOBU, 'N' )
474: *
475: WNTVR = LSAME( JOBV, 'R' )
476: WNTVA = LSAME( JOBV, 'A' ) .OR. LSAME( JOBV, 'V' )
477: RSVEC = WNTVR .OR. WNTVA
478: DNTWV = LSAME( JOBV, 'N' )
479: *
480: ACCLA = LSAME( JOBA, 'A' )
481: ACCLM = LSAME( JOBA, 'M' )
482: CONDA = LSAME( JOBA, 'E' )
483: ACCLH = LSAME( JOBA, 'H' ) .OR. CONDA
484: *
485: ROWPRM = LSAME( JOBP, 'P' )
486: RTRANS = LSAME( JOBR, 'T' )
487: *
488: IF ( ROWPRM ) THEN
489: IMINWRK = MAX( 1, N + M - 1 )
490: RMINWRK = MAX( 2, M, 5*N )
491: ELSE
492: IMINWRK = MAX( 1, N )
493: RMINWRK = MAX( 2, 5*N )
494: END IF
495: LQUERY = (LIWORK .EQ. -1 .OR. LCWORK .EQ. -1 .OR. LRWORK .EQ. -1)
496: INFO = 0
497: IF ( .NOT. ( ACCLA .OR. ACCLM .OR. ACCLH ) ) THEN
498: INFO = -1
499: ELSE IF ( .NOT.( ROWPRM .OR. LSAME( JOBP, 'N' ) ) ) THEN
500: INFO = -2
501: ELSE IF ( .NOT.( RTRANS .OR. LSAME( JOBR, 'N' ) ) ) THEN
502: INFO = -3
503: ELSE IF ( .NOT.( LSVEC .OR. DNTWU ) ) THEN
504: INFO = -4
505: ELSE IF ( WNTUR .AND. WNTVA ) THEN
506: INFO = -5
507: ELSE IF ( .NOT.( RSVEC .OR. DNTWV )) THEN
508: INFO = -5
509: ELSE IF ( M.LT.0 ) THEN
510: INFO = -6
511: ELSE IF ( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
512: INFO = -7
513: ELSE IF ( LDA.LT.MAX( 1, M ) ) THEN
514: INFO = -9
515: ELSE IF ( LDU.LT.1 .OR. ( LSVC0 .AND. LDU.LT.M ) .OR.
516: $ ( WNTUF .AND. LDU.LT.N ) ) THEN
517: INFO = -12
518: ELSE IF ( LDV.LT.1 .OR. ( RSVEC .AND. LDV.LT.N ) .OR.
519: $ ( CONDA .AND. LDV.LT.N ) ) THEN
520: INFO = -14
521: ELSE IF ( LIWORK .LT. IMINWRK .AND. .NOT. LQUERY ) THEN
522: INFO = -17
523: END IF
524: *
525: *
526: IF ( INFO .EQ. 0 ) THEN
527: * .. compute the minimal and the optimal workspace lengths
528: * [[The expressions for computing the minimal and the optimal
529: * values of LCWORK are written with a lot of redundancy and
530: * can be simplified. However, this detailed form is easier for
531: * maintenance and modifications of the code.]]
532: *
533: * .. minimal workspace length for ZGEQP3 of an M x N matrix
534: LWQP3 = N+1
535: * .. minimal workspace length for ZUNMQR to build left singular vectors
536: IF ( WNTUS .OR. WNTUR ) THEN
537: LWUNQ = MAX( N , 1 )
538: ELSE IF ( WNTUA ) THEN
539: LWUNQ = MAX( M , 1 )
540: END IF
541: * .. minimal workspace length for ZPOCON of an N x N matrix
542: LWCON = 2 * N
543: * .. ZGESVD of an N x N matrix
544: LWSVD = MAX( 3 * N, 1 )
545: IF ( LQUERY ) THEN
546: CALL ZGEQP3( M, N, A, LDA, IWORK, CDUMMY, CDUMMY, -1,
547: $ RDUMMY, IERR )
548: LWRK_ZGEQP3 = INT( CDUMMY(1) )
549: IF ( WNTUS .OR. WNTUR ) THEN
550: CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
551: $ LDU, CDUMMY, -1, IERR )
552: LWRK_ZUNMQR = INT( CDUMMY(1) )
553: ELSE IF ( WNTUA ) THEN
554: CALL ZUNMQR( 'L', 'N', M, M, N, A, LDA, CDUMMY, U,
555: $ LDU, CDUMMY, -1, IERR )
556: LWRK_ZUNMQR = INT( CDUMMY(1) )
557: ELSE
558: LWRK_ZUNMQR = 0
559: END IF
560: END IF
561: MINWRK = 2
562: OPTWRK = 2
563: IF ( .NOT. (LSVEC .OR. RSVEC ) ) THEN
564: * .. minimal and optimal sizes of the complex workspace if
565: * only the singular values are requested
566: IF ( CONDA ) THEN
567: MINWRK = MAX( N+LWQP3, LWCON, LWSVD )
568: ELSE
569: MINWRK = MAX( N+LWQP3, LWSVD )
570: END IF
571: IF ( LQUERY ) THEN
572: CALL ZGESVD( 'N', 'N', N, N, A, LDA, S, U, LDU,
573: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
574: LWRK_ZGESVD = INT( CDUMMY(1) )
575: IF ( CONDA ) THEN
576: OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON, LWRK_ZGESVD )
577: ELSE
578: OPTWRK = MAX( N+LWRK_ZGEQP3, LWRK_ZGESVD )
579: END IF
580: END IF
581: ELSE IF ( LSVEC .AND. (.NOT.RSVEC) ) THEN
582: * .. minimal and optimal sizes of the complex workspace if the
583: * singular values and the left singular vectors are requested
584: IF ( CONDA ) THEN
585: MINWRK = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ )
586: ELSE
587: MINWRK = N + MAX( LWQP3, LWSVD, LWUNQ )
588: END IF
589: IF ( LQUERY ) THEN
590: IF ( RTRANS ) THEN
591: CALL ZGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU,
592: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
593: ELSE
594: CALL ZGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU,
595: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
596: END IF
597: LWRK_ZGESVD = INT( CDUMMY(1) )
598: IF ( CONDA ) THEN
599: OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, LWRK_ZGESVD,
600: $ LWRK_ZUNMQR )
601: ELSE
602: OPTWRK = N + MAX( LWRK_ZGEQP3, LWRK_ZGESVD,
603: $ LWRK_ZUNMQR )
604: END IF
605: END IF
606: ELSE IF ( RSVEC .AND. (.NOT.LSVEC) ) THEN
607: * .. minimal and optimal sizes of the complex workspace if the
608: * singular values and the right singular vectors are requested
609: IF ( CONDA ) THEN
610: MINWRK = N + MAX( LWQP3, LWCON, LWSVD )
611: ELSE
612: MINWRK = N + MAX( LWQP3, LWSVD )
613: END IF
614: IF ( LQUERY ) THEN
615: IF ( RTRANS ) THEN
616: CALL ZGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU,
617: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
618: ELSE
619: CALL ZGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU,
620: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
621: END IF
622: LWRK_ZGESVD = INT( CDUMMY(1) )
623: IF ( CONDA ) THEN
624: OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, LWRK_ZGESVD )
625: ELSE
626: OPTWRK = N + MAX( LWRK_ZGEQP3, LWRK_ZGESVD )
627: END IF
628: END IF
629: ELSE
630: * .. minimal and optimal sizes of the complex workspace if the
631: * full SVD is requested
632: IF ( RTRANS ) THEN
633: MINWRK = MAX( LWQP3, LWSVD, LWUNQ )
634: IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON )
635: MINWRK = MINWRK + N
636: IF ( WNTVA ) THEN
637: * .. minimal workspace length for N x N/2 ZGEQRF
638: LWQRF = MAX( N/2, 1 )
1.2 ! bertrand 639: * .. minimal workspace length for N/2 x N/2 ZGESVD
1.1 bertrand 640: LWSVD2 = MAX( 3 * (N/2), 1 )
641: LWUNQ2 = MAX( N, 1 )
642: MINWRK2 = MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2,
643: $ N/2+LWUNQ2, LWUNQ )
644: IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON )
645: MINWRK2 = N + MINWRK2
646: MINWRK = MAX( MINWRK, MINWRK2 )
647: END IF
648: ELSE
649: MINWRK = MAX( LWQP3, LWSVD, LWUNQ )
650: IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON )
651: MINWRK = MINWRK + N
652: IF ( WNTVA ) THEN
653: * .. minimal workspace length for N/2 x N ZGELQF
654: LWLQF = MAX( N/2, 1 )
655: LWSVD2 = MAX( 3 * (N/2), 1 )
656: LWUNLQ = MAX( N , 1 )
657: MINWRK2 = MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2,
658: $ N/2+LWUNLQ, LWUNQ )
659: IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON )
660: MINWRK2 = N + MINWRK2
661: MINWRK = MAX( MINWRK, MINWRK2 )
662: END IF
663: END IF
664: IF ( LQUERY ) THEN
665: IF ( RTRANS ) THEN
666: CALL ZGESVD( 'O', 'A', N, N, A, LDA, S, U, LDU,
667: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
668: LWRK_ZGESVD = INT( CDUMMY(1) )
669: OPTWRK = MAX(LWRK_ZGEQP3,LWRK_ZGESVD,LWRK_ZUNMQR)
670: IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON )
671: OPTWRK = N + OPTWRK
672: IF ( WNTVA ) THEN
673: CALL ZGEQRF(N,N/2,U,LDU,CDUMMY,CDUMMY,-1,IERR)
674: LWRK_ZGEQRF = INT( CDUMMY(1) )
675: CALL ZGESVD( 'S', 'O', N/2,N/2, V,LDV, S, U,LDU,
676: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
677: LWRK_ZGESVD2 = INT( CDUMMY(1) )
678: CALL ZUNMQR( 'R', 'C', N, N, N/2, U, LDU, CDUMMY,
679: $ V, LDV, CDUMMY, -1, IERR )
680: LWRK_ZUNMQR2 = INT( CDUMMY(1) )
681: OPTWRK2 = MAX( LWRK_ZGEQP3, N/2+LWRK_ZGEQRF,
682: $ N/2+LWRK_ZGESVD2, N/2+LWRK_ZUNMQR2 )
683: IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON )
684: OPTWRK2 = N + OPTWRK2
685: OPTWRK = MAX( OPTWRK, OPTWRK2 )
686: END IF
687: ELSE
688: CALL ZGESVD( 'S', 'O', N, N, A, LDA, S, U, LDU,
689: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
690: LWRK_ZGESVD = INT( CDUMMY(1) )
691: OPTWRK = MAX(LWRK_ZGEQP3,LWRK_ZGESVD,LWRK_ZUNMQR)
692: IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON )
693: OPTWRK = N + OPTWRK
694: IF ( WNTVA ) THEN
695: CALL ZGELQF(N/2,N,U,LDU,CDUMMY,CDUMMY,-1,IERR)
696: LWRK_ZGELQF = INT( CDUMMY(1) )
697: CALL ZGESVD( 'S','O', N/2,N/2, V, LDV, S, U, LDU,
698: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
699: LWRK_ZGESVD2 = INT( CDUMMY(1) )
700: CALL ZUNMLQ( 'R', 'N', N, N, N/2, U, LDU, CDUMMY,
701: $ V, LDV, CDUMMY,-1,IERR )
702: LWRK_ZUNMLQ = INT( CDUMMY(1) )
703: OPTWRK2 = MAX( LWRK_ZGEQP3, N/2+LWRK_ZGELQF,
704: $ N/2+LWRK_ZGESVD2, N/2+LWRK_ZUNMLQ )
705: IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON )
706: OPTWRK2 = N + OPTWRK2
707: OPTWRK = MAX( OPTWRK, OPTWRK2 )
708: END IF
709: END IF
710: END IF
711: END IF
712: *
713: MINWRK = MAX( 2, MINWRK )
714: OPTWRK = MAX( 2, OPTWRK )
715: IF ( LCWORK .LT. MINWRK .AND. (.NOT.LQUERY) ) INFO = -19
716: *
717: END IF
718: *
719: IF (INFO .EQ. 0 .AND. LRWORK .LT. RMINWRK .AND. .NOT. LQUERY) THEN
720: INFO = -21
721: END IF
722: IF( INFO.NE.0 ) THEN
723: CALL XERBLA( 'ZGESVDQ', -INFO )
724: RETURN
725: ELSE IF ( LQUERY ) THEN
726: *
727: * Return optimal workspace
728: *
729: IWORK(1) = IMINWRK
730: CWORK(1) = OPTWRK
731: CWORK(2) = MINWRK
732: RWORK(1) = RMINWRK
733: RETURN
734: END IF
735: *
736: * Quick return if the matrix is void.
737: *
738: IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) THEN
739: * .. all output is void.
740: RETURN
741: END IF
742: *
743: BIG = DLAMCH('O')
744: ASCALED = .FALSE.
745: IF ( ROWPRM ) THEN
746: * .. reordering the rows in decreasing sequence in the
747: * ell-infinity norm - this enhances numerical robustness in
748: * the case of differently scaled rows.
749: DO 1904 p = 1, M
750: * RWORK(p) = ABS( A(p,IZAMAX(N,A(p,1),LDA)) )
751: * [[ZLANGE will return NaN if an entry of the p-th row is Nan]]
752: RWORK(p) = ZLANGE( 'M', 1, N, A(p,1), LDA, RDUMMY )
753: * .. check for NaN's and Inf's
754: IF ( ( RWORK(p) .NE. RWORK(p) ) .OR.
755: $ ( (RWORK(p)*ZERO) .NE. ZERO ) ) THEN
756: INFO = -8
757: CALL XERBLA( 'ZGESVDQ', -INFO )
758: RETURN
759: END IF
760: 1904 CONTINUE
761: DO 1952 p = 1, M - 1
762: q = IDAMAX( M-p+1, RWORK(p), 1 ) + p - 1
763: IWORK(N+p) = q
764: IF ( p .NE. q ) THEN
765: RTMP = RWORK(p)
766: RWORK(p) = RWORK(q)
767: RWORK(q) = RTMP
768: END IF
769: 1952 CONTINUE
770: *
771: IF ( RWORK(1) .EQ. ZERO ) THEN
772: * Quick return: A is the M x N zero matrix.
773: NUMRANK = 0
774: CALL DLASET( 'G', N, 1, ZERO, ZERO, S, N )
775: IF ( WNTUS ) CALL ZLASET('G', M, N, CZERO, CONE, U, LDU)
776: IF ( WNTUA ) CALL ZLASET('G', M, M, CZERO, CONE, U, LDU)
777: IF ( WNTVA ) CALL ZLASET('G', N, N, CZERO, CONE, V, LDV)
778: IF ( WNTUF ) THEN
779: CALL ZLASET( 'G', N, 1, CZERO, CZERO, CWORK, N )
780: CALL ZLASET( 'G', M, N, CZERO, CONE, U, LDU )
781: END IF
782: DO 5001 p = 1, N
783: IWORK(p) = p
784: 5001 CONTINUE
785: IF ( ROWPRM ) THEN
786: DO 5002 p = N + 1, N + M - 1
787: IWORK(p) = p - N
788: 5002 CONTINUE
789: END IF
790: IF ( CONDA ) RWORK(1) = -1
791: RWORK(2) = -1
792: RETURN
793: END IF
794: *
795: IF ( RWORK(1) .GT. BIG / SQRT(DBLE(M)) ) THEN
796: * .. to prevent overflow in the QR factorization, scale the
797: * matrix by 1/sqrt(M) if too large entry detected
798: CALL ZLASCL('G',0,0,SQRT(DBLE(M)),ONE, M,N, A,LDA, IERR)
799: ASCALED = .TRUE.
800: END IF
801: CALL ZLASWP( N, A, LDA, 1, M-1, IWORK(N+1), 1 )
802: END IF
803: *
804: * .. At this stage, preemptive scaling is done only to avoid column
805: * norms overflows during the QR factorization. The SVD procedure should
806: * have its own scaling to save the singular values from overflows and
807: * underflows. That depends on the SVD procedure.
808: *
809: IF ( .NOT.ROWPRM ) THEN
810: RTMP = ZLANGE( 'M', M, N, A, LDA, RWORK )
811: IF ( ( RTMP .NE. RTMP ) .OR.
812: $ ( (RTMP*ZERO) .NE. ZERO ) ) THEN
813: INFO = -8
814: CALL XERBLA( 'ZGESVDQ', -INFO )
815: RETURN
816: END IF
817: IF ( RTMP .GT. BIG / SQRT(DBLE(M)) ) THEN
818: * .. to prevent overflow in the QR factorization, scale the
819: * matrix by 1/sqrt(M) if too large entry detected
820: CALL ZLASCL('G',0,0, SQRT(DBLE(M)),ONE, M,N, A,LDA, IERR)
821: ASCALED = .TRUE.
822: END IF
823: END IF
824: *
825: * .. QR factorization with column pivoting
826: *
827: * A * P = Q * [ R ]
828: * [ 0 ]
829: *
830: DO 1963 p = 1, N
831: * .. all columns are free columns
832: IWORK(p) = 0
833: 1963 CONTINUE
834: CALL ZGEQP3( M, N, A, LDA, IWORK, CWORK, CWORK(N+1), LCWORK-N,
835: $ RWORK, IERR )
836: *
837: * If the user requested accuracy level allows truncation in the
838: * computed upper triangular factor, the matrix R is examined and,
839: * if possible, replaced with its leading upper trapezoidal part.
840: *
841: EPSLN = DLAMCH('E')
842: SFMIN = DLAMCH('S')
843: * SMALL = SFMIN / EPSLN
844: NR = N
845: *
846: IF ( ACCLA ) THEN
847: *
848: * Standard absolute error bound suffices. All sigma_i with
849: * sigma_i < N*EPS*||A||_F are flushed to zero. This is an
850: * aggressive enforcement of lower numerical rank by introducing a
851: * backward error of the order of N*EPS*||A||_F.
852: NR = 1
853: RTMP = SQRT(DBLE(N))*EPSLN
854: DO 3001 p = 2, N
855: IF ( ABS(A(p,p)) .LT. (RTMP*ABS(A(1,1))) ) GO TO 3002
856: NR = NR + 1
857: 3001 CONTINUE
858: 3002 CONTINUE
859: *
860: ELSEIF ( ACCLM ) THEN
861: * .. similarly as above, only slightly more gentle (less aggressive).
862: * Sudden drop on the diagonal of R is used as the criterion for being
863: * close-to-rank-deficient. The threshold is set to EPSLN=DLAMCH('E').
864: * [[This can be made more flexible by replacing this hard-coded value
865: * with a user specified threshold.]] Also, the values that underflow
866: * will be truncated.
867: NR = 1
868: DO 3401 p = 2, N
869: IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR.
870: $ ( ABS(A(p,p)) .LT. SFMIN ) ) GO TO 3402
871: NR = NR + 1
872: 3401 CONTINUE
873: 3402 CONTINUE
874: *
875: ELSE
876: * .. RRQR not authorized to determine numerical rank except in the
877: * obvious case of zero pivots.
878: * .. inspect R for exact zeros on the diagonal;
879: * R(i,i)=0 => R(i:N,i:N)=0.
880: NR = 1
881: DO 3501 p = 2, N
882: IF ( ABS(A(p,p)) .EQ. ZERO ) GO TO 3502
883: NR = NR + 1
884: 3501 CONTINUE
885: 3502 CONTINUE
886: *
887: IF ( CONDA ) THEN
888: * Estimate the scaled condition number of A. Use the fact that it is
889: * the same as the scaled condition number of R.
890: * .. V is used as workspace
891: CALL ZLACPY( 'U', N, N, A, LDA, V, LDV )
892: * Only the leading NR x NR submatrix of the triangular factor
893: * is considered. Only if NR=N will this give a reliable error
894: * bound. However, even for NR < N, this can be used on an
895: * expert level and obtain useful information in the sense of
896: * perturbation theory.
897: DO 3053 p = 1, NR
898: RTMP = DZNRM2( p, V(1,p), 1 )
899: CALL ZDSCAL( p, ONE/RTMP, V(1,p), 1 )
900: 3053 CONTINUE
901: IF ( .NOT. ( LSVEC .OR. RSVEC ) ) THEN
902: CALL ZPOCON( 'U', NR, V, LDV, ONE, RTMP,
903: $ CWORK, RWORK, IERR )
904: ELSE
905: CALL ZPOCON( 'U', NR, V, LDV, ONE, RTMP,
906: $ CWORK(N+1), RWORK, IERR )
907: END IF
908: SCONDA = ONE / SQRT(RTMP)
909: * For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1),
910: * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
911: * See the reference [1] for more details.
912: END IF
913: *
914: ENDIF
915: *
916: IF ( WNTUR ) THEN
917: N1 = NR
918: ELSE IF ( WNTUS .OR. WNTUF) THEN
919: N1 = N
920: ELSE IF ( WNTUA ) THEN
921: N1 = M
922: END IF
923: *
924: IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
925: *.......................................................................
926: * .. only the singular values are requested
927: *.......................................................................
928: IF ( RTRANS ) THEN
929: *
930: * .. compute the singular values of R**H = [A](1:NR,1:N)**H
931: * .. set the lower triangle of [A] to [A](1:NR,1:N)**H and
932: * the upper triangle of [A] to zero.
933: DO 1146 p = 1, MIN( N, NR )
934: A(p,p) = CONJG(A(p,p))
935: DO 1147 q = p + 1, N
936: A(q,p) = CONJG(A(p,q))
937: IF ( q .LE. NR ) A(p,q) = CZERO
938: 1147 CONTINUE
939: 1146 CONTINUE
940: *
941: CALL ZGESVD( 'N', 'N', N, NR, A, LDA, S, U, LDU,
942: $ V, LDV, CWORK, LCWORK, RWORK, INFO )
943: *
944: ELSE
945: *
946: * .. compute the singular values of R = [A](1:NR,1:N)
947: *
948: IF ( NR .GT. 1 )
949: $ CALL ZLASET( 'L', NR-1,NR-1, CZERO,CZERO, A(2,1), LDA )
950: CALL ZGESVD( 'N', 'N', NR, N, A, LDA, S, U, LDU,
951: $ V, LDV, CWORK, LCWORK, RWORK, INFO )
952: *
953: END IF
954: *
955: ELSE IF ( LSVEC .AND. ( .NOT. RSVEC) ) THEN
956: *.......................................................................
957: * .. the singular values and the left singular vectors requested
958: *.......................................................................""""""""
959: IF ( RTRANS ) THEN
960: * .. apply ZGESVD to R**H
961: * .. copy R**H into [U] and overwrite [U] with the right singular
962: * vectors of R
963: DO 1192 p = 1, NR
964: DO 1193 q = p, N
965: U(q,p) = CONJG(A(p,q))
966: 1193 CONTINUE
967: 1192 CONTINUE
968: IF ( NR .GT. 1 )
969: $ CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, U(1,2), LDU )
970: * .. the left singular vectors not computed, the NR right singular
971: * vectors overwrite [U](1:NR,1:NR) as conjugate transposed. These
972: * will be pre-multiplied by Q to build the left singular vectors of A.
973: CALL ZGESVD( 'N', 'O', N, NR, U, LDU, S, U, LDU,
974: $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO )
975: *
976: DO 1119 p = 1, NR
977: U(p,p) = CONJG(U(p,p))
978: DO 1120 q = p + 1, NR
979: CTMP = CONJG(U(q,p))
980: U(q,p) = CONJG(U(p,q))
981: U(p,q) = CTMP
982: 1120 CONTINUE
983: 1119 CONTINUE
984: *
985: ELSE
986: * .. apply ZGESVD to R
987: * .. copy R into [U] and overwrite [U] with the left singular vectors
988: CALL ZLACPY( 'U', NR, N, A, LDA, U, LDU )
989: IF ( NR .GT. 1 )
990: $ CALL ZLASET( 'L', NR-1, NR-1, CZERO, CZERO, U(2,1), LDU )
991: * .. the right singular vectors not computed, the NR left singular
992: * vectors overwrite [U](1:NR,1:NR)
993: CALL ZGESVD( 'O', 'N', NR, N, U, LDU, S, U, LDU,
994: $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO )
995: * .. now [U](1:NR,1:NR) contains the NR left singular vectors of
996: * R. These will be pre-multiplied by Q to build the left singular
997: * vectors of A.
998: END IF
999: *
1000: * .. assemble the left singular vector matrix U of dimensions
1001: * (M x NR) or (M x N) or (M x M).
1002: IF ( ( NR .LT. M ) .AND. ( .NOT.WNTUF ) ) THEN
1003: CALL ZLASET('A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU)
1004: IF ( NR .LT. N1 ) THEN
1005: CALL ZLASET( 'A',NR,N1-NR,CZERO,CZERO,U(1,NR+1), LDU )
1006: CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,
1007: $ U(NR+1,NR+1), LDU )
1008: END IF
1009: END IF
1010: *
1011: * The Q matrix from the first QRF is built into the left singular
1012: * vectors matrix U.
1013: *
1014: IF ( .NOT.WNTUF )
1015: $ CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
1016: $ LDU, CWORK(N+1), LCWORK-N, IERR )
1017: IF ( ROWPRM .AND. .NOT.WNTUF )
1018: $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 )
1019: *
1020: ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN
1021: *.......................................................................
1022: * .. the singular values and the right singular vectors requested
1023: *.......................................................................
1024: IF ( RTRANS ) THEN
1025: * .. apply ZGESVD to R**H
1026: * .. copy R**H into V and overwrite V with the left singular vectors
1027: DO 1165 p = 1, NR
1028: DO 1166 q = p, N
1029: V(q,p) = CONJG(A(p,q))
1030: 1166 CONTINUE
1031: 1165 CONTINUE
1032: IF ( NR .GT. 1 )
1033: $ CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV )
1034: * .. the left singular vectors of R**H overwrite V, the right singular
1035: * vectors not computed
1036: IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
1037: CALL ZGESVD( 'O', 'N', N, NR, V, LDV, S, U, LDU,
1038: $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO )
1039: *
1040: DO 1121 p = 1, NR
1041: V(p,p) = CONJG(V(p,p))
1042: DO 1122 q = p + 1, NR
1043: CTMP = CONJG(V(q,p))
1044: V(q,p) = CONJG(V(p,q))
1045: V(p,q) = CTMP
1046: 1122 CONTINUE
1047: 1121 CONTINUE
1048: *
1049: IF ( NR .LT. N ) THEN
1050: DO 1103 p = 1, NR
1051: DO 1104 q = NR + 1, N
1052: V(p,q) = CONJG(V(q,p))
1053: 1104 CONTINUE
1054: 1103 CONTINUE
1055: END IF
1056: CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK )
1057: ELSE
1058: * .. need all N right singular vectors and NR < N
1059: * [!] This is simple implementation that augments [V](1:N,1:NR)
1060: * by padding a zero block. In the case NR << N, a more efficient
1061: * way is to first use the QR factorization. For more details
1062: * how to implement this, see the " FULL SVD " branch.
1063: CALL ZLASET('G', N, N-NR, CZERO, CZERO, V(1,NR+1), LDV)
1064: CALL ZGESVD( 'O', 'N', N, N, V, LDV, S, U, LDU,
1065: $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO )
1066: *
1067: DO 1123 p = 1, N
1068: V(p,p) = CONJG(V(p,p))
1069: DO 1124 q = p + 1, N
1070: CTMP = CONJG(V(q,p))
1071: V(q,p) = CONJG(V(p,q))
1072: V(p,q) = CTMP
1073: 1124 CONTINUE
1074: 1123 CONTINUE
1075: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
1076: END IF
1077: *
1078: ELSE
1079: * .. aply ZGESVD to R
1080: * .. copy R into V and overwrite V with the right singular vectors
1081: CALL ZLACPY( 'U', NR, N, A, LDA, V, LDV )
1082: IF ( NR .GT. 1 )
1083: $ CALL ZLASET( 'L', NR-1, NR-1, CZERO, CZERO, V(2,1), LDV )
1084: * .. the right singular vectors overwrite V, the NR left singular
1085: * vectors stored in U(1:NR,1:NR)
1086: IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
1087: CALL ZGESVD( 'N', 'O', NR, N, V, LDV, S, U, LDU,
1088: $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO )
1089: CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK )
1090: * .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H
1091: ELSE
1092: * .. need all N right singular vectors and NR < N
1093: * [!] This is simple implementation that augments [V](1:NR,1:N)
1094: * by padding a zero block. In the case NR << N, a more efficient
1095: * way is to first use the LQ factorization. For more details
1096: * how to implement this, see the " FULL SVD " branch.
1097: CALL ZLASET('G', N-NR, N, CZERO,CZERO, V(NR+1,1), LDV)
1098: CALL ZGESVD( 'N', 'O', N, N, V, LDV, S, U, LDU,
1099: $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO )
1100: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
1101: END IF
1102: * .. now [V] contains the adjoint of the matrix of the right singular
1103: * vectors of A.
1104: END IF
1105: *
1106: ELSE
1107: *.......................................................................
1108: * .. FULL SVD requested
1109: *.......................................................................
1110: IF ( RTRANS ) THEN
1111: *
1112: * .. apply ZGESVD to R**H [[this option is left for R&D&T]]
1113: *
1114: IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
1115: * .. copy R**H into [V] and overwrite [V] with the left singular
1116: * vectors of R**H
1117: DO 1168 p = 1, NR
1118: DO 1169 q = p, N
1119: V(q,p) = CONJG(A(p,q))
1120: 1169 CONTINUE
1121: 1168 CONTINUE
1122: IF ( NR .GT. 1 )
1123: $ CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV )
1124: *
1125: * .. the left singular vectors of R**H overwrite [V], the NR right
1126: * singular vectors of R**H stored in [U](1:NR,1:NR) as conjugate
1127: * transposed
1128: CALL ZGESVD( 'O', 'A', N, NR, V, LDV, S, V, LDV,
1129: $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO )
1130: * .. assemble V
1131: DO 1115 p = 1, NR
1132: V(p,p) = CONJG(V(p,p))
1133: DO 1116 q = p + 1, NR
1134: CTMP = CONJG(V(q,p))
1135: V(q,p) = CONJG(V(p,q))
1136: V(p,q) = CTMP
1137: 1116 CONTINUE
1138: 1115 CONTINUE
1139: IF ( NR .LT. N ) THEN
1140: DO 1101 p = 1, NR
1141: DO 1102 q = NR+1, N
1142: V(p,q) = CONJG(V(q,p))
1143: 1102 CONTINUE
1144: 1101 CONTINUE
1145: END IF
1146: CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK )
1147: *
1148: DO 1117 p = 1, NR
1149: U(p,p) = CONJG(U(p,p))
1150: DO 1118 q = p + 1, NR
1151: CTMP = CONJG(U(q,p))
1152: U(q,p) = CONJG(U(p,q))
1153: U(p,q) = CTMP
1154: 1118 CONTINUE
1155: 1117 CONTINUE
1156: *
1157: IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
1158: CALL ZLASET('A', M-NR,NR, CZERO,CZERO, U(NR+1,1), LDU)
1159: IF ( NR .LT. N1 ) THEN
1160: CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
1161: CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,
1162: $ U(NR+1,NR+1), LDU )
1163: END IF
1164: END IF
1165: *
1166: ELSE
1167: * .. need all N right singular vectors and NR < N
1168: * .. copy R**H into [V] and overwrite [V] with the left singular
1169: * vectors of R**H
1170: * [[The optimal ratio N/NR for using QRF instead of padding
1171: * with zeros. Here hard coded to 2; it must be at least
1172: * two due to work space constraints.]]
1173: * OPTRATIO = ILAENV(6, 'ZGESVD', 'S' // 'O', NR,N,0,0)
1174: * OPTRATIO = MAX( OPTRATIO, 2 )
1175: OPTRATIO = 2
1176: IF ( OPTRATIO*NR .GT. N ) THEN
1177: DO 1198 p = 1, NR
1178: DO 1199 q = p, N
1179: V(q,p) = CONJG(A(p,q))
1180: 1199 CONTINUE
1181: 1198 CONTINUE
1182: IF ( NR .GT. 1 )
1183: $ CALL ZLASET('U',NR-1,NR-1, CZERO,CZERO, V(1,2),LDV)
1184: *
1185: CALL ZLASET('A',N,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
1186: CALL ZGESVD( 'O', 'A', N, N, V, LDV, S, V, LDV,
1187: $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO )
1188: *
1189: DO 1113 p = 1, N
1190: V(p,p) = CONJG(V(p,p))
1191: DO 1114 q = p + 1, N
1192: CTMP = CONJG(V(q,p))
1193: V(q,p) = CONJG(V(p,q))
1194: V(p,q) = CTMP
1195: 1114 CONTINUE
1196: 1113 CONTINUE
1197: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
1198: * .. assemble the left singular vector matrix U of dimensions
1199: * (M x N1), i.e. (M x N) or (M x M).
1200: *
1201: DO 1111 p = 1, N
1202: U(p,p) = CONJG(U(p,p))
1203: DO 1112 q = p + 1, N
1204: CTMP = CONJG(U(q,p))
1205: U(q,p) = CONJG(U(p,q))
1206: U(p,q) = CTMP
1207: 1112 CONTINUE
1208: 1111 CONTINUE
1209: *
1210: IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN
1211: CALL ZLASET('A',M-N,N,CZERO,CZERO,U(N+1,1),LDU)
1212: IF ( N .LT. N1 ) THEN
1213: CALL ZLASET('A',N,N1-N,CZERO,CZERO,U(1,N+1),LDU)
1214: CALL ZLASET('A',M-N,N1-N,CZERO,CONE,
1215: $ U(N+1,N+1), LDU )
1216: END IF
1217: END IF
1218: ELSE
1219: * .. copy R**H into [U] and overwrite [U] with the right
1220: * singular vectors of R
1221: DO 1196 p = 1, NR
1222: DO 1197 q = p, N
1223: U(q,NR+p) = CONJG(A(p,q))
1224: 1197 CONTINUE
1225: 1196 CONTINUE
1226: IF ( NR .GT. 1 )
1227: $ CALL ZLASET('U',NR-1,NR-1,CZERO,CZERO,U(1,NR+2),LDU)
1228: CALL ZGEQRF( N, NR, U(1,NR+1), LDU, CWORK(N+1),
1229: $ CWORK(N+NR+1), LCWORK-N-NR, IERR )
1230: DO 1143 p = 1, NR
1231: DO 1144 q = 1, N
1232: V(q,p) = CONJG(U(p,NR+q))
1233: 1144 CONTINUE
1234: 1143 CONTINUE
1235: CALL ZLASET('U',NR-1,NR-1,CZERO,CZERO,V(1,2),LDV)
1236: CALL ZGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU,
1237: $ V,LDV, CWORK(N+NR+1),LCWORK-N-NR,RWORK, INFO )
1238: CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV)
1239: CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
1240: CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
1241: CALL ZUNMQR('R','C', N, N, NR, U(1,NR+1), LDU,
1242: $ CWORK(N+1),V,LDV,CWORK(N+NR+1),LCWORK-N-NR,IERR)
1243: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
1244: * .. assemble the left singular vector matrix U of dimensions
1245: * (M x NR) or (M x N) or (M x M).
1246: IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
1247: CALL ZLASET('A',M-NR,NR,CZERO,CZERO,U(NR+1,1),LDU)
1248: IF ( NR .LT. N1 ) THEN
1249: CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
1250: CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,
1251: $ U(NR+1,NR+1),LDU)
1252: END IF
1253: END IF
1254: END IF
1255: END IF
1256: *
1257: ELSE
1258: *
1259: * .. apply ZGESVD to R [[this is the recommended option]]
1260: *
1261: IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
1262: * .. copy R into [V] and overwrite V with the right singular vectors
1263: CALL ZLACPY( 'U', NR, N, A, LDA, V, LDV )
1264: IF ( NR .GT. 1 )
1265: $ CALL ZLASET( 'L', NR-1,NR-1, CZERO,CZERO, V(2,1), LDV )
1266: * .. the right singular vectors of R overwrite [V], the NR left
1267: * singular vectors of R stored in [U](1:NR,1:NR)
1268: CALL ZGESVD( 'S', 'O', NR, N, V, LDV, S, U, LDU,
1269: $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO )
1270: CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK )
1271: * .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H
1272: * .. assemble the left singular vector matrix U of dimensions
1273: * (M x NR) or (M x N) or (M x M).
1274: IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
1275: CALL ZLASET('A', M-NR,NR, CZERO,CZERO, U(NR+1,1), LDU)
1276: IF ( NR .LT. N1 ) THEN
1277: CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
1278: CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,
1279: $ U(NR+1,NR+1), LDU )
1280: END IF
1281: END IF
1282: *
1283: ELSE
1284: * .. need all N right singular vectors and NR < N
1285: * .. the requested number of the left singular vectors
1286: * is then N1 (N or M)
1287: * [[The optimal ratio N/NR for using LQ instead of padding
1288: * with zeros. Here hard coded to 2; it must be at least
1289: * two due to work space constraints.]]
1290: * OPTRATIO = ILAENV(6, 'ZGESVD', 'S' // 'O', NR,N,0,0)
1291: * OPTRATIO = MAX( OPTRATIO, 2 )
1292: OPTRATIO = 2
1293: IF ( OPTRATIO * NR .GT. N ) THEN
1294: CALL ZLACPY( 'U', NR, N, A, LDA, V, LDV )
1295: IF ( NR .GT. 1 )
1296: $ CALL ZLASET('L', NR-1,NR-1, CZERO,CZERO, V(2,1),LDV)
1297: * .. the right singular vectors of R overwrite [V], the NR left
1298: * singular vectors of R stored in [U](1:NR,1:NR)
1299: CALL ZLASET('A', N-NR,N, CZERO,CZERO, V(NR+1,1),LDV)
1300: CALL ZGESVD( 'S', 'O', N, N, V, LDV, S, U, LDU,
1301: $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO )
1302: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
1303: * .. now [V] contains the adjoint of the matrix of the right
1304: * singular vectors of A. The leading N left singular vectors
1305: * are in [U](1:N,1:N)
1306: * .. assemble the left singular vector matrix U of dimensions
1307: * (M x N1), i.e. (M x N) or (M x M).
1308: IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN
1309: CALL ZLASET('A',M-N,N,CZERO,CZERO,U(N+1,1),LDU)
1310: IF ( N .LT. N1 ) THEN
1311: CALL ZLASET('A',N,N1-N,CZERO,CZERO,U(1,N+1),LDU)
1312: CALL ZLASET( 'A',M-N,N1-N,CZERO,CONE,
1313: $ U(N+1,N+1), LDU )
1314: END IF
1315: END IF
1316: ELSE
1317: CALL ZLACPY( 'U', NR, N, A, LDA, U(NR+1,1), LDU )
1318: IF ( NR .GT. 1 )
1319: $ CALL ZLASET('L',NR-1,NR-1,CZERO,CZERO,U(NR+2,1),LDU)
1320: CALL ZGELQF( NR, N, U(NR+1,1), LDU, CWORK(N+1),
1321: $ CWORK(N+NR+1), LCWORK-N-NR, IERR )
1322: CALL ZLACPY('L',NR,NR,U(NR+1,1),LDU,V,LDV)
1323: IF ( NR .GT. 1 )
1324: $ CALL ZLASET('U',NR-1,NR-1,CZERO,CZERO,V(1,2),LDV)
1325: CALL ZGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU,
1326: $ V, LDV, CWORK(N+NR+1), LCWORK-N-NR, RWORK, INFO )
1327: CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV)
1328: CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
1329: CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
1330: CALL ZUNMLQ('R','N',N,N,NR,U(NR+1,1),LDU,CWORK(N+1),
1331: $ V, LDV, CWORK(N+NR+1),LCWORK-N-NR,IERR)
1332: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
1333: * .. assemble the left singular vector matrix U of dimensions
1334: * (M x NR) or (M x N) or (M x M).
1335: IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
1336: CALL ZLASET('A',M-NR,NR,CZERO,CZERO,U(NR+1,1),LDU)
1337: IF ( NR .LT. N1 ) THEN
1338: CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
1339: CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,
1340: $ U(NR+1,NR+1), LDU )
1341: END IF
1342: END IF
1343: END IF
1344: END IF
1345: * .. end of the "R**H or R" branch
1346: END IF
1347: *
1348: * The Q matrix from the first QRF is built into the left singular
1349: * vectors matrix U.
1350: *
1351: IF ( .NOT. WNTUF )
1352: $ CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
1353: $ LDU, CWORK(N+1), LCWORK-N, IERR )
1354: IF ( ROWPRM .AND. .NOT.WNTUF )
1355: $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 )
1356: *
1357: * ... end of the "full SVD" branch
1358: END IF
1359: *
1360: * Check whether some singular values are returned as zeros, e.g.
1361: * due to underflow, and update the numerical rank.
1362: p = NR
1363: DO 4001 q = p, 1, -1
1364: IF ( S(q) .GT. ZERO ) GO TO 4002
1365: NR = NR - 1
1366: 4001 CONTINUE
1367: 4002 CONTINUE
1368: *
1369: * .. if numerical rank deficiency is detected, the truncated
1370: * singular values are set to zero.
1371: IF ( NR .LT. N ) CALL DLASET( 'G', N-NR,1, ZERO,ZERO, S(NR+1), N )
1372: * .. undo scaling; this may cause overflow in the largest singular
1373: * values.
1374: IF ( ASCALED )
1375: $ CALL DLASCL( 'G',0,0, ONE,SQRT(DBLE(M)), NR,1, S, N, IERR )
1376: IF ( CONDA ) RWORK(1) = SCONDA
1377: RWORK(2) = p - NR
1378: * .. p-NR is the number of singular values that are computed as
1379: * exact zeros in ZGESVD() applied to the (possibly truncated)
1380: * full row rank triangular (trapezoidal) factor of A.
1381: NUMRANK = NR
1382: *
1383: RETURN
1384: *
1385: * End of ZGESVDQ
1386: *
1387: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgesvdq.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack