LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dgesvdq.f
Go to the documentation of this file.
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, enclosed 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 gesvdq
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,
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
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgelqf(m, n, a, lda, tau, work, lwork, info)
DGELQF
Definition dgelqf.f:143
subroutine dgeqp3(m, n, a, lda, jpvt, tau, work, lwork, info)
DGEQP3
Definition dgeqp3.f:151
subroutine dgeqrf(m, n, a, lda, tau, work, lwork, info)
DGEQRF
Definition dgeqrf.f:146
subroutine dgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
DGESVD computes the singular value decomposition (SVD) for GE matrices
Definition dgesvd.f:211
subroutine dgesvdq(joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, work, lwork, rwork, lrwork, info)
DGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE...
Definition dgesvdq.f:415
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
subroutine dlapmt(forwrd, m, n, x, ldx, k)
DLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition dlapmt.f:104
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine dlaswp(n, a, lda, k1, k2, ipiv, incx)
DLASWP performs a series of row interchanges on a general rectangular matrix.
Definition dlaswp.f:115
subroutine dpocon(uplo, n, a, lda, anorm, rcond, work, iwork, info)
DPOCON
Definition dpocon.f:121
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dormlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMLQ
Definition dormlq.f:167
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:167