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