LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
cgesvdq.f
Go to the documentation of this file.
1 *> \brief <b> CGESVDQ 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 CGESVDQ + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesvdq.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesvdq.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesvdq.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
22 * S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
23 * CWORK, LCWORK, 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, LCWORK, LRWORK,
29 * INFO
30 * ..
31 * .. Array Arguments ..
32 * COMPLEX A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * )
33 * REAL S( * ), RWORK( * )
34 * INTEGER IWORK( * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> CGESVDQ computes the singular value decomposition (SVD) of a complex
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 unitary 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, CGESVD is applied to
92 *> the adjoint R**H 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 CGESVD. 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**H , 0)**H. 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 COMPLEX 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 CGEQP3. If JOBU = 'F', these Householder
145 *> vectors together with CWORK(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 COMPLEX 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 unitary 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 COMPLEX 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 unitary matrix V**H;
188 *> If JOBV = 'R', V contains the first NUMRANK rows of V**H (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 CGESVD. 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, LCWORK, 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';
231 *> LIWORK >= N if JOBP = 'N'.
232 *>
233 *> If LIWORK = -1, then a workspace query is assumed; the routine
234 *> only calculates and returns the optimal and minimal sizes
235 *> for the CWORK, IWORK, and RWORK arrays, and no error
236 *> message related to LCWORK is issued by XERBLA.
237 *> \endverbatim
238 *>
239 *> \param[out] CWORK
240 *> \verbatim
241 *> CWORK is COMPLEX array, dimension (max(2, LCWORK)), used as a workspace.
242 *> On exit, if, on entry, LCWORK.NE.-1, CWORK(1:N) contains parameters
243 *> needed to recover the Q factor from the QR factorization computed by
244 *> CGEQP3.
245 *>
246 *> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
247 *> CWORK(1) returns the optimal LCWORK, and
248 *> CWORK(2) returns the minimal LCWORK.
249 *> \endverbatim
250 *>
251 *> \param[in,out] LCWORK
252 *> \verbatim
253 *> LCWORK is INTEGER
254 *> The dimension of the array CWORK. It is determined as follows:
255 *> Let LWQP3 = N+1, LWCON = 2*N, and let
256 *> LWUNQ = { MAX( N, 1 ), if JOBU = 'R', 'S', or 'U'
257 *> { MAX( M, 1 ), if JOBU = 'A'
258 *> LWSVD = MAX( 3*N, 1 )
259 *> LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 3*(N/2), 1 ), LWUNLQ = MAX( N, 1 ),
260 *> LWQRF = MAX( N/2, 1 ), LWUNQ2 = MAX( N, 1 )
261 *> Then the minimal value of LCWORK is:
262 *> = MAX( N + LWQP3, LWSVD ) if only the singular values are needed;
263 *> = MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values are needed,
264 *> and a scaled condition estimate requested;
265 *>
266 *> = N + MAX( LWQP3, LWSVD, LWUNQ ) if the singular values and the left
267 *> singular vectors are requested;
268 *> = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the singular values and the left
269 *> singular vectors are requested, and also
270 *> a scaled condition estimate requested;
271 *>
272 *> = N + MAX( LWQP3, LWSVD ) if the singular values and the right
273 *> singular vectors are requested;
274 *> = N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and the right
275 *> singular vectors are requested, and also
276 *> a scaled condition etimate requested;
277 *>
278 *> = N + MAX( LWQP3, LWSVD, LWUNQ ) if the full SVD is requested with JOBV = 'R';
279 *> independent of JOBR;
280 *> = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the full SVD is requested,
281 *> JOBV = 'R' and, also a scaled condition
282 *> estimate requested; independent of JOBR;
283 *> = MAX( N + MAX( LWQP3, LWSVD, LWUNQ ),
284 *> N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ) ) if the
285 *> full SVD is requested with JOBV = 'A' or 'V', and
286 *> JOBR ='N'
287 *> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
288 *> N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ ) )
289 *> if the full SVD is requested with JOBV = 'A' or 'V', and
290 *> JOBR ='N', and also a scaled condition number estimate
291 *> requested.
292 *> = MAX( N + MAX( LWQP3, LWSVD, LWUNQ ),
293 *> N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) ) if the
294 *> full SVD is requested with JOBV = 'A', 'V', and JOBR ='T'
295 *> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
296 *> N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) )
297 *> if the full SVD is requested with JOBV = 'A', 'V' and
298 *> JOBR ='T', and also a scaled condition number estimate
299 *> requested.
300 *> Finally, LCWORK must be at least two: LCWORK = MAX( 2, LCWORK ).
301 *>
302 *> If LCWORK = -1, then a workspace query is assumed; the routine
303 *> only calculates and returns the optimal and minimal sizes
304 *> for the CWORK, IWORK, and RWORK arrays, and no error
305 *> message related to LCWORK is issued by XERBLA.
306 *> \endverbatim
307 *>
308 *> \param[out] RWORK
309 *> \verbatim
310 *> RWORK is REAL array, dimension (max(1, LRWORK)).
311 *> On exit,
312 *> 1. If JOBA = 'E', RWORK(1) contains an estimate of the condition
313 *> number of column scaled A. If A = C * D where D is diagonal and C
314 *> has unit columns in the Euclidean norm, then, assuming full column rank,
315 *> N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1).
316 *> Otherwise, RWORK(1) = -1.
317 *> 2. RWORK(2) contains the number of singular values computed as
318 *> exact zeros in CGESVD applied to the upper triangular or trapeziodal
319 *> R (from the initial QR factorization). In case of early exit (no call to
320 *> CGESVD, such as in the case of zero matrix) RWORK(2) = -1.
321 *>
322 *> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
323 *> RWORK(1) returns the minimal LRWORK.
324 *> \endverbatim
325 *>
326 *> \param[in] LRWORK
327 *> \verbatim
328 *> LRWORK is INTEGER.
329 *> The dimension of the array RWORK.
330 *> If JOBP ='P', then LRWORK >= MAX(2, M, 5*N);
331 *> Otherwise, LRWORK >= MAX(2, 5*N).
332 *>
333 *> If LRWORK = -1, then a workspace query is assumed; the routine
334 *> only calculates and returns the optimal and minimal sizes
335 *> for the CWORK, IWORK, and RWORK arrays, and no error
336 *> message related to LCWORK is issued by XERBLA.
337 *> \endverbatim
338 *>
339 *> \param[out] INFO
340 *> \verbatim
341 *> INFO is INTEGER
342 *> = 0: successful exit.
343 *> < 0: if INFO = -i, the i-th argument had an illegal value.
344 *> > 0: if CBDSQR did not converge, INFO specifies how many superdiagonals
345 *> of an intermediate bidiagonal form B (computed in CGESVD) did not
346 *> converge to zero.
347 *> \endverbatim
348 *
349 *> \par Further Details:
350 * ========================
351 *>
352 *> \verbatim
353 *>
354 *> 1. The data movement (matrix transpose) is coded using simple nested
355 *> DO-loops because BLAS and LAPACK do not provide corresponding subroutines.
356 *> Those DO-loops are easily identified in this source code - by the CONTINUE
357 *> statements labeled with 11**. In an optimized version of this code, the
358 *> nested DO loops should be replaced with calls to an optimized subroutine.
359 *> 2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause
360 *> column norm overflow. This is the minial precaution and it is left to the
361 *> SVD routine (CGESVD) to do its own preemptive scaling if potential over-
362 *> or underflows are detected. To avoid repeated scanning of the array A,
363 *> an optimal implementation would do all necessary scaling before calling
364 *> CGESVD and the scaling in CGESVD can be switched off.
365 *> 3. Other comments related to code optimization are given in comments in the
366 *> code, enlosed in [[double brackets]].
367 *> \endverbatim
368 *
369 *> \par Bugs, examples and comments
370 * ===========================
371 *
372 *> \verbatim
373 *> Please report all bugs and send interesting examples and/or comments to
374 *> drmac@math.hr. Thank you.
375 *> \endverbatim
376 *
377 *> \par References
378 * ===============
379 *
380 *> \verbatim
381 *> [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
382 *> Computing the SVD with High Accuracy. ACM Trans. Math. Softw.
383 *> 44(1): 11:1-11:30 (2017)
384 *>
385 *> SIGMA library, xGESVDQ section updated February 2016.
386 *> Developed and coded by Zlatko Drmac, Department of Mathematics
387 *> University of Zagreb, Croatia, drmac@math.hr
388 *> \endverbatim
389 *
390 *
391 *> \par Contributors:
392 * ==================
393 *>
394 *> \verbatim
395 *> Developed and coded by Zlatko Drmac, Department of Mathematics
396 *> University of Zagreb, Croatia, drmac@math.hr
397 *> \endverbatim
398 *
399 * Authors:
400 * ========
401 *
402 *> \author Univ. of Tennessee
403 *> \author Univ. of California Berkeley
404 *> \author Univ. of Colorado Denver
405 *> \author NAG Ltd.
406 *
407 *> \date November 2018
408 *
409 *> \ingroup complexGEsing
410 *
411 * =====================================================================
412  SUBROUTINE cgesvdq( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
413  $ S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
414  $ CWORK, LCWORK, 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, LCWORK, LRWORK,
419  $ info
420 * ..
421 * .. Array Arguments ..
422  COMPLEX A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * )
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  COMPLEX CZERO, CONE
432  parameter( czero = ( 0.0e0, 0.0e0 ), cone = ( 1.0e0, 0.0e0 ) )
433 * ..
434 * .. Local Scalars ..
435  INTEGER IERR, NR, N1, OPTRATIO, p, q
436  INTEGER LWCON, LWQP3, LWRK_CGELQF, LWRK_CGESVD, LWRK_CGESVD2,
437  $ lwrk_cgeqp3, lwrk_cgeqrf, lwrk_cunmlq, lwrk_cunmqr,
438  $ lwrk_cunmqr2, lwlqf, lwqrf, lwsvd, lwsvd2, lwunq,
439  $ lwunq2, 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  COMPLEX CTMP
446 * ..
447 * .. Local Arrays
448  COMPLEX CDUMMY(1)
449  REAL RDUMMY(1)
450 * ..
451 * .. External Subroutines (BLAS, LAPACK)
452  EXTERNAL cgelqf, cgeqp3, cgeqrf, cgesvd, clacpy, clapmt,
453  $ clascl, claset, claswp, csscal, slaset, slascl,
454  $ cpocon, cunmlq, cunmqr, xerbla
455 * ..
456 * .. External Functions (BLAS, LAPACK)
457  LOGICAL LSAME
458  INTEGER ISAMAX
459  REAL CLANGE, SCNRM2, SLAMCH
460  EXTERNAL clange, lsame, isamax, scnrm2, slamch
461 * ..
462 * .. Intrinsic Functions ..
463  INTRINSIC abs, conjg, max, min, real, sqrt
464 * ..
465 * .. Executable Statements ..
466 *
467 * Test the input arguments
468 *
469  wntus = lsame( jobu, 'S' ) .OR. lsame( jobu, 'U' )
470  wntur = lsame( jobu, 'R' )
471  wntua = lsame( jobu, 'A' )
472  wntuf = lsame( jobu, 'F' )
473  lsvc0 = wntus .OR. wntur .OR. wntua
474  lsvec = lsvc0 .OR. wntuf
475  dntwu = lsame( jobu, 'N' )
476 *
477  wntvr = lsame( jobv, 'R' )
478  wntva = lsame( jobv, 'A' ) .OR. lsame( jobv, 'V' )
479  rsvec = wntvr .OR. wntva
480  dntwv = lsame( jobv, 'N' )
481 *
482  accla = lsame( joba, 'A' )
483  acclm = lsame( joba, 'M' )
484  conda = lsame( joba, 'E' )
485  acclh = lsame( joba, 'H' ) .OR. conda
486 *
487  rowprm = lsame( jobp, 'P' )
488  rtrans = lsame( jobr, 'T' )
489 *
490  IF ( rowprm ) THEN
491  iminwrk = max( 1, n + m - 1 )
492  rminwrk = max( 2, m, 5*n )
493  ELSE
494  iminwrk = max( 1, n )
495  rminwrk = max( 2, 5*n )
496  END IF
497  lquery = (liwork .EQ. -1 .OR. lcwork .EQ. -1 .OR. lrwork .EQ. -1)
498  info = 0
499  IF ( .NOT. ( accla .OR. acclm .OR. acclh ) ) THEN
500  info = -1
501  ELSE IF ( .NOT.( rowprm .OR. lsame( jobp, 'N' ) ) ) THEN
502  info = -2
503  ELSE IF ( .NOT.( rtrans .OR. lsame( jobr, 'N' ) ) ) THEN
504  info = -3
505  ELSE IF ( .NOT.( lsvec .OR. dntwu ) ) THEN
506  info = -4
507  ELSE IF ( wntur .AND. wntva ) THEN
508  info = -5
509  ELSE IF ( .NOT.( rsvec .OR. dntwv )) THEN
510  info = -5
511  ELSE IF ( m.LT.0 ) THEN
512  info = -6
513  ELSE IF ( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
514  info = -7
515  ELSE IF ( lda.LT.max( 1, m ) ) THEN
516  info = -9
517  ELSE IF ( ldu.LT.1 .OR. ( lsvc0 .AND. ldu.LT.m ) .OR.
518  $ ( wntuf .AND. ldu.LT.n ) ) THEN
519  info = -12
520  ELSE IF ( ldv.LT.1 .OR. ( rsvec .AND. ldv.LT.n ) .OR.
521  $ ( conda .AND. ldv.LT.n ) ) THEN
522  info = -14
523  ELSE IF ( liwork .LT. iminwrk .AND. .NOT. lquery ) THEN
524  info = -17
525  END IF
526 *
527 *
528  IF ( info .EQ. 0 ) THEN
529 *
530 * Compute workspace
531 * .. compute the minimal and the optimal workspace lengths
532 * [[The expressions for computing the minimal and the optimal
533 * values of LCWORK are written with a lot of redundancy and
534 * can be simplified. However, this detailed form is easier for
535 * maintenance and modifications of the code.]]
536 *
537 * .. minimal workspace length for CGEQP3 of an M x N matrix
538  lwqp3 = n+1
539 * .. minimal workspace length for CUNMQR to build left singular vectors
540  IF ( wntus .OR. wntur ) THEN
541  lwunq = max( n , 1 )
542  ELSE IF ( wntua ) THEN
543  lwunq = max( m , 1 )
544  END IF
545 * .. minimal workspace length for CPOCON of an N x N matrix
546  lwcon = 2 * n
547 * .. CGESVD of an N x N matrix
548  lwsvd = max( 3 * n, 1 )
549  IF ( lquery ) THEN
550  CALL cgeqp3( m, n, a, lda, iwork, cdummy, cdummy, -1,
551  $ rdummy, ierr )
552  lwrk_cgeqp3 = int( cdummy(1) )
553  IF ( wntus .OR. wntur ) THEN
554  CALL cunmqr( 'L', 'N', m, n, n, a, lda, cdummy, u,
555  $ ldu, cdummy, -1, ierr )
556  lwrk_cunmqr = int( cdummy(1) )
557  ELSE IF ( wntua ) THEN
558  CALL cunmqr( 'L', 'N', m, m, n, a, lda, cdummy, u,
559  $ ldu, cdummy, -1, ierr )
560  lwrk_cunmqr = int( cdummy(1) )
561  ELSE
562  lwrk_cunmqr = 0
563  END IF
564  END IF
565  minwrk = 2
566  optwrk = 2
567  IF ( .NOT. (lsvec .OR. rsvec )) THEN
568 * .. minimal and optimal sizes of the complex workspace if
569 * only the singular values are requested
570  IF ( conda ) THEN
571  minwrk = max( n+lwqp3, lwcon, lwsvd )
572  ELSE
573  minwrk = max( n+lwqp3, lwsvd )
574  END IF
575  IF ( lquery ) THEN
576  CALL cgesvd( 'N', 'N', n, n, a, lda, s, u, ldu,
577  $ v, ldv, cdummy, -1, rdummy, ierr )
578  lwrk_cgesvd = int( cdummy(1) )
579  IF ( conda ) THEN
580  optwrk = max( n+lwrk_cgeqp3, n+lwcon, lwrk_cgesvd )
581  ELSE
582  optwrk = max( n+lwrk_cgeqp3, lwrk_cgesvd )
583  END IF
584  END IF
585  ELSE IF ( lsvec .AND. (.NOT.rsvec) ) THEN
586 * .. minimal and optimal sizes of the complex workspace if the
587 * singular values and the left singular vectors are requested
588  IF ( conda ) THEN
589  minwrk = n + max( lwqp3, lwcon, lwsvd, lwunq )
590  ELSE
591  minwrk = n + max( lwqp3, lwsvd, lwunq )
592  END IF
593  IF ( lquery ) THEN
594  IF ( rtrans ) THEN
595  CALL cgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
596  $ v, ldv, cdummy, -1, rdummy, ierr )
597  ELSE
598  CALL cgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
599  $ v, ldv, cdummy, -1, rdummy, ierr )
600  END IF
601  lwrk_cgesvd = int( cdummy(1) )
602  IF ( conda ) THEN
603  optwrk = n + max( lwrk_cgeqp3, lwcon, lwrk_cgesvd,
604  $ lwrk_cunmqr )
605  ELSE
606  optwrk = n + max( lwrk_cgeqp3, lwrk_cgesvd,
607  $ lwrk_cunmqr )
608  END IF
609  END IF
610  ELSE IF ( rsvec .AND. (.NOT.lsvec) ) THEN
611 * .. minimal and optimal sizes of the complex workspace if the
612 * singular values and the right singular vectors are requested
613  IF ( conda ) THEN
614  minwrk = n + max( lwqp3, lwcon, lwsvd )
615  ELSE
616  minwrk = n + max( lwqp3, lwsvd )
617  END IF
618  IF ( lquery ) THEN
619  IF ( rtrans ) THEN
620  CALL cgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
621  $ v, ldv, cdummy, -1, rdummy, ierr )
622  ELSE
623  CALL cgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
624  $ v, ldv, cdummy, -1, rdummy, ierr )
625  END IF
626  lwrk_cgesvd = int( cdummy(1) )
627  IF ( conda ) THEN
628  optwrk = n + max( lwrk_cgeqp3, lwcon, lwrk_cgesvd )
629  ELSE
630  optwrk = n + max( lwrk_cgeqp3, lwrk_cgesvd )
631  END IF
632  END IF
633  ELSE
634 * .. minimal and optimal sizes of the complex workspace if the
635 * full SVD is requested
636  IF ( rtrans ) THEN
637  minwrk = max( lwqp3, lwsvd, lwunq )
638  IF ( conda ) minwrk = max( minwrk, lwcon )
639  minwrk = minwrk + n
640  IF ( wntva ) THEN
641 * .. minimal workspace length for N x N/2 CGEQRF
642  lwqrf = max( n/2, 1 )
643 * .. minimal workspace lengt for N/2 x N/2 CGESVD
644  lwsvd2 = max( 3 * (n/2), 1 )
645  lwunq2 = max( n, 1 )
646  minwrk2 = max( lwqp3, n/2+lwqrf, n/2+lwsvd2,
647  $ n/2+lwunq2, lwunq )
648  IF ( conda ) minwrk2 = max( minwrk2, lwcon )
649  minwrk2 = n + minwrk2
650  minwrk = max( minwrk, minwrk2 )
651  END IF
652  ELSE
653  minwrk = max( lwqp3, lwsvd, lwunq )
654  IF ( conda ) minwrk = max( minwrk, lwcon )
655  minwrk = minwrk + n
656  IF ( wntva ) THEN
657 * .. minimal workspace length for N/2 x N CGELQF
658  lwlqf = max( n/2, 1 )
659  lwsvd2 = max( 3 * (n/2), 1 )
660  lwunlq = max( n , 1 )
661  minwrk2 = max( lwqp3, n/2+lwlqf, n/2+lwsvd2,
662  $ n/2+lwunlq, lwunq )
663  IF ( conda ) minwrk2 = max( minwrk2, lwcon )
664  minwrk2 = n + minwrk2
665  minwrk = max( minwrk, minwrk2 )
666  END IF
667  END IF
668  IF ( lquery ) THEN
669  IF ( rtrans ) THEN
670  CALL cgesvd( 'O', 'A', n, n, a, lda, s, u, ldu,
671  $ v, ldv, cdummy, -1, rdummy, ierr )
672  lwrk_cgesvd = int( cdummy(1) )
673  optwrk = max(lwrk_cgeqp3,lwrk_cgesvd,lwrk_cunmqr)
674  IF ( conda ) optwrk = max( optwrk, lwcon )
675  optwrk = n + optwrk
676  IF ( wntva ) THEN
677  CALL cgeqrf(n,n/2,u,ldu,cdummy,cdummy,-1,ierr)
678  lwrk_cgeqrf = int( cdummy(1) )
679  CALL cgesvd( 'S', 'O', n/2,n/2, v,ldv, s, u,ldu,
680  $ v, ldv, cdummy, -1, rdummy, ierr )
681  lwrk_cgesvd2 = int( cdummy(1) )
682  CALL cunmqr( 'R', 'C', n, n, n/2, u, ldu, cdummy,
683  $ v, ldv, cdummy, -1, ierr )
684  lwrk_cunmqr2 = int( cdummy(1) )
685  optwrk2 = max( lwrk_cgeqp3, n/2+lwrk_cgeqrf,
686  $ n/2+lwrk_cgesvd2, n/2+lwrk_cunmqr2 )
687  IF ( conda ) optwrk2 = max( optwrk2, lwcon )
688  optwrk2 = n + optwrk2
689  optwrk = max( optwrk, optwrk2 )
690  END IF
691  ELSE
692  CALL cgesvd( 'S', 'O', n, n, a, lda, s, u, ldu,
693  $ v, ldv, cdummy, -1, rdummy, ierr )
694  lwrk_cgesvd = int( cdummy(1) )
695  optwrk = max(lwrk_cgeqp3,lwrk_cgesvd,lwrk_cunmqr)
696  IF ( conda ) optwrk = max( optwrk, lwcon )
697  optwrk = n + optwrk
698  IF ( wntva ) THEN
699  CALL cgelqf(n/2,n,u,ldu,cdummy,cdummy,-1,ierr)
700  lwrk_cgelqf = int( cdummy(1) )
701  CALL cgesvd( 'S','O', n/2,n/2, v, ldv, s, u, ldu,
702  $ v, ldv, cdummy, -1, rdummy, ierr )
703  lwrk_cgesvd2 = int( cdummy(1) )
704  CALL cunmlq( 'R', 'N', n, n, n/2, u, ldu, cdummy,
705  $ v, ldv, cdummy,-1,ierr )
706  lwrk_cunmlq = int( cdummy(1) )
707  optwrk2 = max( lwrk_cgeqp3, n/2+lwrk_cgelqf,
708  $ n/2+lwrk_cgesvd2, n/2+lwrk_cunmlq )
709  IF ( conda ) optwrk2 = max( optwrk2, lwcon )
710  optwrk2 = n + optwrk2
711  optwrk = max( optwrk, optwrk2 )
712  END IF
713  END IF
714  END IF
715  END IF
716 *
717  minwrk = max( 2, minwrk )
718  optwrk = max( 2, optwrk )
719  IF ( lcwork .LT. minwrk .AND. (.NOT.lquery) ) info = -19
720 *
721  END IF
722 *
723  IF (info .EQ. 0 .AND. lrwork .LT. rminwrk .AND. .NOT. lquery) THEN
724  info = -21
725  END IF
726  IF( info.NE.0 ) THEN
727  CALL xerbla( 'CGESVDQ', -info )
728  RETURN
729  ELSE IF ( lquery ) THEN
730 *
731 * Return optimal workspace
732 *
733  iwork(1) = iminwrk
734  cwork(1) = optwrk
735  cwork(2) = minwrk
736  rwork(1) = rminwrk
737  RETURN
738  END IF
739 *
740 * Quick return if the matrix is void.
741 *
742  IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) ) THEN
743 * .. all output is void.
744  RETURN
745  END IF
746 *
747  big = slamch('O')
748  ascaled = .false.
749  IF ( rowprm ) THEN
750 * .. reordering the rows in decreasing sequence in the
751 * ell-infinity norm - this enhances numerical robustness in
752 * the case of differently scaled rows.
753  DO 1904 p = 1, m
754 * RWORK(p) = ABS( A(p,ICAMAX(N,A(p,1),LDA)) )
755 * [[CLANGE will return NaN if an entry of the p-th row is Nan]]
756  rwork(p) = clange( 'M', 1, n, a(p,1), lda, rdummy )
757 * .. check for NaN's and Inf's
758  IF ( ( rwork(p) .NE. rwork(p) ) .OR.
759  $ ( (rwork(p)*zero) .NE. zero ) ) THEN
760  info = - 8
761  CALL xerbla( 'CGESVDQ', -info )
762  RETURN
763  END IF
764  1904 CONTINUE
765  DO 1952 p = 1, m - 1
766  q = isamax( m-p+1, rwork(p), 1 ) + p - 1
767  iwork(n+p) = q
768  IF ( p .NE. q ) THEN
769  rtmp = rwork(p)
770  rwork(p) = rwork(q)
771  rwork(q) = rtmp
772  END IF
773  1952 CONTINUE
774 *
775  IF ( rwork(1) .EQ. zero ) THEN
776 * Quick return: A is the M x N zero matrix.
777  numrank = 0
778  CALL slaset( 'G', n, 1, zero, zero, s, n )
779  IF ( wntus ) CALL claset('G', m, n, czero, cone, u, ldu)
780  IF ( wntua ) CALL claset('G', m, m, czero, cone, u, ldu)
781  IF ( wntva ) CALL claset('G', n, n, czero, cone, v, ldv)
782  IF ( wntuf ) THEN
783  CALL claset( 'G', n, 1, czero, czero, cwork, n )
784  CALL claset( 'G', m, n, czero, cone, u, ldu )
785  END IF
786  DO 5001 p = 1, n
787  iwork(p) = p
788  5001 CONTINUE
789  IF ( rowprm ) THEN
790  DO 5002 p = n + 1, n + m - 1
791  iwork(p) = p - n
792  5002 CONTINUE
793  END IF
794  IF ( conda ) rwork(1) = -1
795  rwork(2) = -1
796  RETURN
797  END IF
798 *
799  IF ( rwork(1) .GT. big / sqrt(real(m)) ) THEN
800 * .. to prevent overflow in the QR factorization, scale the
801 * matrix by 1/sqrt(M) if too large entry detected
802  CALL clascl('G',0,0,sqrt(real(m)),one, m,n, a,lda, ierr)
803  ascaled = .true.
804  END IF
805  CALL claswp( n, a, lda, 1, m-1, iwork(n+1), 1 )
806  END IF
807 *
808 * .. At this stage, preemptive scaling is done only to avoid column
809 * norms overflows during the QR factorization. The SVD procedure should
810 * have its own scaling to save the singular values from overflows and
811 * underflows. That depends on the SVD procedure.
812 *
813  IF ( .NOT.rowprm ) THEN
814  rtmp = clange( 'M', m, n, a, lda, rwork )
815  IF ( ( rtmp .NE. rtmp ) .OR.
816  $ ( (rtmp*zero) .NE. zero ) ) THEN
817  info = - 8
818  CALL xerbla( 'CGESVDQ', -info )
819  RETURN
820  END IF
821  IF ( rtmp .GT. big / sqrt(real(m)) ) THEN
822 * .. to prevent overflow in the QR factorization, scale the
823 * matrix by 1/sqrt(M) if too large entry detected
824  CALL clascl('G',0,0, sqrt(real(m)),one, m,n, a,lda, ierr)
825  ascaled = .true.
826  END IF
827  END IF
828 *
829 * .. QR factorization with column pivoting
830 *
831 * A * P = Q * [ R ]
832 * [ 0 ]
833 *
834  DO 1963 p = 1, n
835 * .. all columns are free columns
836  iwork(p) = 0
837  1963 CONTINUE
838  CALL cgeqp3( m, n, a, lda, iwork, cwork, cwork(n+1), lcwork-n,
839  $ rwork, ierr )
840 *
841 * If the user requested accuracy level allows truncation in the
842 * computed upper triangular factor, the matrix R is examined and,
843 * if possible, replaced with its leading upper trapezoidal part.
844 *
845  epsln = slamch('E')
846  sfmin = slamch('S')
847 * SMALL = SFMIN / EPSLN
848  nr = n
849 *
850  IF ( accla ) THEN
851 *
852 * Standard absolute error bound suffices. All sigma_i with
853 * sigma_i < N*EPS*||A||_F are flushed to zero. This is an
854 * aggressive enforcement of lower numerical rank by introducing a
855 * backward error of the order of N*EPS*||A||_F.
856  nr = 1
857  rtmp = sqrt(real(n))*epsln
858  DO 3001 p = 2, n
859  IF ( abs(a(p,p)) .LT. (rtmp*abs(a(1,1))) ) GO TO 3002
860  nr = nr + 1
861  3001 CONTINUE
862  3002 CONTINUE
863 *
864  ELSEIF ( acclm ) THEN
865 * .. similarly as above, only slightly more gentle (less aggressive).
866 * Sudden drop on the diagonal of R is used as the criterion for being
867 * close-to-rank-deficient. The threshold is set to EPSLN=SLAMCH('E').
868 * [[This can be made more flexible by replacing this hard-coded value
869 * with a user specified threshold.]] Also, the values that underflow
870 * will be truncated.
871  nr = 1
872  DO 3401 p = 2, n
873  IF ( ( abs(a(p,p)) .LT. (epsln*abs(a(p-1,p-1))) ) .OR.
874  $ ( abs(a(p,p)) .LT. sfmin ) ) GO TO 3402
875  nr = nr + 1
876  3401 CONTINUE
877  3402 CONTINUE
878 *
879  ELSE
880 * .. RRQR not authorized to determine numerical rank except in the
881 * obvious case of zero pivots.
882 * .. inspect R for exact zeros on the diagonal;
883 * R(i,i)=0 => R(i:N,i:N)=0.
884  nr = 1
885  DO 3501 p = 2, n
886  IF ( abs(a(p,p)) .EQ. zero ) GO TO 3502
887  nr = nr + 1
888  3501 CONTINUE
889  3502 CONTINUE
890 *
891  IF ( conda ) THEN
892 * Estimate the scaled condition number of A. Use the fact that it is
893 * the same as the scaled condition number of R.
894 * .. V is used as workspace
895  CALL clacpy( 'U', n, n, a, lda, v, ldv )
896 * Only the leading NR x NR submatrix of the triangular factor
897 * is considered. Only if NR=N will this give a reliable error
898 * bound. However, even for NR < N, this can be used on an
899 * expert level and obtain useful information in the sense of
900 * perturbation theory.
901  DO 3053 p = 1, nr
902  rtmp = scnrm2( p, v(1,p), 1 )
903  CALL csscal( p, one/rtmp, v(1,p), 1 )
904  3053 CONTINUE
905  IF ( .NOT. ( lsvec .OR. rsvec ) ) THEN
906  CALL cpocon( 'U', nr, v, ldv, one, rtmp,
907  $ cwork, rwork, ierr )
908  ELSE
909  CALL cpocon( 'U', nr, v, ldv, one, rtmp,
910  $ cwork(n+1), rwork, ierr )
911  END IF
912  sconda = one / sqrt(rtmp)
913 * For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1),
914 * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
915 * See the reference [1] for more details.
916  END IF
917 *
918  ENDIF
919 *
920  IF ( wntur ) THEN
921  n1 = nr
922  ELSE IF ( wntus .OR. wntuf) THEN
923  n1 = n
924  ELSE IF ( wntua ) THEN
925  n1 = m
926  END IF
927 *
928  IF ( .NOT. ( rsvec .OR. lsvec ) ) THEN
929 *.......................................................................
930 * .. only the singular values are requested
931 *.......................................................................
932  IF ( rtrans ) THEN
933 *
934 * .. compute the singular values of R**H = [A](1:NR,1:N)**H
935 * .. set the lower triangle of [A] to [A](1:NR,1:N)**H and
936 * the upper triangle of [A] to zero.
937  DO 1146 p = 1, min( n, nr )
938  a(p,p) = conjg(a(p,p))
939  DO 1147 q = p + 1, n
940  a(q,p) = conjg(a(p,q))
941  IF ( q .LE. nr ) a(p,q) = czero
942  1147 CONTINUE
943  1146 CONTINUE
944 *
945  CALL cgesvd( 'N', 'N', n, nr, a, lda, s, u, ldu,
946  $ v, ldv, cwork, lcwork, rwork, info )
947 *
948  ELSE
949 *
950 * .. compute the singular values of R = [A](1:NR,1:N)
951 *
952  IF ( nr .GT. 1 )
953  $ CALL claset( 'L', nr-1,nr-1, czero,czero, a(2,1), lda )
954  CALL cgesvd( 'N', 'N', nr, n, a, lda, s, u, ldu,
955  $ v, ldv, cwork, lcwork, rwork, info )
956 *
957  END IF
958 *
959  ELSE IF ( lsvec .AND. ( .NOT. rsvec) ) THEN
960 *.......................................................................
961 * .. the singular values and the left singular vectors requested
962 *.......................................................................""""""""
963  IF ( rtrans ) THEN
964 * .. apply CGESVD to R**H
965 * .. copy R**H into [U] and overwrite [U] with the right singular
966 * vectors of R
967  DO 1192 p = 1, nr
968  DO 1193 q = p, n
969  u(q,p) = conjg(a(p,q))
970  1193 CONTINUE
971  1192 CONTINUE
972  IF ( nr .GT. 1 )
973  $ CALL claset( 'U', nr-1,nr-1, czero,czero, u(1,2), ldu )
974 * .. the left singular vectors not computed, the NR right singular
975 * vectors overwrite [U](1:NR,1:NR) as conjugate transposed. These
976 * will be pre-multiplied by Q to build the left singular vectors of A.
977  CALL cgesvd( 'N', 'O', n, nr, u, ldu, s, u, ldu,
978  $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
979 *
980  DO 1119 p = 1, nr
981  u(p,p) = conjg(u(p,p))
982  DO 1120 q = p + 1, nr
983  ctmp = conjg(u(q,p))
984  u(q,p) = conjg(u(p,q))
985  u(p,q) = ctmp
986  1120 CONTINUE
987  1119 CONTINUE
988 *
989  ELSE
990 * .. apply CGESVD to R
991 * .. copy R into [U] and overwrite [U] with the left singular vectors
992  CALL clacpy( 'U', nr, n, a, lda, u, ldu )
993  IF ( nr .GT. 1 )
994  $ CALL claset( 'L', nr-1, nr-1, czero, czero, u(2,1), ldu )
995 * .. the right singular vectors not computed, the NR left singular
996 * vectors overwrite [U](1:NR,1:NR)
997  CALL cgesvd( 'O', 'N', nr, n, u, ldu, s, u, ldu,
998  $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
999 * .. now [U](1:NR,1:NR) contains the NR left singular vectors of
1000 * R. These will be pre-multiplied by Q to build the left singular
1001 * vectors of A.
1002  END IF
1003 *
1004 * .. assemble the left singular vector matrix U of dimensions
1005 * (M x NR) or (M x N) or (M x M).
1006  IF ( ( nr .LT. m ) .AND. ( .NOT.wntuf ) ) THEN
1007  CALL claset('A', m-nr, nr, czero, czero, u(nr+1,1), ldu)
1008  IF ( nr .LT. n1 ) THEN
1009  CALL claset( 'A',nr,n1-nr,czero,czero,u(1,nr+1), ldu )
1010  CALL claset( 'A',m-nr,n1-nr,czero,cone,
1011  $ u(nr+1,nr+1), ldu )
1012  END IF
1013  END IF
1014 *
1015 * The Q matrix from the first QRF is built into the left singular
1016 * vectors matrix U.
1017 *
1018  IF ( .NOT.wntuf )
1019  $ CALL cunmqr( 'L', 'N', m, n1, n, a, lda, cwork, u,
1020  $ ldu, cwork(n+1), lcwork-n, ierr )
1021  IF ( rowprm .AND. .NOT.wntuf )
1022  $ CALL claswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
1023 *
1024  ELSE IF ( rsvec .AND. ( .NOT. lsvec ) ) THEN
1025 *.......................................................................
1026 * .. the singular values and the right singular vectors requested
1027 *.......................................................................
1028  IF ( rtrans ) THEN
1029 * .. apply CGESVD to R**H
1030 * .. copy R**H into V and overwrite V with the left singular vectors
1031  DO 1165 p = 1, nr
1032  DO 1166 q = p, n
1033  v(q,p) = conjg(a(p,q))
1034  1166 CONTINUE
1035  1165 CONTINUE
1036  IF ( nr .GT. 1 )
1037  $ CALL claset( 'U', nr-1,nr-1, czero,czero, v(1,2), ldv )
1038 * .. the left singular vectors of R**H overwrite V, the right singular
1039 * vectors not computed
1040  IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1041  CALL cgesvd( 'O', 'N', n, nr, v, ldv, s, u, ldu,
1042  $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1043 *
1044  DO 1121 p = 1, nr
1045  v(p,p) = conjg(v(p,p))
1046  DO 1122 q = p + 1, nr
1047  ctmp = conjg(v(q,p))
1048  v(q,p) = conjg(v(p,q))
1049  v(p,q) = ctmp
1050  1122 CONTINUE
1051  1121 CONTINUE
1052 *
1053  IF ( nr .LT. n ) THEN
1054  DO 1103 p = 1, nr
1055  DO 1104 q = nr + 1, n
1056  v(p,q) = conjg(v(q,p))
1057  1104 CONTINUE
1058  1103 CONTINUE
1059  END IF
1060  CALL clapmt( .false., nr, n, v, ldv, iwork )
1061  ELSE
1062 * .. need all N right singular vectors and NR < N
1063 * [!] This is simple implementation that augments [V](1:N,1:NR)
1064 * by padding a zero block. In the case NR << N, a more efficient
1065 * way is to first use the QR factorization. For more details
1066 * how to implement this, see the " FULL SVD " branch.
1067  CALL claset('G', n, n-nr, czero, czero, v(1,nr+1), ldv)
1068  CALL cgesvd( 'O', 'N', n, n, v, ldv, s, u, ldu,
1069  $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1070 *
1071  DO 1123 p = 1, n
1072  v(p,p) = conjg(v(p,p))
1073  DO 1124 q = p + 1, n
1074  ctmp = conjg(v(q,p))
1075  v(q,p) = conjg(v(p,q))
1076  v(p,q) = ctmp
1077  1124 CONTINUE
1078  1123 CONTINUE
1079  CALL clapmt( .false., n, n, v, ldv, iwork )
1080  END IF
1081 *
1082  ELSE
1083 * .. aply CGESVD to R
1084 * .. copy R into V and overwrite V with the right singular vectors
1085  CALL clacpy( 'U', nr, n, a, lda, v, ldv )
1086  IF ( nr .GT. 1 )
1087  $ CALL claset( 'L', nr-1, nr-1, czero, czero, 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 cgesvd( 'N', 'O', nr, n, v, ldv, s, u, ldu,
1092  $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1093  CALL clapmt( .false., nr, n, v, ldv, iwork )
1094 * .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H
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 claset('G', n-nr, n, czero,czero, v(nr+1,1), ldv)
1102  CALL cgesvd( 'N', 'O', n, n, v, ldv, s, u, ldu,
1103  $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1104  CALL clapmt( .false., n, n, v, ldv, iwork )
1105  END IF
1106 * .. now [V] contains the adjoint of the 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 CGESVD to R**H [[this option is left for R&D&T]]
1117 *
1118  IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1119 * .. copy R**H into [V] and overwrite [V] with the left singular
1120 * vectors of R**H
1121  DO 1168 p = 1, nr
1122  DO 1169 q = p, n
1123  v(q,p) = conjg(a(p,q))
1124  1169 CONTINUE
1125  1168 CONTINUE
1126  IF ( nr .GT. 1 )
1127  $ CALL claset( 'U', nr-1,nr-1, czero,czero, v(1,2), ldv )
1128 *
1129 * .. the left singular vectors of R**H overwrite [V], the NR right
1130 * singular vectors of R**H stored in [U](1:NR,1:NR) as conjugate
1131 * transposed
1132  CALL cgesvd( 'O', 'A', n, nr, v, ldv, s, v, ldv,
1133  $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1134 * .. assemble V
1135  DO 1115 p = 1, nr
1136  v(p,p) = conjg(v(p,p))
1137  DO 1116 q = p + 1, nr
1138  ctmp = conjg(v(q,p))
1139  v(q,p) = conjg(v(p,q))
1140  v(p,q) = ctmp
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) = conjg(v(q,p))
1147  1102 CONTINUE
1148  1101 CONTINUE
1149  END IF
1150  CALL clapmt( .false., nr, n, v, ldv, iwork )
1151 *
1152  DO 1117 p = 1, nr
1153  u(p,p) = conjg(u(p,p))
1154  DO 1118 q = p + 1, nr
1155  ctmp = conjg(u(q,p))
1156  u(q,p) = conjg(u(p,q))
1157  u(p,q) = ctmp
1158  1118 CONTINUE
1159  1117 CONTINUE
1160 *
1161  IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1162  CALL claset('A', m-nr,nr, czero,czero, u(nr+1,1), ldu)
1163  IF ( nr .LT. n1 ) THEN
1164  CALL claset('A',nr,n1-nr,czero,czero,u(1,nr+1),ldu)
1165  CALL claset( 'A',m-nr,n1-nr,czero,cone,
1166  $ u(nr+1,nr+1), ldu )
1167  END IF
1168  END IF
1169 *
1170  ELSE
1171 * .. need all N right singular vectors and NR < N
1172 * .. copy R**H into [V] and overwrite [V] with the left singular
1173 * vectors of R**H
1174 * [[The optimal ratio N/NR for using QRF instead of padding
1175 * with zeros. Here hard coded to 2; it must be at least
1176 * two due to work space constraints.]]
1177 * OPTRATIO = ILAENV(6, 'CGESVD', 'S' // 'O', NR,N,0,0)
1178 * OPTRATIO = MAX( OPTRATIO, 2 )
1179  optratio = 2
1180  IF ( optratio*nr .GT. n ) THEN
1181  DO 1198 p = 1, nr
1182  DO 1199 q = p, n
1183  v(q,p) = conjg(a(p,q))
1184  1199 CONTINUE
1185  1198 CONTINUE
1186  IF ( nr .GT. 1 )
1187  $ CALL claset('U',nr-1,nr-1, czero,czero, v(1,2),ldv)
1188 *
1189  CALL claset('A',n,n-nr,czero,czero,v(1,nr+1),ldv)
1190  CALL cgesvd( 'O', 'A', n, n, v, ldv, s, v, ldv,
1191  $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1192 *
1193  DO 1113 p = 1, n
1194  v(p,p) = conjg(v(p,p))
1195  DO 1114 q = p + 1, n
1196  ctmp = conjg(v(q,p))
1197  v(q,p) = conjg(v(p,q))
1198  v(p,q) = ctmp
1199  1114 CONTINUE
1200  1113 CONTINUE
1201  CALL clapmt( .false., n, n, v, ldv, iwork )
1202 * .. assemble the left singular vector matrix U of dimensions
1203 * (M x N1), i.e. (M x N) or (M x M).
1204 *
1205  DO 1111 p = 1, n
1206  u(p,p) = conjg(u(p,p))
1207  DO 1112 q = p + 1, n
1208  ctmp = conjg(u(q,p))
1209  u(q,p) = conjg(u(p,q))
1210  u(p,q) = ctmp
1211  1112 CONTINUE
1212  1111 CONTINUE
1213 *
1214  IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) THEN
1215  CALL claset('A',m-n,n,czero,czero,u(n+1,1),ldu)
1216  IF ( n .LT. n1 ) THEN
1217  CALL claset('A',n,n1-n,czero,czero,u(1,n+1),ldu)
1218  CALL claset('A',m-n,n1-n,czero,cone,
1219  $ u(n+1,n+1), ldu )
1220  END IF
1221  END IF
1222  ELSE
1223 * .. copy R**H into [U] and overwrite [U] with the right
1224 * singular vectors of R
1225  DO 1196 p = 1, nr
1226  DO 1197 q = p, n
1227  u(q,nr+p) = conjg(a(p,q))
1228  1197 CONTINUE
1229  1196 CONTINUE
1230  IF ( nr .GT. 1 )
1231  $ CALL claset('U',nr-1,nr-1,czero,czero,u(1,nr+2),ldu)
1232  CALL cgeqrf( n, nr, u(1,nr+1), ldu, cwork(n+1),
1233  $ cwork(n+nr+1), lcwork-n-nr, ierr )
1234  DO 1143 p = 1, nr
1235  DO 1144 q = 1, n
1236  v(q,p) = conjg(u(p,nr+q))
1237  1144 CONTINUE
1238  1143 CONTINUE
1239  CALL claset('U',nr-1,nr-1,czero,czero,v(1,2),ldv)
1240  CALL cgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
1241  $ v,ldv, cwork(n+nr+1),lcwork-n-nr,rwork, info )
1242  CALL claset('A',n-nr,nr,czero,czero,v(nr+1,1),ldv)
1243  CALL claset('A',nr,n-nr,czero,czero,v(1,nr+1),ldv)
1244  CALL claset('A',n-nr,n-nr,czero,cone,v(nr+1,nr+1),ldv)
1245  CALL cunmqr('R','C', n, n, nr, u(1,nr+1), ldu,
1246  $ cwork(n+1),v,ldv,cwork(n+nr+1),lcwork-n-nr,ierr)
1247  CALL clapmt( .false., n, n, v, ldv, iwork )
1248 * .. assemble the left singular vector matrix U of dimensions
1249 * (M x NR) or (M x N) or (M x M).
1250  IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1251  CALL claset('A',m-nr,nr,czero,czero,u(nr+1,1),ldu)
1252  IF ( nr .LT. n1 ) THEN
1253  CALL claset('A',nr,n1-nr,czero,czero,u(1,nr+1),ldu)
1254  CALL claset( 'A',m-nr,n1-nr,czero,cone,
1255  $ u(nr+1,nr+1),ldu)
1256  END IF
1257  END IF
1258  END IF
1259  END IF
1260 *
1261  ELSE
1262 *
1263 * .. apply CGESVD to R [[this is the recommended option]]
1264 *
1265  IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1266 * .. copy R into [V] and overwrite V with the right singular vectors
1267  CALL clacpy( 'U', nr, n, a, lda, v, ldv )
1268  IF ( nr .GT. 1 )
1269  $ CALL claset( 'L', nr-1,nr-1, czero,czero, v(2,1), ldv )
1270 * .. the right singular vectors of R overwrite [V], the NR left
1271 * singular vectors of R stored in [U](1:NR,1:NR)
1272  CALL cgesvd( 'S', 'O', nr, n, v, ldv, s, u, ldu,
1273  $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1274  CALL clapmt( .false., nr, n, v, ldv, iwork )
1275 * .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H
1276 * .. assemble the left singular vector matrix U of dimensions
1277 * (M x NR) or (M x N) or (M x M).
1278  IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1279  CALL claset('A', m-nr,nr, czero,czero, u(nr+1,1), ldu)
1280  IF ( nr .LT. n1 ) THEN
1281  CALL claset('A',nr,n1-nr,czero,czero,u(1,nr+1),ldu)
1282  CALL claset( 'A',m-nr,n1-nr,czero,cone,
1283  $ u(nr+1,nr+1), ldu )
1284  END IF
1285  END IF
1286 *
1287  ELSE
1288 * .. need all N right singular vectors and NR < N
1289 * .. the requested number of the left singular vectors
1290 * is then N1 (N or M)
1291 * [[The optimal ratio N/NR for using LQ instead of padding
1292 * with zeros. Here hard coded to 2; it must be at least
1293 * two due to work space constraints.]]
1294 * OPTRATIO = ILAENV(6, 'CGESVD', 'S' // 'O', NR,N,0,0)
1295 * OPTRATIO = MAX( OPTRATIO, 2 )
1296  optratio = 2
1297  IF ( optratio * nr .GT. n ) THEN
1298  CALL clacpy( 'U', nr, n, a, lda, v, ldv )
1299  IF ( nr .GT. 1 )
1300  $ CALL claset('L', nr-1,nr-1, czero,czero, v(2,1),ldv)
1301 * .. the right singular vectors of R overwrite [V], the NR left
1302 * singular vectors of R stored in [U](1:NR,1:NR)
1303  CALL claset('A', n-nr,n, czero,czero, v(nr+1,1),ldv)
1304  CALL cgesvd( 'S', 'O', n, n, v, ldv, s, u, ldu,
1305  $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1306  CALL clapmt( .false., n, n, v, ldv, iwork )
1307 * .. now [V] contains the adjoint of the matrix of the right
1308 * singular vectors of A. The leading N left singular vectors
1309 * are in [U](1:N,1:N)
1310 * .. assemble the left singular vector matrix U of dimensions
1311 * (M x N1), i.e. (M x N) or (M x M).
1312  IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) THEN
1313  CALL claset('A',m-n,n,czero,czero,u(n+1,1),ldu)
1314  IF ( n .LT. n1 ) THEN
1315  CALL claset('A',n,n1-n,czero,czero,u(1,n+1),ldu)
1316  CALL claset( 'A',m-n,n1-n,czero,cone,
1317  $ u(n+1,n+1), ldu )
1318  END IF
1319  END IF
1320  ELSE
1321  CALL clacpy( 'U', nr, n, a, lda, u(nr+1,1), ldu )
1322  IF ( nr .GT. 1 )
1323  $ CALL claset('L',nr-1,nr-1,czero,czero,u(nr+2,1),ldu)
1324  CALL cgelqf( nr, n, u(nr+1,1), ldu, cwork(n+1),
1325  $ cwork(n+nr+1), lcwork-n-nr, ierr )
1326  CALL clacpy('L',nr,nr,u(nr+1,1),ldu,v,ldv)
1327  IF ( nr .GT. 1 )
1328  $ CALL claset('U',nr-1,nr-1,czero,czero,v(1,2),ldv)
1329  CALL cgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
1330  $ v, ldv, cwork(n+nr+1), lcwork-n-nr, rwork, info )
1331  CALL claset('A',n-nr,nr,czero,czero,v(nr+1,1),ldv)
1332  CALL claset('A',nr,n-nr,czero,czero,v(1,nr+1),ldv)
1333  CALL claset('A',n-nr,n-nr,czero,cone,v(nr+1,nr+1),ldv)
1334  CALL cunmlq('R','N',n,n,nr,u(nr+1,1),ldu,cwork(n+1),
1335  $ v, ldv, cwork(n+nr+1),lcwork-n-nr,ierr)
1336  CALL clapmt( .false., n, n, v, ldv, iwork )
1337 * .. assemble the left singular vector matrix U of dimensions
1338 * (M x NR) or (M x N) or (M x M).
1339  IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1340  CALL claset('A',m-nr,nr,czero,czero,u(nr+1,1),ldu)
1341  IF ( nr .LT. n1 ) THEN
1342  CALL claset('A',nr,n1-nr,czero,czero,u(1,nr+1),ldu)
1343  CALL claset( 'A',m-nr,n1-nr,czero,cone,
1344  $ u(nr+1,nr+1), ldu )
1345  END IF
1346  END IF
1347  END IF
1348  END IF
1349 * .. end of the "R**H or R" branch
1350  END IF
1351 *
1352 * The Q matrix from the first QRF is built into the left singular
1353 * vectors matrix U.
1354 *
1355  IF ( .NOT. wntuf )
1356  $ CALL cunmqr( 'L', 'N', m, n1, n, a, lda, cwork, u,
1357  $ ldu, cwork(n+1), lcwork-n, ierr )
1358  IF ( rowprm .AND. .NOT.wntuf )
1359  $ CALL claswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
1360 *
1361 * ... end of the "full SVD" branch
1362  END IF
1363 *
1364 * Check whether some singular values are returned as zeros, e.g.
1365 * due to underflow, and update the numerical rank.
1366  p = nr
1367  DO 4001 q = p, 1, -1
1368  IF ( s(q) .GT. zero ) GO TO 4002
1369  nr = nr - 1
1370  4001 CONTINUE
1371  4002 CONTINUE
1372 *
1373 * .. if numerical rank deficiency is detected, the truncated
1374 * singular values are set to zero.
1375  IF ( nr .LT. n ) CALL slaset( 'G', n-nr,1, zero,zero, s(nr+1), n )
1376 * .. undo scaling; this may cause overflow in the largest singular
1377 * values.
1378  IF ( ascaled )
1379  $ CALL slascl( 'G',0,0, one,sqrt(real(m)), nr,1, s, n, ierr )
1380  IF ( conda ) rwork(1) = sconda
1381  rwork(2) = p - nr
1382 * .. p-NR is the number of singular values that are computed as
1383 * exact zeros in CGESVD() applied to the (possibly truncated)
1384 * full row rank triangular (trapezoidal) factor of A.
1385  numrank = nr
1386 *
1387  RETURN
1388 *
1389 * End of CGESVDQ
1390 *
1391  END
csscal
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:80
cunmqr
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:170
cunmlq
subroutine cunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMLQ
Definition: cunmlq.f:170
cgeqp3
subroutine cgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, INFO)
CGEQP3
Definition: cgeqp3.f:161
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
claswp
subroutine claswp(N, A, LDA, K1, K2, IPIV, INCX)
CLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: claswp.f:117
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
clapmt
subroutine clapmt(FORWRD, M, N, X, LDX, K)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: clapmt.f:106
cgelqf
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
Definition: cgelqf.f:145
cgesvd
subroutine cgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
CGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: cgesvd.f:216
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
claset
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:108
cpocon
subroutine cpocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
CPOCON
Definition: cpocon.f:123
cgesvdq
subroutine cgesvdq(JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA, S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK, CWORK, LCWORK, RWORK, LRWORK, INFO)
CGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE...
Definition: cgesvdq.f:415
cgeqrf
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:147
clascl
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:145