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