LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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, 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,
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 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 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 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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
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 dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
Definition: dgelqf.f:143
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 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 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 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 dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:167
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
Definition: dormlq.f:167
subroutine dpocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
DPOCON
Definition: dpocon.f:121