LAPACK  3.9.0
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 *> LIWORK(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 trapeziodal
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 *> \date November 2018
410 *
411 *> \ingroup doubleGEsing
412 *
413 * =====================================================================
414  SUBROUTINE dgesvdq( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
415  $ S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
416  $ WORK, LWORK, RWORK, LRWORK, INFO )
417 * .. Scalar Arguments ..
418  IMPLICIT NONE
419  CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV
420  INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LWORK, LRWORK,
421  $ info
422 * ..
423 * .. Array Arguments ..
424  DOUBLE PRECISION A( LDA, * ), U( LDU, * ), V( LDV, * ), WORK( * )
425  DOUBLE PRECISION S( * ), RWORK( * )
426  INTEGER IWORK( * )
427 *
428 * =====================================================================
429 *
430 * .. Parameters ..
431  DOUBLE PRECISION ZERO, ONE
432  PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
433 * .. Local Scalars ..
434  INTEGER IERR, IWOFF, NR, N1, OPTRATIO, p, q
435  INTEGER LWCON, LWQP3, LWRK_DGELQF, LWRK_DGESVD, LWRK_DGESVD2,
436  $ lwrk_dgeqp3, lwrk_dgeqrf, lwrk_dormlq, lwrk_dormqr,
437  $ lwrk_dormqr2, lwlqf, lwqrf, lwsvd, lwsvd2, lworq,
438  $ lworq2, lworlq, minwrk, minwrk2, optwrk, optwrk2,
439  $ iminwrk, rminwrk
440  LOGICAL ACCLA, ACCLM, ACCLH, ASCALED, CONDA, DNTWU, DNTWV,
441  $ LQUERY, LSVC0, LSVEC, ROWPRM, RSVEC, RTRANS, WNTUA,
442  $ wntuf, wntur, wntus, wntva, wntvr
443  DOUBLE PRECISION BIG, EPSLN, RTMP, SCONDA, SFMIN
444 * .. Local Arrays
445  DOUBLE PRECISION RDUMMY(1)
446 * ..
447 * .. External Subroutines (BLAS, LAPACK)
448  EXTERNAL dgelqf, dgeqp3, dgeqrf, dgesvd, dlacpy, dlapmt,
449  $ dlascl, dlaset, dlaswp, dscal, dpocon, dormlq,
450  $ dormqr, xerbla
451 * ..
452 * .. External Functions (BLAS, LAPACK)
453  LOGICAL LSAME
454  INTEGER IDAMAX
455  DOUBLE PRECISION DLANGE, DNRM2, DLAMCH
456  EXTERNAL dlange, lsame, idamax, dnrm2, dlamch
457 * ..
458 * .. Intrinsic Functions ..
459 *
460  INTRINSIC abs, max, min, dble, sqrt
461 *
462 * Test the input arguments
463 *
464  wntus = lsame( jobu, 'S' ) .OR. lsame( jobu, 'U' )
465  wntur = lsame( jobu, 'R' )
466  wntua = lsame( jobu, 'A' )
467  wntuf = lsame( jobu, 'F' )
468  lsvc0 = wntus .OR. wntur .OR. wntua
469  lsvec = lsvc0 .OR. wntuf
470  dntwu = lsame( jobu, 'N' )
471 *
472  wntvr = lsame( jobv, 'R' )
473  wntva = lsame( jobv, 'A' ) .OR. lsame( jobv, 'V' )
474  rsvec = wntvr .OR. wntva
475  dntwv = lsame( jobv, 'N' )
476 *
477  accla = lsame( joba, 'A' )
478  acclm = lsame( joba, 'M' )
479  conda = lsame( joba, 'E' )
480  acclh = lsame( joba, 'H' ) .OR. conda
481 *
482  rowprm = lsame( jobp, 'P' )
483  rtrans = lsame( jobr, 'T' )
484 *
485  IF ( rowprm ) THEN
486  IF ( conda ) THEN
487  iminwrk = max( 1, n + m - 1 + n )
488  ELSE
489  iminwrk = max( 1, n + m - 1 )
490  END IF
491  rminwrk = max( 2, m )
492  ELSE
493  IF ( conda ) THEN
494  iminwrk = max( 1, n + n )
495  ELSE
496  iminwrk = max( 1, n )
497  END IF
498  rminwrk = 2
499  END IF
500  lquery = (liwork .EQ. -1 .OR. lwork .EQ. -1 .OR. lrwork .EQ. -1)
501  info = 0
502  IF ( .NOT. ( accla .OR. acclm .OR. acclh ) ) THEN
503  info = -1
504  ELSE IF ( .NOT.( rowprm .OR. lsame( jobp, 'N' ) ) ) THEN
505  info = -2
506  ELSE IF ( .NOT.( rtrans .OR. lsame( jobr, 'N' ) ) ) THEN
507  info = -3
508  ELSE IF ( .NOT.( lsvec .OR. dntwu ) ) THEN
509  info = -4
510  ELSE IF ( wntur .AND. wntva ) THEN
511  info = -5
512  ELSE IF ( .NOT.( rsvec .OR. dntwv )) THEN
513  info = -5
514  ELSE IF ( m.LT.0 ) THEN
515  info = -6
516  ELSE IF ( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
517  info = -7
518  ELSE IF ( lda.LT.max( 1, m ) ) THEN
519  info = -9
520  ELSE IF ( ldu.LT.1 .OR. ( lsvc0 .AND. ldu.LT.m ) .OR.
521  $ ( wntuf .AND. ldu.LT.n ) ) THEN
522  info = -12
523  ELSE IF ( ldv.LT.1 .OR. ( rsvec .AND. ldv.LT.n ) .OR.
524  $ ( conda .AND. ldv.LT.n ) ) THEN
525  info = -14
526  ELSE IF ( liwork .LT. iminwrk .AND. .NOT. lquery ) THEN
527  info = -17
528  END IF
529 *
530 *
531  IF ( info .EQ. 0 ) THEN
532 * .. compute the minimal and the optimal workspace lengths
533 * [[The expressions for computing the minimal and the optimal
534 * values of LWORK are written with a lot of redundancy and
535 * can be simplified. However, this detailed form is easier for
536 * maintenance and modifications of the code.]]
537 *
538 * .. minimal workspace length for DGEQP3 of an M x N matrix
539  lwqp3 = 3 * n + 1
540 * .. minimal workspace length for DORMQR to build left singular vectors
541  IF ( wntus .OR. wntur ) THEN
542  lworq = max( n , 1 )
543  ELSE IF ( wntua ) THEN
544  lworq = max( m , 1 )
545  END IF
546 * .. minimal workspace length for DPOCON of an N x N matrix
547  lwcon = 3 * n
548 * .. DGESVD of an N x N matrix
549  lwsvd = max( 5 * n, 1 )
550  IF ( lquery ) THEN
551  CALL dgeqp3( m, n, a, lda, iwork, rdummy, rdummy, -1,
552  $ ierr )
553  lwrk_dgeqp3 = int( rdummy(1) )
554  IF ( wntus .OR. wntur ) THEN
555  CALL dormqr( 'L', 'N', m, n, n, a, lda, rdummy, u,
556  $ ldu, rdummy, -1, ierr )
557  lwrk_dormqr = int( rdummy(1) )
558  ELSE IF ( wntua ) THEN
559  CALL dormqr( 'L', 'N', m, m, n, a, lda, rdummy, u,
560  $ ldu, rdummy, -1, ierr )
561  lwrk_dormqr = int( rdummy(1) )
562  ELSE
563  lwrk_dormqr = 0
564  END IF
565  END IF
566  minwrk = 2
567  optwrk = 2
568  IF ( .NOT. (lsvec .OR. rsvec )) THEN
569 * .. minimal and optimal sizes of the workspace if
570 * only the singular values are requested
571  IF ( conda ) THEN
572  minwrk = max( n+lwqp3, lwcon, lwsvd )
573  ELSE
574  minwrk = max( n+lwqp3, lwsvd )
575  END IF
576  IF ( lquery ) THEN
577  CALL dgesvd( 'N', 'N', n, n, a, lda, s, u, ldu,
578  $ v, ldv, rdummy, -1, ierr )
579  lwrk_dgesvd = int( rdummy(1) )
580  IF ( conda ) THEN
581  optwrk = max( n+lwrk_dgeqp3, n+lwcon, lwrk_dgesvd )
582  ELSE
583  optwrk = max( n+lwrk_dgeqp3, lwrk_dgesvd )
584  END IF
585  END IF
586  ELSE IF ( lsvec .AND. (.NOT.rsvec) ) THEN
587 * .. minimal and optimal sizes of the workspace if the
588 * singular values and the left singular vectors are requested
589  IF ( conda ) THEN
590  minwrk = n + max( lwqp3, lwcon, lwsvd, lworq )
591  ELSE
592  minwrk = n + max( lwqp3, lwsvd, lworq )
593  END IF
594  IF ( lquery ) THEN
595  IF ( rtrans ) THEN
596  CALL dgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
597  $ v, ldv, rdummy, -1, ierr )
598  ELSE
599  CALL dgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
600  $ v, ldv, rdummy, -1, ierr )
601  END IF
602  lwrk_dgesvd = int( rdummy(1) )
603  IF ( conda ) THEN
604  optwrk = n + max( lwrk_dgeqp3, lwcon, lwrk_dgesvd,
605  $ lwrk_dormqr )
606  ELSE
607  optwrk = n + max( lwrk_dgeqp3, lwrk_dgesvd,
608  $ lwrk_dormqr )
609  END IF
610  END IF
611  ELSE IF ( rsvec .AND. (.NOT.lsvec) ) THEN
612 * .. minimal and optimal sizes of the workspace if the
613 * singular values and the right singular vectors are requested
614  IF ( conda ) THEN
615  minwrk = n + max( lwqp3, lwcon, lwsvd )
616  ELSE
617  minwrk = n + max( lwqp3, lwsvd )
618  END IF
619  IF ( lquery ) THEN
620  IF ( rtrans ) THEN
621  CALL dgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
622  $ v, ldv, rdummy, -1, ierr )
623  ELSE
624  CALL dgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
625  $ v, ldv, rdummy, -1, ierr )
626  END IF
627  lwrk_dgesvd = int( rdummy(1) )
628  IF ( conda ) THEN
629  optwrk = n + max( lwrk_dgeqp3, lwcon, lwrk_dgesvd )
630  ELSE
631  optwrk = n + max( lwrk_dgeqp3, lwrk_dgesvd )
632  END IF
633  END IF
634  ELSE
635 * .. minimal and optimal sizes of the workspace if the
636 * full SVD is requested
637  IF ( rtrans ) THEN
638  minwrk = max( lwqp3, lwsvd, lworq )
639  IF ( conda ) minwrk = max( minwrk, lwcon )
640  minwrk = minwrk + n
641  IF ( wntva ) THEN
642 * .. minimal workspace length for N x N/2 DGEQRF
643  lwqrf = max( n/2, 1 )
644 * .. minimal workspace lengt for N/2 x N/2 DGESVD
645  lwsvd2 = max( 5 * (n/2), 1 )
646  lworq2 = max( n, 1 )
647  minwrk2 = max( lwqp3, n/2+lwqrf, n/2+lwsvd2,
648  $ n/2+lworq2, lworq )
649  IF ( conda ) minwrk2 = max( minwrk2, lwcon )
650  minwrk2 = n + minwrk2
651  minwrk = max( minwrk, minwrk2 )
652  END IF
653  ELSE
654  minwrk = max( lwqp3, lwsvd, lworq )
655  IF ( conda ) minwrk = max( minwrk, lwcon )
656  minwrk = minwrk + n
657  IF ( wntva ) THEN
658 * .. minimal workspace length for N/2 x N DGELQF
659  lwlqf = max( n/2, 1 )
660  lwsvd2 = max( 5 * (n/2), 1 )
661  lworlq = max( n , 1 )
662  minwrk2 = max( lwqp3, n/2+lwlqf, n/2+lwsvd2,
663  $ n/2+lworlq, lworq )
664  IF ( conda ) minwrk2 = max( minwrk2, lwcon )
665  minwrk2 = n + minwrk2
666  minwrk = max( minwrk, minwrk2 )
667  END IF
668  END IF
669  IF ( lquery ) THEN
670  IF ( rtrans ) THEN
671  CALL dgesvd( 'O', 'A', n, n, a, lda, s, u, ldu,
672  $ v, ldv, rdummy, -1, ierr )
673  lwrk_dgesvd = int( rdummy(1) )
674  optwrk = max(lwrk_dgeqp3,lwrk_dgesvd,lwrk_dormqr)
675  IF ( conda ) optwrk = max( optwrk, lwcon )
676  optwrk = n + optwrk
677  IF ( wntva ) THEN
678  CALL dgeqrf(n,n/2,u,ldu,rdummy,rdummy,-1,ierr)
679  lwrk_dgeqrf = int( rdummy(1) )
680  CALL dgesvd( 'S', 'O', n/2,n/2, v,ldv, s, u,ldu,
681  $ v, ldv, rdummy, -1, ierr )
682  lwrk_dgesvd2 = int( rdummy(1) )
683  CALL dormqr( 'R', 'C', n, n, n/2, u, ldu, rdummy,
684  $ v, ldv, rdummy, -1, ierr )
685  lwrk_dormqr2 = int( rdummy(1) )
686  optwrk2 = max( lwrk_dgeqp3, n/2+lwrk_dgeqrf,
687  $ n/2+lwrk_dgesvd2, n/2+lwrk_dormqr2 )
688  IF ( conda ) optwrk2 = max( optwrk2, lwcon )
689  optwrk2 = n + optwrk2
690  optwrk = max( optwrk, optwrk2 )
691  END IF
692  ELSE
693  CALL dgesvd( 'S', 'O', n, n, a, lda, s, u, ldu,
694  $ v, ldv, rdummy, -1, ierr )
695  lwrk_dgesvd = int( rdummy(1) )
696  optwrk = max(lwrk_dgeqp3,lwrk_dgesvd,lwrk_dormqr)
697  IF ( conda ) optwrk = max( optwrk, lwcon )
698  optwrk = n + optwrk
699  IF ( wntva ) THEN
700  CALL dgelqf(n/2,n,u,ldu,rdummy,rdummy,-1,ierr)
701  lwrk_dgelqf = int( rdummy(1) )
702  CALL dgesvd( 'S','O', n/2,n/2, v, ldv, s, u, ldu,
703  $ v, ldv, rdummy, -1, ierr )
704  lwrk_dgesvd2 = int( rdummy(1) )
705  CALL dormlq( 'R', 'N', n, n, n/2, u, ldu, rdummy,
706  $ v, ldv, rdummy,-1,ierr )
707  lwrk_dormlq = int( rdummy(1) )
708  optwrk2 = max( lwrk_dgeqp3, n/2+lwrk_dgelqf,
709  $ n/2+lwrk_dgesvd2, n/2+lwrk_dormlq )
710  IF ( conda ) optwrk2 = max( optwrk2, lwcon )
711  optwrk2 = n + optwrk2
712  optwrk = max( optwrk, optwrk2 )
713  END IF
714  END IF
715  END IF
716  END IF
717 *
718  minwrk = max( 2, minwrk )
719  optwrk = max( 2, optwrk )
720  IF ( lwork .LT. minwrk .AND. (.NOT.lquery) ) info = -19
721 *
722  END IF
723 *
724  IF (info .EQ. 0 .AND. lrwork .LT. rminwrk .AND. .NOT. lquery) THEN
725  info = -21
726  END IF
727  IF( info.NE.0 ) THEN
728  CALL xerbla( 'DGESVDQ', -info )
729  RETURN
730  ELSE IF ( lquery ) THEN
731 *
732 * Return optimal workspace
733 *
734  iwork(1) = iminwrk
735  work(1) = optwrk
736  work(2) = minwrk
737  rwork(1) = rminwrk
738  RETURN
739  END IF
740 *
741 * Quick return if the matrix is void.
742 *
743  IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) ) THEN
744 * .. all output is void.
745  RETURN
746  END IF
747 *
748  big = dlamch('O')
749  ascaled = .false.
750  iwoff = 1
751  IF ( rowprm ) THEN
752  iwoff = m
753 * .. reordering the rows in decreasing sequence in the
754 * ell-infinity norm - this enhances numerical robustness in
755 * the case of differently scaled rows.
756  DO 1904 p = 1, m
757 * RWORK(p) = ABS( A(p,ICAMAX(N,A(p,1),LDA)) )
758 * [[DLANGE will return NaN if an entry of the p-th row is Nan]]
759  rwork(p) = dlange( 'M', 1, n, a(p,1), lda, rdummy )
760 * .. check for NaN's and Inf's
761  IF ( ( rwork(p) .NE. rwork(p) ) .OR.
762  $ ( (rwork(p)*zero) .NE. zero ) ) THEN
763  info = -8
764  CALL xerbla( 'DGESVDQ', -info )
765  RETURN
766  END IF
767  1904 CONTINUE
768  DO 1952 p = 1, m - 1
769  q = idamax( m-p+1, rwork(p), 1 ) + p - 1
770  iwork(n+p) = q
771  IF ( p .NE. q ) THEN
772  rtmp = rwork(p)
773  rwork(p) = rwork(q)
774  rwork(q) = rtmp
775  END IF
776  1952 CONTINUE
777 *
778  IF ( rwork(1) .EQ. zero ) THEN
779 * Quick return: A is the M x N zero matrix.
780  numrank = 0
781  CALL dlaset( 'G', n, 1, zero, zero, s, n )
782  IF ( wntus ) CALL dlaset('G', m, n, zero, one, u, ldu)
783  IF ( wntua ) CALL dlaset('G', m, m, zero, one, u, ldu)
784  IF ( wntva ) CALL dlaset('G', n, n, zero, one, v, ldv)
785  IF ( wntuf ) THEN
786  CALL dlaset( 'G', n, 1, zero, zero, work, n )
787  CALL dlaset( 'G', m, n, zero, one, u, ldu )
788  END IF
789  DO 5001 p = 1, n
790  iwork(p) = p
791  5001 CONTINUE
792  IF ( rowprm ) THEN
793  DO 5002 p = n + 1, n + m - 1
794  iwork(p) = p - n
795  5002 CONTINUE
796  END IF
797  IF ( conda ) rwork(1) = -1
798  rwork(2) = -1
799  RETURN
800  END IF
801 *
802  IF ( rwork(1) .GT. big / sqrt(dble(m)) ) THEN
803 * .. to prevent overflow in the QR factorization, scale the
804 * matrix by 1/sqrt(M) if too large entry detected
805  CALL dlascl('G',0,0,sqrt(dble(m)),one, m,n, a,lda, ierr)
806  ascaled = .true.
807  END IF
808  CALL dlaswp( n, a, lda, 1, m-1, iwork(n+1), 1 )
809  END IF
810 *
811 * .. At this stage, preemptive scaling is done only to avoid column
812 * norms overflows during the QR factorization. The SVD procedure should
813 * have its own scaling to save the singular values from overflows and
814 * underflows. That depends on the SVD procedure.
815 *
816  IF ( .NOT.rowprm ) THEN
817  rtmp = dlange( 'M', m, n, a, lda, rdummy )
818  IF ( ( rtmp .NE. rtmp ) .OR.
819  $ ( (rtmp*zero) .NE. zero ) ) THEN
820  info = -8
821  CALL xerbla( 'DGESVDQ', -info )
822  RETURN
823  END IF
824  IF ( rtmp .GT. big / sqrt(dble(m)) ) THEN
825 * .. to prevent overflow in the QR factorization, scale the
826 * matrix by 1/sqrt(M) if too large entry detected
827  CALL dlascl('G',0,0, sqrt(dble(m)),one, m,n, a,lda, ierr)
828  ascaled = .true.
829  END IF
830  END IF
831 *
832 * .. QR factorization with column pivoting
833 *
834 * A * P = Q * [ R ]
835 * [ 0 ]
836 *
837  DO 1963 p = 1, n
838 * .. all columns are free columns
839  iwork(p) = 0
840  1963 CONTINUE
841  CALL dgeqp3( m, n, a, lda, iwork, work, work(n+1), lwork-n,
842  $ ierr )
843 *
844 * If the user requested accuracy level allows truncation in the
845 * computed upper triangular factor, the matrix R is examined and,
846 * if possible, replaced with its leading upper trapezoidal part.
847 *
848  epsln = dlamch('E')
849  sfmin = dlamch('S')
850 * SMALL = SFMIN / EPSLN
851  nr = n
852 *
853  IF ( accla ) THEN
854 *
855 * Standard absolute error bound suffices. All sigma_i with
856 * sigma_i < N*EPS*||A||_F are flushed to zero. This is an
857 * aggressive enforcement of lower numerical rank by introducing a
858 * backward error of the order of N*EPS*||A||_F.
859  nr = 1
860  rtmp = sqrt(dble(n))*epsln
861  DO 3001 p = 2, n
862  IF ( abs(a(p,p)) .LT. (rtmp*abs(a(1,1))) ) GO TO 3002
863  nr = nr + 1
864  3001 CONTINUE
865  3002 CONTINUE
866 *
867  ELSEIF ( acclm ) THEN
868 * .. similarly as above, only slightly more gentle (less aggressive).
869 * Sudden drop on the diagonal of R is used as the criterion for being
870 * close-to-rank-deficient. The threshold is set to EPSLN=DLAMCH('E').
871 * [[This can be made more flexible by replacing this hard-coded value
872 * with a user specified threshold.]] Also, the values that underflow
873 * will be truncated.
874  nr = 1
875  DO 3401 p = 2, n
876  IF ( ( abs(a(p,p)) .LT. (epsln*abs(a(p-1,p-1))) ) .OR.
877  $ ( abs(a(p,p)) .LT. sfmin ) ) GO TO 3402
878  nr = nr + 1
879  3401 CONTINUE
880  3402 CONTINUE
881 *
882  ELSE
883 * .. RRQR not authorized to determine numerical rank except in the
884 * obvious case of zero pivots.
885 * .. inspect R for exact zeros on the diagonal;
886 * R(i,i)=0 => R(i:N,i:N)=0.
887  nr = 1
888  DO 3501 p = 2, n
889  IF ( abs(a(p,p)) .EQ. zero ) GO TO 3502
890  nr = nr + 1
891  3501 CONTINUE
892  3502 CONTINUE
893 *
894  IF ( conda ) THEN
895 * Estimate the scaled condition number of A. Use the fact that it is
896 * the same as the scaled condition number of R.
897 * .. V is used as workspace
898  CALL dlacpy( 'U', n, n, a, lda, v, ldv )
899 * Only the leading NR x NR submatrix of the triangular factor
900 * is considered. Only if NR=N will this give a reliable error
901 * bound. However, even for NR < N, this can be used on an
902 * expert level and obtain useful information in the sense of
903 * perturbation theory.
904  DO 3053 p = 1, nr
905  rtmp = dnrm2( p, v(1,p), 1 )
906  CALL dscal( p, one/rtmp, v(1,p), 1 )
907  3053 CONTINUE
908  IF ( .NOT. ( lsvec .OR. rsvec ) ) THEN
909  CALL dpocon( 'U', nr, v, ldv, one, rtmp,
910  $ work, iwork(n+iwoff), ierr )
911  ELSE
912  CALL dpocon( 'U', nr, v, ldv, one, rtmp,
913  $ work(n+1), iwork(n+iwoff), ierr )
914  END IF
915  sconda = one / sqrt(rtmp)
916 * For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1),
917 * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
918 * See the reference [1] for more details.
919  END IF
920 *
921  ENDIF
922 *
923  IF ( wntur ) THEN
924  n1 = nr
925  ELSE IF ( wntus .OR. wntuf) THEN
926  n1 = n
927  ELSE IF ( wntua ) THEN
928  n1 = m
929  END IF
930 *
931  IF ( .NOT. ( rsvec .OR. lsvec ) ) THEN
932 *.......................................................................
933 * .. only the singular values are requested
934 *.......................................................................
935  IF ( rtrans ) THEN
936 *
937 * .. compute the singular values of R**T = [A](1:NR,1:N)**T
938 * .. set the lower triangle of [A] to [A](1:NR,1:N)**T and
939 * the upper triangle of [A] to zero.
940  DO 1146 p = 1, min( n, nr )
941  DO 1147 q = p + 1, n
942  a(q,p) = a(p,q)
943  IF ( q .LE. nr ) a(p,q) = zero
944  1147 CONTINUE
945  1146 CONTINUE
946 *
947  CALL dgesvd( 'N', 'N', n, nr, a, lda, s, u, ldu,
948  $ v, ldv, work, lwork, info )
949 *
950  ELSE
951 *
952 * .. compute the singular values of R = [A](1:NR,1:N)
953 *
954  IF ( nr .GT. 1 )
955  $ CALL dlaset( 'L', nr-1,nr-1, zero,zero, a(2,1), lda )
956  CALL dgesvd( 'N', 'N', nr, n, a, lda, s, u, ldu,
957  $ v, ldv, work, lwork, info )
958 *
959  END IF
960 *
961  ELSE IF ( lsvec .AND. ( .NOT. rsvec) ) THEN
962 *.......................................................................
963 * .. the singular values and the left singular vectors requested
964 *.......................................................................""""""""
965  IF ( rtrans ) THEN
966 * .. apply DGESVD to R**T
967 * .. copy R**T into [U] and overwrite [U] with the right singular
968 * vectors of R
969  DO 1192 p = 1, nr
970  DO 1193 q = p, n
971  u(q,p) = a(p,q)
972  1193 CONTINUE
973  1192 CONTINUE
974  IF ( nr .GT. 1 )
975  $ CALL dlaset( 'U', nr-1,nr-1, zero,zero, u(1,2), ldu )
976 * .. the left singular vectors not computed, the NR right singular
977 * vectors overwrite [U](1:NR,1:NR) as transposed. These
978 * will be pre-multiplied by Q to build the left singular vectors of A.
979  CALL dgesvd( 'N', 'O', n, nr, u, ldu, s, u, ldu,
980  $ u, ldu, work(n+1), lwork-n, info )
981 *
982  DO 1119 p = 1, nr
983  DO 1120 q = p + 1, nr
984  rtmp = u(q,p)
985  u(q,p) = u(p,q)
986  u(p,q) = rtmp
987  1120 CONTINUE
988  1119 CONTINUE
989 *
990  ELSE
991 * .. apply DGESVD to R
992 * .. copy R into [U] and overwrite [U] with the left singular vectors
993  CALL dlacpy( 'U', nr, n, a, lda, u, ldu )
994  IF ( nr .GT. 1 )
995  $ CALL dlaset( 'L', nr-1, nr-1, zero, zero, u(2,1), ldu )
996 * .. the right singular vectors not computed, the NR left singular
997 * vectors overwrite [U](1:NR,1:NR)
998  CALL dgesvd( 'O', 'N', nr, n, u, ldu, s, u, ldu,
999  $ v, ldv, work(n+1), lwork-n, info )
1000 * .. now [U](1:NR,1:NR) contains the NR left singular vectors of
1001 * R. These will be pre-multiplied by Q to build the left singular
1002 * vectors of A.
1003  END IF
1004 *
1005 * .. assemble the left singular vector matrix U of dimensions
1006 * (M x NR) or (M x N) or (M x M).
1007  IF ( ( nr .LT. m ) .AND. ( .NOT.wntuf ) ) THEN
1008  CALL dlaset('A', m-nr, nr, zero, zero, u(nr+1,1), ldu)
1009  IF ( nr .LT. n1 ) THEN
1010  CALL dlaset( 'A',nr,n1-nr,zero,zero,u(1,nr+1), ldu )
1011  CALL dlaset( 'A',m-nr,n1-nr,zero,one,
1012  $ u(nr+1,nr+1), ldu )
1013  END IF
1014  END IF
1015 *
1016 * The Q matrix from the first QRF is built into the left singular
1017 * vectors matrix U.
1018 *
1019  IF ( .NOT.wntuf )
1020  $ CALL dormqr( 'L', 'N', m, n1, n, a, lda, work, u,
1021  $ ldu, work(n+1), lwork-n, ierr )
1022  IF ( rowprm .AND. .NOT.wntuf )
1023  $ CALL dlaswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
1024 *
1025  ELSE IF ( rsvec .AND. ( .NOT. lsvec ) ) THEN
1026 *.......................................................................
1027 * .. the singular values and the right singular vectors requested
1028 *.......................................................................
1029  IF ( rtrans ) THEN
1030 * .. apply DGESVD to R**T
1031 * .. copy R**T into V and overwrite V with the left singular vectors
1032  DO 1165 p = 1, nr
1033  DO 1166 q = p, n
1034  v(q,p) = (a(p,q))
1035  1166 CONTINUE
1036  1165 CONTINUE
1037  IF ( nr .GT. 1 )
1038  $ CALL dlaset( 'U', nr-1,nr-1, zero,zero, v(1,2), ldv )
1039 * .. the left singular vectors of R**T overwrite V, the right singular
1040 * vectors not computed
1041  IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1042  CALL dgesvd( 'O', 'N', n, nr, v, ldv, s, u, ldu,
1043  $ u, ldu, work(n+1), lwork-n, info )
1044 *
1045  DO 1121 p = 1, nr
1046  DO 1122 q = p + 1, nr
1047  rtmp = v(q,p)
1048  v(q,p) = v(p,q)
1049  v(p,q) = rtmp
1050  1122 CONTINUE
1051  1121 CONTINUE
1052 *
1053  IF ( nr .LT. n ) THEN
1054  DO 1103 p = 1, nr
1055  DO 1104 q = nr + 1, n
1056  v(p,q) = v(q,p)
1057  1104 CONTINUE
1058  1103 CONTINUE
1059  END IF
1060  CALL dlapmt( .false., nr, n, v, ldv, iwork )
1061  ELSE
1062 * .. need all N right singular vectors and NR < N
1063 * [!] This is simple implementation that augments [V](1:N,1:NR)
1064 * by padding a zero block. In the case NR << N, a more efficient
1065 * way is to first use the QR factorization. For more details
1066 * how to implement this, see the " FULL SVD " branch.
1067  CALL dlaset('G', n, n-nr, zero, zero, v(1,nr+1), ldv)
1068  CALL dgesvd( 'O', 'N', n, n, v, ldv, s, u, ldu,
1069  $ u, ldu, work(n+1), lwork-n, info )
1070 *
1071  DO 1123 p = 1, n
1072  DO 1124 q = p + 1, n
1073  rtmp = v(q,p)
1074  v(q,p) = v(p,q)
1075  v(p,q) = rtmp
1076  1124 CONTINUE
1077  1123 CONTINUE
1078  CALL dlapmt( .false., n, n, v, ldv, iwork )
1079  END IF
1080 *
1081  ELSE
1082 * .. aply DGESVD to R
1083 * .. copy R into V and overwrite V with the right singular vectors
1084  CALL dlacpy( 'U', nr, n, a, lda, v, ldv )
1085  IF ( nr .GT. 1 )
1086  $ CALL dlaset( 'L', nr-1, nr-1, zero, zero, v(2,1), ldv )
1087 * .. the right singular vectors overwrite V, the NR left singular
1088 * vectors stored in U(1:NR,1:NR)
1089  IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1090  CALL dgesvd( 'N', 'O', nr, n, v, ldv, s, u, ldu,
1091  $ v, ldv, work(n+1), lwork-n, info )
1092  CALL dlapmt( .false., nr, n, v, ldv, iwork )
1093 * .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T
1094  ELSE
1095 * .. need all N right singular vectors and NR < N
1096 * [!] This is simple implementation that augments [V](1:NR,1:N)
1097 * by padding a zero block. In the case NR << N, a more efficient
1098 * way is to first use the LQ factorization. For more details
1099 * how to implement this, see the " FULL SVD " branch.
1100  CALL dlaset('G', n-nr, n, zero,zero, v(nr+1,1), ldv)
1101  CALL dgesvd( 'N', 'O', n, n, v, ldv, s, u, ldu,
1102  $ v, ldv, work(n+1), lwork-n, info )
1103  CALL dlapmt( .false., n, n, v, ldv, iwork )
1104  END IF
1105 * .. now [V] contains the transposed matrix of the right singular
1106 * vectors of A.
1107  END IF
1108 *
1109  ELSE
1110 *.......................................................................
1111 * .. FULL SVD requested
1112 *.......................................................................
1113  IF ( rtrans ) THEN
1114 *
1115 * .. apply DGESVD to R**T [[this option is left for R&D&T]]
1116 *
1117  IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1118 * .. copy R**T into [V] and overwrite [V] with the left singular
1119 * vectors of R**T
1120  DO 1168 p = 1, nr
1121  DO 1169 q = p, n
1122  v(q,p) = a(p,q)
1123  1169 CONTINUE
1124  1168 CONTINUE
1125  IF ( nr .GT. 1 )
1126  $ CALL dlaset( 'U', nr-1,nr-1, zero,zero, v(1,2), ldv )
1127 *
1128 * .. the left singular vectors of R**T overwrite [V], the NR right
1129 * singular vectors of R**T stored in [U](1:NR,1:NR) as transposed
1130  CALL dgesvd( 'O', 'A', n, nr, v, ldv, s, v, ldv,
1131  $ u, ldu, work(n+1), lwork-n, info )
1132 * .. assemble V
1133  DO 1115 p = 1, nr
1134  DO 1116 q = p + 1, nr
1135  rtmp = v(q,p)
1136  v(q,p) = v(p,q)
1137  v(p,q) = rtmp
1138  1116 CONTINUE
1139  1115 CONTINUE
1140  IF ( nr .LT. n ) THEN
1141  DO 1101 p = 1, nr
1142  DO 1102 q = nr+1, n
1143  v(p,q) = v(q,p)
1144  1102 CONTINUE
1145  1101 CONTINUE
1146  END IF
1147  CALL dlapmt( .false., nr, n, v, ldv, iwork )
1148 *
1149  DO 1117 p = 1, nr
1150  DO 1118 q = p + 1, nr
1151  rtmp = u(q,p)
1152  u(q,p) = u(p,q)
1153  u(p,q) = rtmp
1154  1118 CONTINUE
1155  1117 CONTINUE
1156 *
1157  IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1158  CALL dlaset('A', m-nr,nr, zero,zero, u(nr+1,1), ldu)
1159  IF ( nr .LT. n1 ) THEN
1160  CALL dlaset('A',nr,n1-nr,zero,zero,u(1,nr+1),ldu)
1161  CALL dlaset( 'A',m-nr,n1-nr,zero,one,
1162  $ u(nr+1,nr+1), ldu )
1163  END IF
1164  END IF
1165 *
1166  ELSE
1167 * .. need all N right singular vectors and NR < N
1168 * .. copy R**T into [V] and overwrite [V] with the left singular
1169 * vectors of R**T
1170 * [[The optimal ratio N/NR for using QRF instead of padding
1171 * with zeros. Here hard coded to 2; it must be at least
1172 * two due to work space constraints.]]
1173 * OPTRATIO = ILAENV(6, 'DGESVD', 'S' // 'O', NR,N,0,0)
1174 * OPTRATIO = MAX( OPTRATIO, 2 )
1175  optratio = 2
1176  IF ( optratio*nr .GT. n ) THEN
1177  DO 1198 p = 1, nr
1178  DO 1199 q = p, n
1179  v(q,p) = a(p,q)
1180  1199 CONTINUE
1181  1198 CONTINUE
1182  IF ( nr .GT. 1 )
1183  $ CALL dlaset('U',nr-1,nr-1, zero,zero, v(1,2),ldv)
1184 *
1185  CALL dlaset('A',n,n-nr,zero,zero,v(1,nr+1),ldv)
1186  CALL dgesvd( 'O', 'A', n, n, v, ldv, s, v, ldv,
1187  $ u, ldu, work(n+1), lwork-n, info )
1188 *
1189  DO 1113 p = 1, n
1190  DO 1114 q = p + 1, n
1191  rtmp = v(q,p)
1192  v(q,p) = v(p,q)
1193  v(p,q) = rtmp
1194  1114 CONTINUE
1195  1113 CONTINUE
1196  CALL dlapmt( .false., n, n, v, ldv, iwork )
1197 * .. assemble the left singular vector matrix U of dimensions
1198 * (M x N1), i.e. (M x N) or (M x M).
1199 *
1200  DO 1111 p = 1, n
1201  DO 1112 q = p + 1, n
1202  rtmp = u(q,p)
1203  u(q,p) = u(p,q)
1204  u(p,q) = rtmp
1205  1112 CONTINUE
1206  1111 CONTINUE
1207 *
1208  IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) THEN
1209  CALL dlaset('A',m-n,n,zero,zero,u(n+1,1),ldu)
1210  IF ( n .LT. n1 ) THEN
1211  CALL dlaset('A',n,n1-n,zero,zero,u(1,n+1),ldu)
1212  CALL dlaset('A',m-n,n1-n,zero,one,
1213  $ u(n+1,n+1), ldu )
1214  END IF
1215  END IF
1216  ELSE
1217 * .. copy R**T into [U] and overwrite [U] with the right
1218 * singular vectors of R
1219  DO 1196 p = 1, nr
1220  DO 1197 q = p, n
1221  u(q,nr+p) = a(p,q)
1222  1197 CONTINUE
1223  1196 CONTINUE
1224  IF ( nr .GT. 1 )
1225  $ CALL dlaset('U',nr-1,nr-1,zero,zero,u(1,nr+2),ldu)
1226  CALL dgeqrf( n, nr, u(1,nr+1), ldu, work(n+1),
1227  $ work(n+nr+1), lwork-n-nr, ierr )
1228  DO 1143 p = 1, nr
1229  DO 1144 q = 1, n
1230  v(q,p) = u(p,nr+q)
1231  1144 CONTINUE
1232  1143 CONTINUE
1233  CALL dlaset('U',nr-1,nr-1,zero,zero,v(1,2),ldv)
1234  CALL dgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
1235  $ v,ldv, work(n+nr+1),lwork-n-nr, info )
1236  CALL dlaset('A',n-nr,nr,zero,zero,v(nr+1,1),ldv)
1237  CALL dlaset('A',nr,n-nr,zero,zero,v(1,nr+1),ldv)
1238  CALL dlaset('A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv)
1239  CALL dormqr('R','C', n, n, nr, u(1,nr+1), ldu,
1240  $ work(n+1),v,ldv,work(n+nr+1),lwork-n-nr,ierr)
1241  CALL dlapmt( .false., n, n, v, ldv, iwork )
1242 * .. assemble the left singular vector matrix U of dimensions
1243 * (M x NR) or (M x N) or (M x M).
1244  IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1245  CALL dlaset('A',m-nr,nr,zero,zero,u(nr+1,1),ldu)
1246  IF ( nr .LT. n1 ) THEN
1247  CALL dlaset('A',nr,n1-nr,zero,zero,u(1,nr+1),ldu)
1248  CALL dlaset( 'A',m-nr,n1-nr,zero,one,
1249  $ u(nr+1,nr+1),ldu)
1250  END IF
1251  END IF
1252  END IF
1253  END IF
1254 *
1255  ELSE
1256 *
1257 * .. apply DGESVD to R [[this is the recommended option]]
1258 *
1259  IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1260 * .. copy R into [V] and overwrite V with the right singular vectors
1261  CALL dlacpy( 'U', nr, n, a, lda, v, ldv )
1262  IF ( nr .GT. 1 )
1263  $ CALL dlaset( 'L', nr-1,nr-1, zero,zero, v(2,1), ldv )
1264 * .. the right singular vectors of R overwrite [V], the NR left
1265 * singular vectors of R stored in [U](1:NR,1:NR)
1266  CALL dgesvd( 'S', 'O', nr, n, v, ldv, s, u, ldu,
1267  $ v, ldv, work(n+1), lwork-n, info )
1268  CALL dlapmt( .false., nr, n, v, ldv, iwork )
1269 * .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T
1270 * .. assemble the left singular vector matrix U of dimensions
1271 * (M x NR) or (M x N) or (M x M).
1272  IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1273  CALL dlaset('A', m-nr,nr, zero,zero, u(nr+1,1), ldu)
1274  IF ( nr .LT. n1 ) THEN
1275  CALL dlaset('A',nr,n1-nr,zero,zero,u(1,nr+1),ldu)
1276  CALL dlaset( 'A',m-nr,n1-nr,zero,one,
1277  $ u(nr+1,nr+1), ldu )
1278  END IF
1279  END IF
1280 *
1281  ELSE
1282 * .. need all N right singular vectors and NR < N
1283 * .. the requested number of the left singular vectors
1284 * is then N1 (N or M)
1285 * [[The optimal ratio N/NR for using LQ instead of padding
1286 * with zeros. Here hard coded to 2; it must be at least
1287 * two due to work space constraints.]]
1288 * OPTRATIO = ILAENV(6, 'DGESVD', 'S' // 'O', NR,N,0,0)
1289 * OPTRATIO = MAX( OPTRATIO, 2 )
1290  optratio = 2
1291  IF ( optratio * nr .GT. n ) THEN
1292  CALL dlacpy( 'U', nr, n, a, lda, v, ldv )
1293  IF ( nr .GT. 1 )
1294  $ CALL dlaset('L', nr-1,nr-1, zero,zero, v(2,1),ldv)
1295 * .. the right singular vectors of R overwrite [V], the NR left
1296 * singular vectors of R stored in [U](1:NR,1:NR)
1297  CALL dlaset('A', n-nr,n, zero,zero, v(nr+1,1),ldv)
1298  CALL dgesvd( 'S', 'O', n, n, v, ldv, s, u, ldu,
1299  $ v, ldv, work(n+1), lwork-n, info )
1300  CALL dlapmt( .false., n, n, v, ldv, iwork )
1301 * .. now [V] contains the transposed matrix of the right
1302 * singular vectors of A. The leading N left singular vectors
1303 * are in [U](1:N,1:N)
1304 * .. assemble the left singular vector matrix U of dimensions
1305 * (M x N1), i.e. (M x N) or (M x M).
1306  IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) THEN
1307  CALL dlaset('A',m-n,n,zero,zero,u(n+1,1),ldu)
1308  IF ( n .LT. n1 ) THEN
1309  CALL dlaset('A',n,n1-n,zero,zero,u(1,n+1),ldu)
1310  CALL dlaset( 'A',m-n,n1-n,zero,one,
1311  $ u(n+1,n+1), ldu )
1312  END IF
1313  END IF
1314  ELSE
1315  CALL dlacpy( 'U', nr, n, a, lda, u(nr+1,1), ldu )
1316  IF ( nr .GT. 1 )
1317  $ CALL dlaset('L',nr-1,nr-1,zero,zero,u(nr+2,1),ldu)
1318  CALL dgelqf( nr, n, u(nr+1,1), ldu, work(n+1),
1319  $ work(n+nr+1), lwork-n-nr, ierr )
1320  CALL dlacpy('L',nr,nr,u(nr+1,1),ldu,v,ldv)
1321  IF ( nr .GT. 1 )
1322  $ CALL dlaset('U',nr-1,nr-1,zero,zero,v(1,2),ldv)
1323  CALL dgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
1324  $ v, ldv, work(n+nr+1), lwork-n-nr, info )
1325  CALL dlaset('A',n-nr,nr,zero,zero,v(nr+1,1),ldv)
1326  CALL dlaset('A',nr,n-nr,zero,zero,v(1,nr+1),ldv)
1327  CALL dlaset('A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv)
1328  CALL dormlq('R','N',n,n,nr,u(nr+1,1),ldu,work(n+1),
1329  $ v, ldv, work(n+nr+1),lwork-n-nr,ierr)
1330  CALL dlapmt( .false., n, n, v, ldv, iwork )
1331 * .. assemble the left singular vector matrix U of dimensions
1332 * (M x NR) or (M x N) or (M x M).
1333  IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1334  CALL dlaset('A',m-nr,nr,zero,zero,u(nr+1,1),ldu)
1335  IF ( nr .LT. n1 ) THEN
1336  CALL dlaset('A',nr,n1-nr,zero,zero,u(1,nr+1),ldu)
1337  CALL dlaset( 'A',m-nr,n1-nr,zero,one,
1338  $ u(nr+1,nr+1), ldu )
1339  END IF
1340  END IF
1341  END IF
1342  END IF
1343 * .. end of the "R**T or R" branch
1344  END IF
1345 *
1346 * The Q matrix from the first QRF is built into the left singular
1347 * vectors matrix U.
1348 *
1349  IF ( .NOT. wntuf )
1350  $ CALL dormqr( 'L', 'N', m, n1, n, a, lda, work, u,
1351  $ ldu, work(n+1), lwork-n, ierr )
1352  IF ( rowprm .AND. .NOT.wntuf )
1353  $ CALL dlaswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
1354 *
1355 * ... end of the "full SVD" branch
1356  END IF
1357 *
1358 * Check whether some singular values are returned as zeros, e.g.
1359 * due to underflow, and update the numerical rank.
1360  p = nr
1361  DO 4001 q = p, 1, -1
1362  IF ( s(q) .GT. zero ) GO TO 4002
1363  nr = nr - 1
1364  4001 CONTINUE
1365  4002 CONTINUE
1366 *
1367 * .. if numerical rank deficiency is detected, the truncated
1368 * singular values are set to zero.
1369  IF ( nr .LT. n ) CALL dlaset( 'G', n-nr,1, zero,zero, s(nr+1), n )
1370 * .. undo scaling; this may cause overflow in the largest singular
1371 * values.
1372  IF ( ascaled )
1373  $ CALL dlascl( 'G',0,0, one,sqrt(dble(m)), nr,1, s, n, ierr )
1374  IF ( conda ) rwork(1) = sconda
1375  rwork(2) = p - nr
1376 * .. p-NR is the number of singular values that are computed as
1377 * exact zeros in DGESVD() applied to the (possibly truncated)
1378 * full row rank triangular (trapezoidal) factor of A.
1379  numrank = nr
1380 *
1381  RETURN
1382 *
1383 * End of DGESVDQ
1384 *
1385  END
dlascl
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:145
dormlq
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
Definition: dormlq.f:169
dlapmt
subroutine dlapmt(FORWRD, M, N, X, LDX, K)
DLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: dlapmt.f:106
dgelqf
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
Definition: dgelqf.f:145
dlaswp
subroutine dlaswp(N, A, LDA, K1, K2, IPIV, INCX)
DLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: dlaswp.f:117
dgesvd
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:213
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dgesvdq
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:417
dgeqp3
subroutine dgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
DGEQP3
Definition: dgeqp3.f:153
dlaset
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:112
dormqr
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:169
dgeqrf
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:147
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81
dpocon
subroutine dpocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
DPOCON
Definition: dpocon.f:123
dlacpy
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105