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