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