LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
zgejsv.f
Go to the documentation of this file.
1 *> \brief \b ZGEJSV
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZGEJSV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgejsv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgejsv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgejsv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
22 * M, N, A, LDA, SVA, U, LDU, V, LDV,
23 * CWORK, LWORK, RWORK, LRWORK, IWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * IMPLICIT NONE
27 * INTEGER INFO, LDA, LDU, LDV, LWORK, M, N
28 * ..
29 * .. Array Arguments ..
30 * COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( LWORK )
31 * DOUBLE PRECISION SVA( N ), RWORK( LRWORK )
32 * INTEGER IWORK( * )
33 * CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> ZGEJSV computes the singular value decomposition (SVD) of a complex M-by-N
43 *> matrix [A], where M >= N. The SVD of [A] is written as
44 *>
45 *> [A] = [U] * [SIGMA] * [V]^*,
46 *>
47 *> where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N
48 *> diagonal elements, [U] is an M-by-N (or M-by-M) unitary matrix, and
49 *> [V] is an N-by-N unitary matrix. The diagonal elements of [SIGMA] are
50 *> the singular values of [A]. The columns of [U] and [V] are the left and
51 *> the right singular vectors of [A], respectively. The matrices [U] and [V]
52 *> are computed and stored in the arrays U and V, respectively. The diagonal
53 *> of [SIGMA] is computed and stored in the array SVA.
54 *> \endverbatim
55 *>
56 *> Arguments:
57 *> ==========
58 *>
59 *> \param[in] JOBA
60 *> \verbatim
61 *> JOBA is CHARACTER*1
62 *> Specifies the level of accuracy:
63 *> = 'C': This option works well (high relative accuracy) if A = B * D,
64 *> with well-conditioned B and arbitrary diagonal matrix D.
65 *> The accuracy cannot be spoiled by COLUMN scaling. The
66 *> accuracy of the computed output depends on the condition of
67 *> B, and the procedure aims at the best theoretical accuracy.
68 *> The relative error max_{i=1:N}|d sigma_i| / sigma_i is
69 *> bounded by f(M,N)*epsilon* cond(B), independent of D.
70 *> The input matrix is preprocessed with the QRF with column
71 *> pivoting. This initial preprocessing and preconditioning by
72 *> a rank revealing QR factorization is common for all values of
73 *> JOBA. Additional actions are specified as follows:
74 *> = 'E': Computation as with 'C' with an additional estimate of the
75 *> condition number of B. It provides a realistic error bound.
76 *> = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings
77 *> D1, D2, and well-conditioned matrix C, this option gives
78 *> higher accuracy than the 'C' option. If the structure of the
79 *> input matrix is not known, and relative accuracy is
80 *> desirable, then this option is advisable. The input matrix A
81 *> is preprocessed with QR factorization with FULL (row and
82 *> column) pivoting.
83 *> = 'G': Computation as with 'F' with an additional estimate of the
84 *> condition number of B, where A=B*D. If A has heavily weighted
85 *> rows, then using this condition number gives too pessimistic
86 *> error bound.
87 *> = 'A': Small singular values are not well determined by the data
88 *> and are considered as noisy; the matrix is treated as
89 *> numerically rank deficient. The error in the computed
90 *> singular values is bounded by f(m,n)*epsilon*||A||.
91 *> The computed SVD A = U * S * V^* restores A up to
92 *> f(m,n)*epsilon*||A||.
93 *> This gives the procedure the licence to discard (set to zero)
94 *> all singular values below N*epsilon*||A||.
95 *> = 'R': Similar as in 'A'. Rank revealing property of the initial
96 *> QR factorization is used do reveal (using triangular factor)
97 *> a gap sigma_{r+1} < epsilon * sigma_r in which case the
98 *> numerical RANK is declared to be r. The SVD is computed with
99 *> absolute error bounds, but more accurately than with 'A'.
100 *> \endverbatim
101 *>
102 *> \param[in] JOBU
103 *> \verbatim
104 *> JOBU is CHARACTER*1
105 *> Specifies whether to compute the columns of U:
106 *> = 'U': N columns of U are returned in the array U.
107 *> = 'F': full set of M left sing. vectors is returned in the array U.
108 *> = 'W': U may be used as workspace of length M*N. See the description
109 *> of U.
110 *> = 'N': U is not computed.
111 *> \endverbatim
112 *>
113 *> \param[in] JOBV
114 *> \verbatim
115 *> JOBV is CHARACTER*1
116 *> Specifies whether to compute the matrix V:
117 *> = 'V': N columns of V are returned in the array V; Jacobi rotations
118 *> are not explicitly accumulated.
119 *> = 'J': N columns of V are returned in the array V, but they are
120 *> computed as the product of Jacobi rotations, if JOBT = 'N'.
121 *> = 'W': V may be used as workspace of length N*N. See the description
122 *> of V.
123 *> = 'N': V is not computed.
124 *> \endverbatim
125 *>
126 *> \param[in] JOBR
127 *> \verbatim
128 *> JOBR is CHARACTER*1
129 *> Specifies the RANGE for the singular values. Issues the licence to
130 *> set to zero small positive singular values if they are outside
131 *> specified range. If A .NE. 0 is scaled so that the largest singular
132 *> value of c*A is around SQRT(BIG), BIG=DLAMCH('O'), then JOBR issues
133 *> the licence to kill columns of A whose norm in c*A is less than
134 *> SQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN,
135 *> where SFMIN=DLAMCH('S'), EPSLN=DLAMCH('E').
136 *> = 'N': Do not kill small columns of c*A. This option assumes that
137 *> BLAS and QR factorizations and triangular solvers are
138 *> implemented to work in that range. If the condition of A
139 *> is greater than BIG, use ZGESVJ.
140 *> = 'R': RESTRICTED range for sigma(c*A) is [SQRT(SFMIN), SQRT(BIG)]
141 *> (roughly, as described above). This option is recommended.
142 *> ===========================
143 *> For computing the singular values in the FULL range [SFMIN,BIG]
144 *> use ZGESVJ.
145 *> \endverbatim
146 *>
147 *> \param[in] JOBT
148 *> \verbatim
149 *> JOBT is CHARACTER*1
150 *> If the matrix is square then the procedure may determine to use
151 *> transposed A if A^* seems to be better with respect to convergence.
152 *> If the matrix is not square, JOBT is ignored.
153 *> The decision is based on two values of entropy over the adjoint
154 *> orbit of A^* * A. See the descriptions of WORK(6) and WORK(7).
155 *> = 'T': transpose if entropy test indicates possibly faster
156 *> convergence of Jacobi process if A^* is taken as input. If A is
157 *> replaced with A^*, then the row pivoting is included automatically.
158 *> = 'N': do not speculate.
159 *> The option 'T' can be used to compute only the singular values, or
160 *> the full SVD (U, SIGMA and V). For only one set of singular vectors
161 *> (U or V), the caller should provide both U and V, as one of the
162 *> matrices is used as workspace if the matrix A is transposed.
163 *> The implementer can easily remove this constraint and make the
164 *> code more complicated. See the descriptions of U and V.
165 *> In general, this option is considered experimental, and 'N'; should
166 *> be preferred. This is subject to changes in the future.
167 *> \endverbatim
168 *>
169 *> \param[in] JOBP
170 *> \verbatim
171 *> JOBP is CHARACTER*1
172 *> Issues the licence to introduce structured perturbations to drown
173 *> denormalized numbers. This licence should be active if the
174 *> denormals are poorly implemented, causing slow computation,
175 *> especially in cases of fast convergence (!). For details see [1,2].
176 *> For the sake of simplicity, this perturbations are included only
177 *> when the full SVD or only the singular values are requested. The
178 *> implementer/user can easily add the perturbation for the cases of
179 *> computing one set of singular vectors.
180 *> = 'P': introduce perturbation
181 *> = 'N': do not perturb
182 *> \endverbatim
183 *>
184 *> \param[in] M
185 *> \verbatim
186 *> M is INTEGER
187 *> The number of rows of the input matrix A. M >= 0.
188 *> \endverbatim
189 *>
190 *> \param[in] N
191 *> \verbatim
192 *> N is INTEGER
193 *> The number of columns of the input matrix A. M >= N >= 0.
194 *> \endverbatim
195 *>
196 *> \param[in,out] A
197 *> \verbatim
198 *> A is COMPLEX*16 array, dimension (LDA,N)
199 *> On entry, the M-by-N matrix A.
200 *> \endverbatim
201 *>
202 *> \param[in] LDA
203 *> \verbatim
204 *> LDA is INTEGER
205 *> The leading dimension of the array A. LDA >= max(1,M).
206 *> \endverbatim
207 *>
208 *> \param[out] SVA
209 *> \verbatim
210 *> SVA is DOUBLE PRECISION array, dimension (N)
211 *> On exit,
212 *> - For WORK(1)/WORK(2) = ONE: The singular values of A. During the
213 *> computation SVA contains Euclidean column norms of the
214 *> iterated matrices in the array A.
215 *> - For WORK(1) .NE. WORK(2): The singular values of A are
216 *> (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if
217 *> sigma_max(A) overflows or if small singular values have been
218 *> saved from underflow by scaling the input matrix A.
219 *> - If JOBR='R' then some of the singular values may be returned
220 *> as exact zeros obtained by "set to zero" because they are
221 *> below the numerical rank threshold or are denormalized numbers.
222 *> \endverbatim
223 *>
224 *> \param[out] U
225 *> \verbatim
226 *> U is COMPLEX*16 array, dimension ( LDU, N )
227 *> If JOBU = 'U', then U contains on exit the M-by-N matrix of
228 *> the left singular vectors.
229 *> If JOBU = 'F', then U contains on exit the M-by-M matrix of
230 *> the left singular vectors, including an ONB
231 *> of the orthogonal complement of the Range(A).
232 *> If JOBU = 'W' .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N),
233 *> then U is used as workspace if the procedure
234 *> replaces A with A^*. In that case, [V] is computed
235 *> in U as left singular vectors of A^* and then
236 *> copied back to the V array. This 'W' option is just
237 *> a reminder to the caller that in this case U is
238 *> reserved as workspace of length N*N.
239 *> If JOBU = 'N' U is not referenced, unless JOBT='T'.
240 *> \endverbatim
241 *>
242 *> \param[in] LDU
243 *> \verbatim
244 *> LDU is INTEGER
245 *> The leading dimension of the array U, LDU >= 1.
246 *> IF JOBU = 'U' or 'F' or 'W', then LDU >= M.
247 *> \endverbatim
248 *>
249 *> \param[out] V
250 *> \verbatim
251 *> V is COMPLEX*16 array, dimension ( LDV, N )
252 *> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
253 *> the right singular vectors;
254 *> If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N),
255 *> then V is used as workspace if the pprocedure
256 *> replaces A with A^*. In that case, [U] is computed
257 *> in V as right singular vectors of A^* and then
258 *> copied back to the U array. This 'W' option is just
259 *> a reminder to the caller that in this case V is
260 *> reserved as workspace of length N*N.
261 *> If JOBV = 'N' V is not referenced, unless JOBT='T'.
262 *> \endverbatim
263 *>
264 *> \param[in] LDV
265 *> \verbatim
266 *> LDV is INTEGER
267 *> The leading dimension of the array V, LDV >= 1.
268 *> If JOBV = 'V' or 'J' or 'W', then LDV >= N.
269 *> \endverbatim
270 *>
271 *> \param[out] CWORK
272 *> \verbatim
273 *> CWORK is COMPLEX*16 array, dimension (MAX(2,LWORK))
274 *> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or
275 *> LRWORK=-1), then on exit CWORK(1) contains the required length of
276 *> CWORK for the job parameters used in the call.
277 *> \endverbatim
278 *>
279 *> \param[in] LWORK
280 *> \verbatim
281 *> LWORK is INTEGER
282 *> Length of CWORK to confirm proper allocation of workspace.
283 *> LWORK depends on the job:
284 *>
285 *> 1. If only SIGMA is needed ( JOBU = 'N', JOBV = 'N' ) and
286 *> 1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'):
287 *> LWORK >= 2*N+1. This is the minimal requirement.
288 *> ->> For optimal performance (blocked code) the optimal value
289 *> is LWORK >= N + (N+1)*NB. Here NB is the optimal
290 *> block size for ZGEQP3 and ZGEQRF.
291 *> In general, optimal LWORK is computed as
292 *> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ)).
293 *> 1.2. .. an estimate of the scaled condition number of A is
294 *> required (JOBA='E', or 'G'). In this case, LWORK the minimal
295 *> requirement is LWORK >= N*N + 2*N.
296 *> ->> For optimal performance (blocked code) the optimal value
297 *> is LWORK >= max(N+(N+1)*NB, N*N+2*N)=N**2+2*N.
298 *> In general, the optimal length LWORK is computed as
299 *> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ),
300 *> N*N+LWORK(ZPOCON)).
301 *> 2. If SIGMA and the right singular vectors are needed (JOBV = 'V'),
302 *> (JOBU = 'N')
303 *> 2.1 .. no scaled condition estimate requested (JOBE = 'N'):
304 *> -> the minimal requirement is LWORK >= 3*N.
305 *> -> For optimal performance,
306 *> LWORK >= max(N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
307 *> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQ,
308 *> ZUNMLQ. In general, the optimal length LWORK is computed as
309 *> LWORK >= max(N+LWORK(ZGEQP3), N+LWORK(ZGESVJ),
310 *> N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)).
311 *> 2.2 .. an estimate of the scaled condition number of A is
312 *> required (JOBA='E', or 'G').
313 *> -> the minimal requirement is LWORK >= 3*N.
314 *> -> For optimal performance,
315 *> LWORK >= max(N+(N+1)*NB, 2*N,2*N+N*NB)=2*N+N*NB,
316 *> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQ,
317 *> ZUNMLQ. In general, the optimal length LWORK is computed as
318 *> LWORK >= max(N+LWORK(ZGEQP3), LWORK(ZPOCON), N+LWORK(ZGESVJ),
319 *> N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)).
320 *> 3. If SIGMA and the left singular vectors are needed
321 *> 3.1 .. no scaled condition estimate requested (JOBE = 'N'):
322 *> -> the minimal requirement is LWORK >= 3*N.
323 *> -> For optimal performance:
324 *> if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
325 *> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR.
326 *> In general, the optimal length LWORK is computed as
327 *> LWORK >= max(N+LWORK(ZGEQP3), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)).
328 *> 3.2 .. an estimate of the scaled condition number of A is
329 *> required (JOBA='E', or 'G').
330 *> -> the minimal requirement is LWORK >= 3*N.
331 *> -> For optimal performance:
332 *> if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
333 *> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR.
334 *> In general, the optimal length LWORK is computed as
335 *> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZPOCON),
336 *> 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)).
337 *> 4. If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and
338 *> 4.1. if JOBV = 'V'
339 *> the minimal requirement is LWORK >= 5*N+2*N*N.
340 *> 4.2. if JOBV = 'J' the minimal requirement is
341 *> LWORK >= 4*N+N*N.
342 *> In both cases, the allocated CWORK can accommodate blocked runs
343 *> of ZGEQP3, ZGEQRF, ZGELQF, SUNMQR, ZUNMLQ.
344 *>
345 *> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or
346 *> LRWORK=-1), then on exit CWORK(1) contains the optimal and CWORK(2) contains the
347 *> minimal length of CWORK for the job parameters used in the call.
348 *> \endverbatim
349 *>
350 *> \param[out] RWORK
351 *> \verbatim
352 *> RWORK is DOUBLE PRECISION array, dimension (MAX(7,LWORK))
353 *> On exit,
354 *> RWORK(1) = Determines the scaling factor SCALE = RWORK(2) / RWORK(1)
355 *> such that SCALE*SVA(1:N) are the computed singular values
356 *> of A. (See the description of SVA().)
357 *> RWORK(2) = See the description of RWORK(1).
358 *> RWORK(3) = SCONDA is an estimate for the condition number of
359 *> column equilibrated A. (If JOBA = 'E' or 'G')
360 *> SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1).
361 *> It is computed using SPOCON. It holds
362 *> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
363 *> where R is the triangular factor from the QRF of A.
364 *> However, if R is truncated and the numerical rank is
365 *> determined to be strictly smaller than N, SCONDA is
366 *> returned as -1, thus indicating that the smallest
367 *> singular values might be lost.
368 *>
369 *> If full SVD is needed, the following two condition numbers are
370 *> useful for the analysis of the algorithm. They are provied for
371 *> a developer/implementer who is familiar with the details of
372 *> the method.
373 *>
374 *> RWORK(4) = an estimate of the scaled condition number of the
375 *> triangular factor in the first QR factorization.
376 *> RWORK(5) = an estimate of the scaled condition number of the
377 *> triangular factor in the second QR factorization.
378 *> The following two parameters are computed if JOBT = 'T'.
379 *> They are provided for a developer/implementer who is familiar
380 *> with the details of the method.
381 *> RWORK(6) = the entropy of A^* * A :: this is the Shannon entropy
382 *> of diag(A^* * A) / Trace(A^* * A) taken as point in the
383 *> probability simplex.
384 *> RWORK(7) = the entropy of A * A^*. (See the description of RWORK(6).)
385 *> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or
386 *> LRWORK=-1), then on exit RWORK(1) contains the required length of
387 *> RWORK for the job parameters used in the call.
388 *> \endverbatim
389 *>
390 *> \param[in] LRWORK
391 *> \verbatim
392 *> LRWORK is INTEGER
393 *> Length of RWORK to confirm proper allocation of workspace.
394 *> LRWORK depends on the job:
395 *>
396 *> 1. If only the singular values are requested i.e. if
397 *> LSAME(JOBU,'N') .AND. LSAME(JOBV,'N')
398 *> then:
399 *> 1.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
400 *> then: LRWORK = max( 7, 2 * M ).
401 *> 1.2. Otherwise, LRWORK = max( 7, N ).
402 *> 2. If singular values with the right singular vectors are requested
403 *> i.e. if
404 *> (LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) .AND.
405 *> .NOT.(LSAME(JOBU,'U').OR.LSAME(JOBU,'F'))
406 *> then:
407 *> 2.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
408 *> then LRWORK = max( 7, 2 * M ).
409 *> 2.2. Otherwise, LRWORK = max( 7, N ).
410 *> 3. If singular values with the left singular vectors are requested, i.e. if
411 *> (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND.
412 *> .NOT.(LSAME(JOBV,'V').OR.LSAME(JOBV,'J'))
413 *> then:
414 *> 3.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
415 *> then LRWORK = max( 7, 2 * M ).
416 *> 3.2. Otherwise, LRWORK = max( 7, N ).
417 *> 4. If singular values with both the left and the right singular vectors
418 *> are requested, i.e. if
419 *> (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND.
420 *> (LSAME(JOBV,'V').OR.LSAME(JOBV,'J'))
421 *> then:
422 *> 4.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
423 *> then LRWORK = max( 7, 2 * M ).
424 *> 4.2. Otherwise, LRWORK = max( 7, N ).
425 *>
426 *> If, on entry, LRWORK = -1 or LWORK=-1, a workspace query is assumed and
427 *> the length of RWORK is returned in RWORK(1).
428 *> \endverbatim
429 *>
430 *> \param[out] IWORK
431 *> \verbatim
432 *> IWORK is INTEGER array, of dimension at least 4, that further depends
433 *> on the job:
434 *>
435 *> 1. If only the singular values are requested then:
436 *> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
437 *> then the length of IWORK is N+M; otherwise the length of IWORK is N.
438 *> 2. If the singular values and the right singular vectors are requested then:
439 *> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
440 *> then the length of IWORK is N+M; otherwise the length of IWORK is N.
441 *> 3. If the singular values and the left singular vectors are requested then:
442 *> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
443 *> then the length of IWORK is N+M; otherwise the length of IWORK is N.
444 *> 4. If the singular values with both the left and the right singular vectors
445 *> are requested, then:
446 *> 4.1. If LSAME(JOBV,'J') the length of IWORK is determined as follows:
447 *> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
448 *> then the length of IWORK is N+M; otherwise the length of IWORK is N.
449 *> 4.2. If LSAME(JOBV,'V') the length of IWORK is determined as follows:
450 *> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
451 *> then the length of IWORK is 2*N+M; otherwise the length of IWORK is 2*N.
452 *>
453 *> On exit,
454 *> IWORK(1) = the numerical rank determined after the initial
455 *> QR factorization with pivoting. See the descriptions
456 *> of JOBA and JOBR.
457 *> IWORK(2) = the number of the computed nonzero singular values
458 *> IWORK(3) = if nonzero, a warning message:
459 *> If IWORK(3) = 1 then some of the column norms of A
460 *> were denormalized floats. The requested high accuracy
461 *> is not warranted by the data.
462 *> IWORK(4) = 1 or -1. If IWORK(4) = 1, then the procedure used A^* to
463 *> do the job as specified by the JOB parameters.
464 *> If the call to ZGEJSV is a workspace query (indicated by LWORK = -1 or
465 *> LRWORK = -1), then on exit IWORK(1) contains the required length of
466 *> IWORK for the job parameters used in the call.
467 *> \endverbatim
468 *>
469 *> \param[out] INFO
470 *> \verbatim
471 *> INFO is INTEGER
472 *> < 0: if INFO = -i, then the i-th argument had an illegal value.
473 *> = 0: successful exit;
474 *> > 0: ZGEJSV did not converge in the maximal allowed number
475 *> of sweeps. The computed values may be inaccurate.
476 *> \endverbatim
477 *
478 * Authors:
479 * ========
480 *
481 *> \author Univ. of Tennessee
482 *> \author Univ. of California Berkeley
483 *> \author Univ. of Colorado Denver
484 *> \author NAG Ltd.
485 *
486 *> \date June 2016
487 *
488 *> \ingroup complex16GEsing
489 *
490 *> \par Further Details:
491 * =====================
492 *>
493 *> \verbatim
494 *>
495 *> ZGEJSV implements a preconditioned Jacobi SVD algorithm. It uses ZGEQP3,
496 *> ZGEQRF, and ZGELQF as preprocessors and preconditioners. Optionally, an
497 *> additional row pivoting can be used as a preprocessor, which in some
498 *> cases results in much higher accuracy. An example is matrix A with the
499 *> structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
500 *> diagonal matrices and C is well-conditioned matrix. In that case, complete
501 *> pivoting in the first QR factorizations provides accuracy dependent on the
502 *> condition number of C, and independent of D1, D2. Such higher accuracy is
503 *> not completely understood theoretically, but it works well in practice.
504 *> Further, if A can be written as A = B*D, with well-conditioned B and some
505 *> diagonal D, then the high accuracy is guaranteed, both theoretically and
506 *> in software, independent of D. For more details see [1], [2].
507 *> The computational range for the singular values can be the full range
508 *> ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
509 *> & LAPACK routines called by ZGEJSV are implemented to work in that range.
510 *> If that is not the case, then the restriction for safe computation with
511 *> the singular values in the range of normalized IEEE numbers is that the
512 *> spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
513 *> overflow. This code (ZGEJSV) is best used in this restricted range,
514 *> meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are
515 *> returned as zeros. See JOBR for details on this.
516 *> Further, this implementation is somewhat slower than the one described
517 *> in [1,2] due to replacement of some non-LAPACK components, and because
518 *> the choice of some tuning parameters in the iterative part (ZGESVJ) is
519 *> left to the implementer on a particular machine.
520 *> The rank revealing QR factorization (in this code: ZGEQP3) should be
521 *> implemented as in [3]. We have a new version of ZGEQP3 under development
522 *> that is more robust than the current one in LAPACK, with a cleaner cut in
523 *> rank deficient cases. It will be available in the SIGMA library [4].
524 *> If M is much larger than N, it is obvious that the initial QRF with
525 *> column pivoting can be preprocessed by the QRF without pivoting. That
526 *> well known trick is not used in ZGEJSV because in some cases heavy row
527 *> weighting can be treated with complete pivoting. The overhead in cases
528 *> M much larger than N is then only due to pivoting, but the benefits in
529 *> terms of accuracy have prevailed. The implementer/user can incorporate
530 *> this extra QRF step easily. The implementer can also improve data movement
531 *> (matrix transpose, matrix copy, matrix transposed copy) - this
532 *> implementation of ZGEJSV uses only the simplest, naive data movement.
533 *> \endverbatim
534 *
535 *> \par Contributor:
536 * ==================
537 *>
538 *> Zlatko Drmac, Department of Mathematics, Faculty of Science,
539 *> University of Zagreb (Zagreb, Croatia); drmac@math.hr
540 *
541 *> \par References:
542 * ================
543 *>
544 *> \verbatim
545 *>
546 *> [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
547 *> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
548 *> LAPACK Working note 169.
549 *> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
550 *> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
551 *> LAPACK Working note 170.
552 *> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
553 *> factorization software - a case study.
554 *> ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
555 *> LAPACK Working note 176.
556 *> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
557 *> QSVD, (H,K)-SVD computations.
558 *> Department of Mathematics, University of Zagreb, 2008, 2016.
559 *> \endverbatim
560 *
561 *> \par Bugs, examples and comments:
562 * =================================
563 *>
564 *> Please report all bugs and send interesting examples and/or comments to
565 *> drmac@math.hr. Thank you.
566 *>
567 * =====================================================================
568  SUBROUTINE zgejsv( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
569  $ M, N, A, LDA, SVA, U, LDU, V, LDV,
570  $ CWORK, LWORK, RWORK, LRWORK, IWORK, INFO )
571 *
572 * -- LAPACK computational routine (version 3.7.1) --
573 * -- LAPACK is a software package provided by Univ. of Tennessee, --
574 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
575 * June 2017
576 *
577 * .. Scalar Arguments ..
578  IMPLICIT NONE
579  INTEGER INFO, LDA, LDU, LDV, LWORK, LRWORK, M, N
580 * ..
581 * .. Array Arguments ..
582  COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ),
583  $ CWORK( LWORK )
584  DOUBLE PRECISION SVA( N ), RWORK( LRWORK )
585  INTEGER IWORK( * )
586  CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
587 * ..
588 *
589 * ===========================================================================
590 *
591 * .. Local Parameters ..
592  DOUBLE PRECISION ZERO, ONE
593  PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
594  COMPLEX*16 CZERO, CONE
595  parameter( czero = ( 0.0d0, 0.0d0 ), cone = ( 1.0d0, 0.0d0 ) )
596 * ..
597 * .. Local Scalars ..
598  COMPLEX*16 CTEMP
599  DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1,
600  $ cond_ok, condr1, condr2, entra, entrat, epsln,
601  $ maxprj, scalem, sconda, sfmin, small, temp1,
602  $ uscal1, uscal2, xsc
603  INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING
604  LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LQUERY,
605  $ lsvec, l2aber, l2kill, l2pert, l2rank, l2tran, noscal,
606  $ rowpiv, rsvec, transp
607 *
608  INTEGER OPTWRK, MINWRK, MINRWRK, MINIWRK
609  INTEGER LWCON, LWLQF, LWQP3, LWQRF, LWUNMLQ, LWUNMQR, LWUNMQRM,
610  $ lwsvdj, lwsvdjv, lrwqp3, lrwcon, lrwsvdj, iwoff
611  INTEGER LWRK_ZGELQF, LWRK_ZGEQP3, LWRK_ZGEQP3N, LWRK_ZGEQRF,
612  $ LWRK_ZGESVJ, LWRK_ZGESVJV, LWRK_ZGESVJU, LWRK_ZUNMLQ,
613  $ lwrk_zunmqr, lwrk_zunmqrm
614 * ..
615 * .. Local Arrays
616  COMPLEX*16 CDUMMY(1)
617  DOUBLE PRECISION RDUMMY(1)
618 *
619 * .. Intrinsic Functions ..
620  INTRINSIC abs, dcmplx, conjg, dlog, max, min, dble, nint, sqrt
621 * ..
622 * .. External Functions ..
623  DOUBLE PRECISION DLAMCH, DZNRM2
624  INTEGER IDAMAX, IZAMAX
625  LOGICAL LSAME
626  EXTERNAL idamax, izamax, lsame, dlamch, dznrm2
627 * ..
628 * .. External Subroutines ..
629  EXTERNAL dlassq, zcopy, zgelqf, zgeqp3, zgeqrf, zlacpy, zlapmr,
630  $ zlascl, dlascl, zlaset, zlassq, zlaswp, zungqr, zunmlq,
631  $ zunmqr, zpocon, dscal, zdscal, zswap, ztrsm, zlacgv,
632  $ xerbla
633 *
634  EXTERNAL zgesvj
635 * ..
636 *
637 * Test the input arguments
638 *
639  lsvec = lsame( jobu, 'U' ) .OR. lsame( jobu, 'F' )
640  jracc = lsame( jobv, 'J' )
641  rsvec = lsame( jobv, 'V' ) .OR. jracc
642  rowpiv = lsame( joba, 'F' ) .OR. lsame( joba, 'G' )
643  l2rank = lsame( joba, 'R' )
644  l2aber = lsame( joba, 'A' )
645  errest = lsame( joba, 'E' ) .OR. lsame( joba, 'G' )
646  l2tran = lsame( jobt, 'T' ) .AND. ( m .EQ. n )
647  l2kill = lsame( jobr, 'R' )
648  defr = lsame( jobr, 'N' )
649  l2pert = lsame( jobp, 'P' )
650 *
651  lquery = ( lwork .EQ. -1 ) .OR. ( lrwork .EQ. -1 )
652 *
653  IF ( .NOT.(rowpiv .OR. l2rank .OR. l2aber .OR.
654  $ errest .OR. lsame( joba, 'C' ) )) THEN
655  info = - 1
656  ELSE IF ( .NOT.( lsvec .OR. lsame( jobu, 'N' ) .OR.
657  $ ( lsame( jobu, 'W' ) .AND. rsvec .AND. l2tran ) ) ) THEN
658  info = - 2
659  ELSE IF ( .NOT.( rsvec .OR. lsame( jobv, 'N' ) .OR.
660  $ ( lsame( jobv, 'W' ) .AND. lsvec .AND. l2tran ) ) ) THEN
661  info = - 3
662  ELSE IF ( .NOT. ( l2kill .OR. defr ) ) THEN
663  info = - 4
664  ELSE IF ( .NOT. ( lsame(jobt,'T') .OR. lsame(jobt,'N') ) ) THEN
665  info = - 5
666  ELSE IF ( .NOT. ( l2pert .OR. lsame( jobp, 'N' ) ) ) THEN
667  info = - 6
668  ELSE IF ( m .LT. 0 ) THEN
669  info = - 7
670  ELSE IF ( ( n .LT. 0 ) .OR. ( n .GT. m ) ) THEN
671  info = - 8
672  ELSE IF ( lda .LT. m ) THEN
673  info = - 10
674  ELSE IF ( lsvec .AND. ( ldu .LT. m ) ) THEN
675  info = - 13
676  ELSE IF ( rsvec .AND. ( ldv .LT. n ) ) THEN
677  info = - 15
678  ELSE
679 * #:)
680  info = 0
681  END IF
682 *
683  IF ( info .EQ. 0 ) THEN
684 * .. compute the minimal and the optimal workspace lengths
685 * [[The expressions for computing the minimal and the optimal
686 * values of LCWORK, LRWORK are written with a lot of redundancy and
687 * can be simplified. However, this verbose form is useful for
688 * maintenance and modifications of the code.]]
689 *
690 * .. minimal workspace length for ZGEQP3 of an M x N matrix,
691 * ZGEQRF of an N x N matrix, ZGELQF of an N x N matrix,
692 * ZUNMLQ for computing N x N matrix, ZUNMQR for computing N x N
693 * matrix, ZUNMQR for computing M x N matrix, respectively.
694  lwqp3 = n+1
695  lwqrf = max( 1, n )
696  lwlqf = max( 1, n )
697  lwunmlq = max( 1, n )
698  lwunmqr = max( 1, n )
699  lwunmqrm = max( 1, m )
700 * .. minimal workspace length for ZPOCON of an N x N matrix
701  lwcon = 2 * n
702 * .. minimal workspace length for ZGESVJ of an N x N matrix,
703 * without and with explicit accumulation of Jacobi rotations
704  lwsvdj = max( 2 * n, 1 )
705  lwsvdjv = max( 2 * n, 1 )
706 * .. minimal REAL workspace length for ZGEQP3, ZPOCON, ZGESVJ
707  lrwqp3 = 2 * n
708  lrwcon = n
709  lrwsvdj = n
710  IF ( lquery ) THEN
711  CALL zgeqp3( m, n, a, lda, iwork, cdummy, cdummy, -1,
712  $ rdummy, ierr )
713  lwrk_zgeqp3 = cdummy(1)
714  CALL zgeqrf( n, n, a, lda, cdummy, cdummy,-1, ierr )
715  lwrk_zgeqrf = cdummy(1)
716  CALL zgelqf( n, n, a, lda, cdummy, cdummy,-1, ierr )
717  lwrk_zgelqf = cdummy(1)
718  END IF
719  minwrk = 2
720  optwrk = 2
721  miniwrk = n
722  IF ( .NOT. (lsvec .OR. rsvec ) ) THEN
723 * .. minimal and optimal sizes of the complex workspace if
724 * only the singular values are requested
725  IF ( errest ) THEN
726  minwrk = max( n+lwqp3, n**2+lwcon, n+lwqrf, lwsvdj )
727  ELSE
728  minwrk = max( n+lwqp3, n+lwqrf, lwsvdj )
729  END IF
730  IF ( lquery ) THEN
731  CALL zgesvj( 'L', 'N', 'N', n, n, a, lda, sva, n, v,
732  $ ldv, cdummy, -1, rdummy, -1, ierr )
733  lwrk_zgesvj = cdummy(1)
734  IF ( errest ) THEN
735  optwrk = max( n+lwrk_zgeqp3, n**2+lwcon,
736  $ n+lwrk_zgeqrf, lwrk_zgesvj )
737  ELSE
738  optwrk = max( n+lwrk_zgeqp3, n+lwrk_zgeqrf,
739  $ lwrk_zgesvj )
740  END IF
741  END IF
742  IF ( l2tran .OR. rowpiv ) THEN
743  IF ( errest ) THEN
744  minrwrk = max( 7, 2*m, lrwqp3, lrwcon, lrwsvdj )
745  ELSE
746  minrwrk = max( 7, 2*m, lrwqp3, lrwsvdj )
747  END IF
748  ELSE
749  IF ( errest ) THEN
750  minrwrk = max( 7, lrwqp3, lrwcon, lrwsvdj )
751  ELSE
752  minrwrk = max( 7, lrwqp3, lrwsvdj )
753  END IF
754  END IF
755  IF ( rowpiv .OR. l2tran ) miniwrk = miniwrk + m
756  ELSE IF ( rsvec .AND. (.NOT.lsvec) ) THEN
757 * .. minimal and optimal sizes of the complex workspace if the
758 * singular values and the right singular vectors are requested
759  IF ( errest ) THEN
760  minwrk = max( n+lwqp3, lwcon, lwsvdj, n+lwlqf,
761  $ 2*n+lwqrf, n+lwsvdj, n+lwunmlq )
762  ELSE
763  minwrk = max( n+lwqp3, lwsvdj, n+lwlqf, 2*n+lwqrf,
764  $ n+lwsvdj, n+lwunmlq )
765  END IF
766  IF ( lquery ) THEN
767  CALL zgesvj( 'L', 'U', 'N', n,n, u, ldu, sva, n, a,
768  $ lda, cdummy, -1, rdummy, -1, ierr )
769  lwrk_zgesvj = cdummy(1)
770  CALL zunmlq( 'L', 'C', n, n, n, a, lda, cdummy,
771  $ v, ldv, cdummy, -1, ierr )
772  lwrk_zunmlq = cdummy(1)
773  IF ( errest ) THEN
774  optwrk = max( n+lwrk_zgeqp3, lwcon, lwrk_zgesvj,
775  $ n+lwrk_zgelqf, 2*n+lwrk_zgeqrf,
776  $ n+lwrk_zgesvj, n+lwrk_zunmlq )
777  ELSE
778  optwrk = max( n+lwrk_zgeqp3, lwrk_zgesvj,n+lwrk_zgelqf,
779  $ 2*n+lwrk_zgeqrf, n+lwrk_zgesvj,
780  $ n+lwrk_zunmlq )
781  END IF
782  END IF
783  IF ( l2tran .OR. rowpiv ) THEN
784  IF ( errest ) THEN
785  minrwrk = max( 7, 2*m, lrwqp3, lrwsvdj, lrwcon )
786  ELSE
787  minrwrk = max( 7, 2*m, lrwqp3, lrwsvdj )
788  END IF
789  ELSE
790  IF ( errest ) THEN
791  minrwrk = max( 7, lrwqp3, lrwsvdj, lrwcon )
792  ELSE
793  minrwrk = max( 7, lrwqp3, lrwsvdj )
794  END IF
795  END IF
796  IF ( rowpiv .OR. l2tran ) miniwrk = miniwrk + m
797  ELSE IF ( lsvec .AND. (.NOT.rsvec) ) THEN
798 * .. minimal and optimal sizes of the complex workspace if the
799 * singular values and the left singular vectors are requested
800  IF ( errest ) THEN
801  minwrk = n + max( lwqp3,lwcon,n+lwqrf,lwsvdj,lwunmqrm )
802  ELSE
803  minwrk = n + max( lwqp3, n+lwqrf, lwsvdj, lwunmqrm )
804  END IF
805  IF ( lquery ) THEN
806  CALL zgesvj( 'L', 'U', 'N', n,n, u, ldu, sva, n, a,
807  $ lda, cdummy, -1, rdummy, -1, ierr )
808  lwrk_zgesvj = cdummy(1)
809  CALL zunmqr( 'L', 'N', m, n, n, a, lda, cdummy, u,
810  $ ldu, cdummy, -1, ierr )
811  lwrk_zunmqrm = cdummy(1)
812  IF ( errest ) THEN
813  optwrk = n + max( lwrk_zgeqp3, lwcon, n+lwrk_zgeqrf,
814  $ lwrk_zgesvj, lwrk_zunmqrm )
815  ELSE
816  optwrk = n + max( lwrk_zgeqp3, n+lwrk_zgeqrf,
817  $ lwrk_zgesvj, lwrk_zunmqrm )
818  END IF
819  END IF
820  IF ( l2tran .OR. rowpiv ) THEN
821  IF ( errest ) THEN
822  minrwrk = max( 7, 2*m, lrwqp3, lrwsvdj, lrwcon )
823  ELSE
824  minrwrk = max( 7, 2*m, lrwqp3, lrwsvdj )
825  END IF
826  ELSE
827  IF ( errest ) THEN
828  minrwrk = max( 7, lrwqp3, lrwsvdj, lrwcon )
829  ELSE
830  minrwrk = max( 7, lrwqp3, lrwsvdj )
831  END IF
832  END IF
833  IF ( rowpiv .OR. l2tran ) miniwrk = miniwrk + m
834  ELSE
835 * .. minimal and optimal sizes of the complex workspace if the
836 * full SVD is requested
837  IF ( .NOT. jracc ) THEN
838  IF ( errest ) THEN
839  minwrk = max( n+lwqp3, n+lwcon, 2*n+n**2+lwcon,
840  $ 2*n+lwqrf, 2*n+lwqp3,
841  $ 2*n+n**2+n+lwlqf, 2*n+n**2+n+n**2+lwcon,
842  $ 2*n+n**2+n+lwsvdj, 2*n+n**2+n+lwsvdjv,
843  $ 2*n+n**2+n+lwunmqr,2*n+n**2+n+lwunmlq,
844  $ n+n**2+lwsvdj, n+lwunmqrm )
845  ELSE
846  minwrk = max( n+lwqp3, 2*n+n**2+lwcon,
847  $ 2*n+lwqrf, 2*n+lwqp3,
848  $ 2*n+n**2+n+lwlqf, 2*n+n**2+n+n**2+lwcon,
849  $ 2*n+n**2+n+lwsvdj, 2*n+n**2+n+lwsvdjv,
850  $ 2*n+n**2+n+lwunmqr,2*n+n**2+n+lwunmlq,
851  $ n+n**2+lwsvdj, n+lwunmqrm )
852  END IF
853  miniwrk = miniwrk + n
854  IF ( rowpiv .OR. l2tran ) miniwrk = miniwrk + m
855  ELSE
856  IF ( errest ) THEN
857  minwrk = max( n+lwqp3, n+lwcon, 2*n+lwqrf,
858  $ 2*n+n**2+lwsvdjv, 2*n+n**2+n+lwunmqr,
859  $ n+lwunmqrm )
860  ELSE
861  minwrk = max( n+lwqp3, 2*n+lwqrf,
862  $ 2*n+n**2+lwsvdjv, 2*n+n**2+n+lwunmqr,
863  $ n+lwunmqrm )
864  END IF
865  IF ( rowpiv .OR. l2tran ) miniwrk = miniwrk + m
866  END IF
867  IF ( lquery ) THEN
868  CALL zunmqr( 'L', 'N', m, n, n, a, lda, cdummy, u,
869  $ ldu, cdummy, -1, ierr )
870  lwrk_zunmqrm = cdummy(1)
871  CALL zunmqr( 'L', 'N', n, n, n, a, lda, cdummy, u,
872  $ ldu, cdummy, -1, ierr )
873  lwrk_zunmqr = cdummy(1)
874  IF ( .NOT. jracc ) THEN
875  CALL zgeqp3( n,n, a, lda, iwork, cdummy,cdummy, -1,
876  $ rdummy, ierr )
877  lwrk_zgeqp3n = cdummy(1)
878  CALL zgesvj( 'L', 'U', 'N', n, n, u, ldu, sva,
879  $ n, v, ldv, cdummy, -1, rdummy, -1, ierr )
880  lwrk_zgesvj = cdummy(1)
881  CALL zgesvj( 'U', 'U', 'N', n, n, u, ldu, sva,
882  $ n, v, ldv, cdummy, -1, rdummy, -1, ierr )
883  lwrk_zgesvju = cdummy(1)
884  CALL zgesvj( 'L', 'U', 'V', n, n, u, ldu, sva,
885  $ n, v, ldv, cdummy, -1, rdummy, -1, ierr )
886  lwrk_zgesvjv = cdummy(1)
887  CALL zunmlq( 'L', 'C', n, n, n, a, lda, cdummy,
888  $ v, ldv, cdummy, -1, ierr )
889  lwrk_zunmlq = cdummy(1)
890  IF ( errest ) THEN
891  optwrk = max( n+lwrk_zgeqp3, n+lwcon,
892  $ 2*n+n**2+lwcon, 2*n+lwrk_zgeqrf,
893  $ 2*n+lwrk_zgeqp3n,
894  $ 2*n+n**2+n+lwrk_zgelqf,
895  $ 2*n+n**2+n+n**2+lwcon,
896  $ 2*n+n**2+n+lwrk_zgesvj,
897  $ 2*n+n**2+n+lwrk_zgesvjv,
898  $ 2*n+n**2+n+lwrk_zunmqr,
899  $ 2*n+n**2+n+lwrk_zunmlq,
900  $ n+n**2+lwrk_zgesvju,
901  $ n+lwrk_zunmqrm )
902  ELSE
903  optwrk = max( n+lwrk_zgeqp3,
904  $ 2*n+n**2+lwcon, 2*n+lwrk_zgeqrf,
905  $ 2*n+lwrk_zgeqp3n,
906  $ 2*n+n**2+n+lwrk_zgelqf,
907  $ 2*n+n**2+n+n**2+lwcon,
908  $ 2*n+n**2+n+lwrk_zgesvj,
909  $ 2*n+n**2+n+lwrk_zgesvjv,
910  $ 2*n+n**2+n+lwrk_zunmqr,
911  $ 2*n+n**2+n+lwrk_zunmlq,
912  $ n+n**2+lwrk_zgesvju,
913  $ n+lwrk_zunmqrm )
914  END IF
915  ELSE
916  CALL zgesvj( 'L', 'U', 'V', n, n, u, ldu, sva,
917  $ n, v, ldv, cdummy, -1, rdummy, -1, ierr )
918  lwrk_zgesvjv = cdummy(1)
919  CALL zunmqr( 'L', 'N', n, n, n, cdummy, n, cdummy,
920  $ v, ldv, cdummy, -1, ierr )
921  lwrk_zunmqr = cdummy(1)
922  CALL zunmqr( 'L', 'N', m, n, n, a, lda, cdummy, u,
923  $ ldu, cdummy, -1, ierr )
924  lwrk_zunmqrm = cdummy(1)
925  IF ( errest ) THEN
926  optwrk = max( n+lwrk_zgeqp3, n+lwcon,
927  $ 2*n+lwrk_zgeqrf, 2*n+n**2,
928  $ 2*n+n**2+lwrk_zgesvjv,
929  $ 2*n+n**2+n+lwrk_zunmqr,n+lwrk_zunmqrm )
930  ELSE
931  optwrk = max( n+lwrk_zgeqp3, 2*n+lwrk_zgeqrf,
932  $ 2*n+n**2, 2*n+n**2+lwrk_zgesvjv,
933  $ 2*n+n**2+n+lwrk_zunmqr,
934  $ n+lwrk_zunmqrm )
935  END IF
936  END IF
937  END IF
938  IF ( l2tran .OR. rowpiv ) THEN
939  minrwrk = max( 7, 2*m, lrwqp3, lrwsvdj, lrwcon )
940  ELSE
941  minrwrk = max( 7, lrwqp3, lrwsvdj, lrwcon )
942  END IF
943  END IF
944  minwrk = max( 2, minwrk )
945  optwrk = max( minwrk, optwrk )
946  IF ( lwork .LT. minwrk .AND. (.NOT.lquery) ) info = - 17
947  IF ( lrwork .LT. minrwrk .AND. (.NOT.lquery) ) info = - 19
948  END IF
949 *
950  IF ( info .NE. 0 ) THEN
951 * #:(
952  CALL xerbla( 'ZGEJSV', - info )
953  RETURN
954  ELSE IF ( lquery ) THEN
955  cwork(1) = optwrk
956  cwork(2) = minwrk
957  rwork(1) = minrwrk
958  iwork(1) = max( 4, miniwrk )
959  RETURN
960  END IF
961 *
962 * Quick return for void matrix (Y3K safe)
963 * #:)
964  IF ( ( m .EQ. 0 ) .OR. ( n .EQ. 0 ) ) THEN
965  iwork(1:4) = 0
966  rwork(1:7) = 0
967  RETURN
968  ENDIF
969 *
970 * Determine whether the matrix U should be M x N or M x M
971 *
972  IF ( lsvec ) THEN
973  n1 = n
974  IF ( lsame( jobu, 'F' ) ) n1 = m
975  END IF
976 *
977 * Set numerical parameters
978 *
979 *! NOTE: Make sure DLAMCH() does not fail on the target architecture.
980 *
981  epsln = dlamch('Epsilon')
982  sfmin = dlamch('SafeMinimum')
983  small = sfmin / epsln
984  big = dlamch('O')
985 * BIG = ONE / SFMIN
986 *
987 * Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
988 *
989 *(!) If necessary, scale SVA() to protect the largest norm from
990 * overflow. It is possible that this scaling pushes the smallest
991 * column norm left from the underflow threshold (extreme case).
992 *
993  scalem = one / sqrt(dble(m)*dble(n))
994  noscal = .true.
995  goscal = .true.
996  DO 1874 p = 1, n
997  aapp = zero
998  aaqq = one
999  CALL zlassq( m, a(1,p), 1, aapp, aaqq )
1000  IF ( aapp .GT. big ) THEN
1001  info = - 9
1002  CALL xerbla( 'ZGEJSV', -info )
1003  RETURN
1004  END IF
1005  aaqq = sqrt(aaqq)
1006  IF ( ( aapp .LT. (big / aaqq) ) .AND. noscal ) THEN
1007  sva(p) = aapp * aaqq
1008  ELSE
1009  noscal = .false.
1010  sva(p) = aapp * ( aaqq * scalem )
1011  IF ( goscal ) THEN
1012  goscal = .false.
1013  CALL dscal( p-1, scalem, sva, 1 )
1014  END IF
1015  END IF
1016  1874 CONTINUE
1017 *
1018  IF ( noscal ) scalem = one
1019 *
1020  aapp = zero
1021  aaqq = big
1022  DO 4781 p = 1, n
1023  aapp = max( aapp, sva(p) )
1024  IF ( sva(p) .NE. zero ) aaqq = min( aaqq, sva(p) )
1025  4781 CONTINUE
1026 *
1027 * Quick return for zero M x N matrix
1028 * #:)
1029  IF ( aapp .EQ. zero ) THEN
1030  IF ( lsvec ) CALL zlaset( 'G', m, n1, czero, cone, u, ldu )
1031  IF ( rsvec ) CALL zlaset( 'G', n, n, czero, cone, v, ldv )
1032  rwork(1) = one
1033  rwork(2) = one
1034  IF ( errest ) rwork(3) = one
1035  IF ( lsvec .AND. rsvec ) THEN
1036  rwork(4) = one
1037  rwork(5) = one
1038  END IF
1039  IF ( l2tran ) THEN
1040  rwork(6) = zero
1041  rwork(7) = zero
1042  END IF
1043  iwork(1) = 0
1044  iwork(2) = 0
1045  iwork(3) = 0
1046  iwork(4) = -1
1047  RETURN
1048  END IF
1049 *
1050 * Issue warning if denormalized column norms detected. Override the
1051 * high relative accuracy request. Issue licence to kill nonzero columns
1052 * (set them to zero) whose norm is less than sigma_max / BIG (roughly).
1053 * #:(
1054  warning = 0
1055  IF ( aaqq .LE. sfmin ) THEN
1056  l2rank = .true.
1057  l2kill = .true.
1058  warning = 1
1059  END IF
1060 *
1061 * Quick return for one-column matrix
1062 * #:)
1063  IF ( n .EQ. 1 ) THEN
1064 *
1065  IF ( lsvec ) THEN
1066  CALL zlascl( 'G',0,0,sva(1),scalem, m,1,a(1,1),lda,ierr )
1067  CALL zlacpy( 'A', m, 1, a, lda, u, ldu )
1068 * computing all M left singular vectors of the M x 1 matrix
1069  IF ( n1 .NE. n ) THEN
1070  CALL zgeqrf( m, n, u,ldu, cwork, cwork(n+1),lwork-n,ierr )
1071  CALL zungqr( m,n1,1, u,ldu,cwork,cwork(n+1),lwork-n,ierr )
1072  CALL zcopy( m, a(1,1), 1, u(1,1), 1 )
1073  END IF
1074  END IF
1075  IF ( rsvec ) THEN
1076  v(1,1) = cone
1077  END IF
1078  IF ( sva(1) .LT. (big*scalem) ) THEN
1079  sva(1) = sva(1) / scalem
1080  scalem = one
1081  END IF
1082  rwork(1) = one / scalem
1083  rwork(2) = one
1084  IF ( sva(1) .NE. zero ) THEN
1085  iwork(1) = 1
1086  IF ( ( sva(1) / scalem) .GE. sfmin ) THEN
1087  iwork(2) = 1
1088  ELSE
1089  iwork(2) = 0
1090  END IF
1091  ELSE
1092  iwork(1) = 0
1093  iwork(2) = 0
1094  END IF
1095  iwork(3) = 0
1096  iwork(4) = -1
1097  IF ( errest ) rwork(3) = one
1098  IF ( lsvec .AND. rsvec ) THEN
1099  rwork(4) = one
1100  rwork(5) = one
1101  END IF
1102  IF ( l2tran ) THEN
1103  rwork(6) = zero
1104  rwork(7) = zero
1105  END IF
1106  RETURN
1107 *
1108  END IF
1109 *
1110  transp = .false.
1111 *
1112  aatmax = -one
1113  aatmin = big
1114  IF ( rowpiv .OR. l2tran ) THEN
1115 *
1116 * Compute the row norms, needed to determine row pivoting sequence
1117 * (in the case of heavily row weighted A, row pivoting is strongly
1118 * advised) and to collect information needed to compare the
1119 * structures of A * A^* and A^* * A (in the case L2TRAN.EQ..TRUE.).
1120 *
1121  IF ( l2tran ) THEN
1122  DO 1950 p = 1, m
1123  xsc = zero
1124  temp1 = one
1125  CALL zlassq( n, a(p,1), lda, xsc, temp1 )
1126 * ZLASSQ gets both the ell_2 and the ell_infinity norm
1127 * in one pass through the vector
1128  rwork(m+p) = xsc * scalem
1129  rwork(p) = xsc * (scalem*sqrt(temp1))
1130  aatmax = max( aatmax, rwork(p) )
1131  IF (rwork(p) .NE. zero)
1132  $ aatmin = min(aatmin,rwork(p))
1133  1950 CONTINUE
1134  ELSE
1135  DO 1904 p = 1, m
1136  rwork(m+p) = scalem*abs( a(p,izamax(n,a(p,1),lda)) )
1137  aatmax = max( aatmax, rwork(m+p) )
1138  aatmin = min( aatmin, rwork(m+p) )
1139  1904 CONTINUE
1140  END IF
1141 *
1142  END IF
1143 *
1144 * For square matrix A try to determine whether A^* would be better
1145 * input for the preconditioned Jacobi SVD, with faster convergence.
1146 * The decision is based on an O(N) function of the vector of column
1147 * and row norms of A, based on the Shannon entropy. This should give
1148 * the right choice in most cases when the difference actually matters.
1149 * It may fail and pick the slower converging side.
1150 *
1151  entra = zero
1152  entrat = zero
1153  IF ( l2tran ) THEN
1154 *
1155  xsc = zero
1156  temp1 = one
1157  CALL dlassq( n, sva, 1, xsc, temp1 )
1158  temp1 = one / temp1
1159 *
1160  entra = zero
1161  DO 1113 p = 1, n
1162  big1 = ( ( sva(p) / xsc )**2 ) * temp1
1163  IF ( big1 .NE. zero ) entra = entra + big1 * dlog(big1)
1164  1113 CONTINUE
1165  entra = - entra / dlog(dble(n))
1166 *
1167 * Now, SVA().^2/Trace(A^* * A) is a point in the probability simplex.
1168 * It is derived from the diagonal of A^* * A. Do the same with the
1169 * diagonal of A * A^*, compute the entropy of the corresponding
1170 * probability distribution. Note that A * A^* and A^* * A have the
1171 * same trace.
1172 *
1173  entrat = zero
1174  DO 1114 p = 1, m
1175  big1 = ( ( rwork(p) / xsc )**2 ) * temp1
1176  IF ( big1 .NE. zero ) entrat = entrat + big1 * dlog(big1)
1177  1114 CONTINUE
1178  entrat = - entrat / dlog(dble(m))
1179 *
1180 * Analyze the entropies and decide A or A^*. Smaller entropy
1181 * usually means better input for the algorithm.
1182 *
1183  transp = ( entrat .LT. entra )
1184 *
1185 * If A^* is better than A, take the adjoint of A. This is allowed
1186 * only for square matrices, M=N.
1187  IF ( transp ) THEN
1188 * In an optimal implementation, this trivial transpose
1189 * should be replaced with faster transpose.
1190  DO 1115 p = 1, n - 1
1191  a(p,p) = conjg(a(p,p))
1192  DO 1116 q = p + 1, n
1193  ctemp = conjg(a(q,p))
1194  a(q,p) = conjg(a(p,q))
1195  a(p,q) = ctemp
1196  1116 CONTINUE
1197  1115 CONTINUE
1198  a(n,n) = conjg(a(n,n))
1199  DO 1117 p = 1, n
1200  rwork(m+p) = sva(p)
1201  sva(p) = rwork(p)
1202 * previously computed row 2-norms are now column 2-norms
1203 * of the transposed matrix
1204  1117 CONTINUE
1205  temp1 = aapp
1206  aapp = aatmax
1207  aatmax = temp1
1208  temp1 = aaqq
1209  aaqq = aatmin
1210  aatmin = temp1
1211  kill = lsvec
1212  lsvec = rsvec
1213  rsvec = kill
1214  IF ( lsvec ) n1 = n
1215 *
1216  rowpiv = .true.
1217  END IF
1218 *
1219  END IF
1220 * END IF L2TRAN
1221 *
1222 * Scale the matrix so that its maximal singular value remains less
1223 * than SQRT(BIG) -- the matrix is scaled so that its maximal column
1224 * has Euclidean norm equal to SQRT(BIG/N). The only reason to keep
1225 * SQRT(BIG) instead of BIG is the fact that ZGEJSV uses LAPACK and
1226 * BLAS routines that, in some implementations, are not capable of
1227 * working in the full interval [SFMIN,BIG] and that they may provoke
1228 * overflows in the intermediate results. If the singular values spread
1229 * from SFMIN to BIG, then ZGESVJ will compute them. So, in that case,
1230 * one should use ZGESVJ instead of ZGEJSV.
1231 * >> change in the April 2016 update: allow bigger range, i.e. the
1232 * largest column is allowed up to BIG/N and ZGESVJ will do the rest.
1233  big1 = sqrt( big )
1234  temp1 = sqrt( big / dble(n) )
1235 * TEMP1 = BIG/DBLE(N)
1236 *
1237  CALL dlascl( 'G', 0, 0, aapp, temp1, n, 1, sva, n, ierr )
1238  IF ( aaqq .GT. (aapp * sfmin) ) THEN
1239  aaqq = ( aaqq / aapp ) * temp1
1240  ELSE
1241  aaqq = ( aaqq * temp1 ) / aapp
1242  END IF
1243  temp1 = temp1 * scalem
1244  CALL zlascl( 'G', 0, 0, aapp, temp1, m, n, a, lda, ierr )
1245 *
1246 * To undo scaling at the end of this procedure, multiply the
1247 * computed singular values with USCAL2 / USCAL1.
1248 *
1249  uscal1 = temp1
1250  uscal2 = aapp
1251 *
1252  IF ( l2kill ) THEN
1253 * L2KILL enforces computation of nonzero singular values in
1254 * the restricted range of condition number of the initial A,
1255 * sigma_max(A) / sigma_min(A) approx. SQRT(BIG)/SQRT(SFMIN).
1256  xsc = sqrt( sfmin )
1257  ELSE
1258  xsc = small
1259 *
1260 * Now, if the condition number of A is too big,
1261 * sigma_max(A) / sigma_min(A) .GT. SQRT(BIG/N) * EPSLN / SFMIN,
1262 * as a precaution measure, the full SVD is computed using ZGESVJ
1263 * with accumulated Jacobi rotations. This provides numerically
1264 * more robust computation, at the cost of slightly increased run
1265 * time. Depending on the concrete implementation of BLAS and LAPACK
1266 * (i.e. how they behave in presence of extreme ill-conditioning) the
1267 * implementor may decide to remove this switch.
1268  IF ( ( aaqq.LT.sqrt(sfmin) ) .AND. lsvec .AND. rsvec ) THEN
1269  jracc = .true.
1270  END IF
1271 *
1272  END IF
1273  IF ( aaqq .LT. xsc ) THEN
1274  DO 700 p = 1, n
1275  IF ( sva(p) .LT. xsc ) THEN
1276  CALL zlaset( 'A', m, 1, czero, czero, a(1,p), lda )
1277  sva(p) = zero
1278  END IF
1279  700 CONTINUE
1280  END IF
1281 *
1282 * Preconditioning using QR factorization with pivoting
1283 *
1284  IF ( rowpiv ) THEN
1285 * Optional row permutation (Bjoerck row pivoting):
1286 * A result by Cox and Higham shows that the Bjoerck's
1287 * row pivoting combined with standard column pivoting
1288 * has similar effect as Powell-Reid complete pivoting.
1289 * The ell-infinity norms of A are made nonincreasing.
1290  IF ( ( lsvec .AND. rsvec ) .AND. .NOT.( jracc ) ) THEN
1291  iwoff = 2*n
1292  ELSE
1293  iwoff = n
1294  END IF
1295  DO 1952 p = 1, m - 1
1296  q = idamax( m-p+1, rwork(m+p), 1 ) + p - 1
1297  iwork(iwoff+p) = q
1298  IF ( p .NE. q ) THEN
1299  temp1 = rwork(m+p)
1300  rwork(m+p) = rwork(m+q)
1301  rwork(m+q) = temp1
1302  END IF
1303  1952 CONTINUE
1304  CALL zlaswp( n, a, lda, 1, m-1, iwork(iwoff+1), 1 )
1305  END IF
1306 *
1307 * End of the preparation phase (scaling, optional sorting and
1308 * transposing, optional flushing of small columns).
1309 *
1310 * Preconditioning
1311 *
1312 * If the full SVD is needed, the right singular vectors are computed
1313 * from a matrix equation, and for that we need theoretical analysis
1314 * of the Businger-Golub pivoting. So we use ZGEQP3 as the first RR QRF.
1315 * In all other cases the first RR QRF can be chosen by other criteria
1316 * (eg speed by replacing global with restricted window pivoting, such
1317 * as in xGEQPX from TOMS # 782). Good results will be obtained using
1318 * xGEQPX with properly (!) chosen numerical parameters.
1319 * Any improvement of ZGEQP3 improves overal performance of ZGEJSV.
1320 *
1321 * A * P1 = Q1 * [ R1^* 0]^*:
1322  DO 1963 p = 1, n
1323 * .. all columns are free columns
1324  iwork(p) = 0
1325  1963 CONTINUE
1326  CALL zgeqp3( m, n, a, lda, iwork, cwork, cwork(n+1), lwork-n,
1327  $ rwork, ierr )
1328 *
1329 * The upper triangular matrix R1 from the first QRF is inspected for
1330 * rank deficiency and possibilities for deflation, or possible
1331 * ill-conditioning. Depending on the user specified flag L2RANK,
1332 * the procedure explores possibilities to reduce the numerical
1333 * rank by inspecting the computed upper triangular factor. If
1334 * L2RANK or L2ABER are up, then ZGEJSV will compute the SVD of
1335 * A + dA, where ||dA|| <= f(M,N)*EPSLN.
1336 *
1337  nr = 1
1338  IF ( l2aber ) THEN
1339 * Standard absolute error bound suffices. All sigma_i with
1340 * sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
1341 * aggressive enforcement of lower numerical rank by introducing a
1342 * backward error of the order of N*EPSLN*||A||.
1343  temp1 = sqrt(dble(n))*epsln
1344  DO 3001 p = 2, n
1345  IF ( abs(a(p,p)) .GE. (temp1*abs(a(1,1))) ) THEN
1346  nr = nr + 1
1347  ELSE
1348  GO TO 3002
1349  END IF
1350  3001 CONTINUE
1351  3002 CONTINUE
1352  ELSE IF ( l2rank ) THEN
1353 * .. similarly as above, only slightly more gentle (less aggressive).
1354 * Sudden drop on the diagonal of R1 is used as the criterion for
1355 * close-to-rank-deficient.
1356  temp1 = sqrt(sfmin)
1357  DO 3401 p = 2, n
1358  IF ( ( abs(a(p,p)) .LT. (epsln*abs(a(p-1,p-1))) ) .OR.
1359  $ ( abs(a(p,p)) .LT. small ) .OR.
1360  $ ( l2kill .AND. (abs(a(p,p)) .LT. temp1) ) ) GO TO 3402
1361  nr = nr + 1
1362  3401 CONTINUE
1363  3402 CONTINUE
1364 *
1365  ELSE
1366 * The goal is high relative accuracy. However, if the matrix
1367 * has high scaled condition number the relative accuracy is in
1368 * general not feasible. Later on, a condition number estimator
1369 * will be deployed to estimate the scaled condition number.
1370 * Here we just remove the underflowed part of the triangular
1371 * factor. This prevents the situation in which the code is
1372 * working hard to get the accuracy not warranted by the data.
1373  temp1 = sqrt(sfmin)
1374  DO 3301 p = 2, n
1375  IF ( ( abs(a(p,p)) .LT. small ) .OR.
1376  $ ( l2kill .AND. (abs(a(p,p)) .LT. temp1) ) ) GO TO 3302
1377  nr = nr + 1
1378  3301 CONTINUE
1379  3302 CONTINUE
1380 *
1381  END IF
1382 *
1383  almort = .false.
1384  IF ( nr .EQ. n ) THEN
1385  maxprj = one
1386  DO 3051 p = 2, n
1387  temp1 = abs(a(p,p)) / sva(iwork(p))
1388  maxprj = min( maxprj, temp1 )
1389  3051 CONTINUE
1390  IF ( maxprj**2 .GE. one - dble(n)*epsln ) almort = .true.
1391  END IF
1392 *
1393 *
1394  sconda = - one
1395  condr1 = - one
1396  condr2 = - one
1397 *
1398  IF ( errest ) THEN
1399  IF ( n .EQ. nr ) THEN
1400  IF ( rsvec ) THEN
1401 * .. V is available as workspace
1402  CALL zlacpy( 'U', n, n, a, lda, v, ldv )
1403  DO 3053 p = 1, n
1404  temp1 = sva(iwork(p))
1405  CALL zdscal( p, one/temp1, v(1,p), 1 )
1406  3053 CONTINUE
1407  IF ( lsvec )THEN
1408  CALL zpocon( 'U', n, v, ldv, one, temp1,
1409  $ cwork(n+1), rwork, ierr )
1410  ELSE
1411  CALL zpocon( 'U', n, v, ldv, one, temp1,
1412  $ cwork, rwork, ierr )
1413  END IF
1414 *
1415  ELSE IF ( lsvec ) THEN
1416 * .. U is available as workspace
1417  CALL zlacpy( 'U', n, n, a, lda, u, ldu )
1418  DO 3054 p = 1, n
1419  temp1 = sva(iwork(p))
1420  CALL zdscal( p, one/temp1, u(1,p), 1 )
1421  3054 CONTINUE
1422  CALL zpocon( 'U', n, u, ldu, one, temp1,
1423  $ cwork(n+1), rwork, ierr )
1424  ELSE
1425  CALL zlacpy( 'U', n, n, a, lda, cwork, n )
1426 *[] CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N )
1427 * Change: here index shifted by N to the left, CWORK(1:N)
1428 * not needed for SIGMA only computation
1429  DO 3052 p = 1, n
1430  temp1 = sva(iwork(p))
1431 *[] CALL ZDSCAL( p, ONE/TEMP1, CWORK(N+(p-1)*N+1), 1 )
1432  CALL zdscal( p, one/temp1, cwork((p-1)*n+1), 1 )
1433  3052 CONTINUE
1434 * .. the columns of R are scaled to have unit Euclidean lengths.
1435 *[] CALL ZPOCON( 'U', N, CWORK(N+1), N, ONE, TEMP1,
1436 *[] $ CWORK(N+N*N+1), RWORK, IERR )
1437  CALL zpocon( 'U', n, cwork, n, one, temp1,
1438  $ cwork(n*n+1), rwork, ierr )
1439 *
1440  END IF
1441  IF ( temp1 .NE. zero ) THEN
1442  sconda = one / sqrt(temp1)
1443  ELSE
1444  sconda = - one
1445  END IF
1446 * SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1).
1447 * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
1448  ELSE
1449  sconda = - one
1450  END IF
1451  END IF
1452 *
1453  l2pert = l2pert .AND. ( abs( a(1,1)/a(nr,nr) ) .GT. sqrt(big1) )
1454 * If there is no violent scaling, artificial perturbation is not needed.
1455 *
1456 * Phase 3:
1457 *
1458  IF ( .NOT. ( rsvec .OR. lsvec ) ) THEN
1459 *
1460 * Singular Values only
1461 *
1462 * .. transpose A(1:NR,1:N)
1463  DO 1946 p = 1, min( n-1, nr )
1464  CALL zcopy( n-p, a(p,p+1), lda, a(p+1,p), 1 )
1465  CALL zlacgv( n-p+1, a(p,p), 1 )
1466  1946 CONTINUE
1467  IF ( nr .EQ. n ) a(n,n) = conjg(a(n,n))
1468 *
1469 * The following two DO-loops introduce small relative perturbation
1470 * into the strict upper triangle of the lower triangular matrix.
1471 * Small entries below the main diagonal are also changed.
1472 * This modification is useful if the computing environment does not
1473 * provide/allow FLUSH TO ZERO underflow, for it prevents many
1474 * annoying denormalized numbers in case of strongly scaled matrices.
1475 * The perturbation is structured so that it does not introduce any
1476 * new perturbation of the singular values, and it does not destroy
1477 * the job done by the preconditioner.
1478 * The licence for this perturbation is in the variable L2PERT, which
1479 * should be .FALSE. if FLUSH TO ZERO underflow is active.
1480 *
1481  IF ( .NOT. almort ) THEN
1482 *
1483  IF ( l2pert ) THEN
1484 * XSC = SQRT(SMALL)
1485  xsc = epsln / dble(n)
1486  DO 4947 q = 1, nr
1487  ctemp = dcmplx(xsc*abs(a(q,q)),zero)
1488  DO 4949 p = 1, n
1489  IF ( ( (p.GT.q) .AND. (abs(a(p,q)).LE.temp1) )
1490  $ .OR. ( p .LT. q ) )
1491 * $ A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) )
1492  $ a(p,q) = ctemp
1493  4949 CONTINUE
1494  4947 CONTINUE
1495  ELSE
1496  CALL zlaset( 'U', nr-1,nr-1, czero,czero, a(1,2),lda )
1497  END IF
1498 *
1499 * .. second preconditioning using the QR factorization
1500 *
1501  CALL zgeqrf( n,nr, a,lda, cwork, cwork(n+1),lwork-n, ierr )
1502 *
1503 * .. and transpose upper to lower triangular
1504  DO 1948 p = 1, nr - 1
1505  CALL zcopy( nr-p, a(p,p+1), lda, a(p+1,p), 1 )
1506  CALL zlacgv( nr-p+1, a(p,p), 1 )
1507  1948 CONTINUE
1508 *
1509  END IF
1510 *
1511 * Row-cyclic Jacobi SVD algorithm with column pivoting
1512 *
1513 * .. again some perturbation (a "background noise") is added
1514 * to drown denormals
1515  IF ( l2pert ) THEN
1516 * XSC = SQRT(SMALL)
1517  xsc = epsln / dble(n)
1518  DO 1947 q = 1, nr
1519  ctemp = dcmplx(xsc*abs(a(q,q)),zero)
1520  DO 1949 p = 1, nr
1521  IF ( ( (p.GT.q) .AND. (abs(a(p,q)).LE.temp1) )
1522  $ .OR. ( p .LT. q ) )
1523 * $ A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) )
1524  $ a(p,q) = ctemp
1525  1949 CONTINUE
1526  1947 CONTINUE
1527  ELSE
1528  CALL zlaset( 'U', nr-1, nr-1, czero, czero, a(1,2), lda )
1529  END IF
1530 *
1531 * .. and one-sided Jacobi rotations are started on a lower
1532 * triangular matrix (plus perturbation which is ignored in
1533 * the part which destroys triangular form (confusing?!))
1534 *
1535  CALL zgesvj( 'L', 'N', 'N', nr, nr, a, lda, sva,
1536  $ n, v, ldv, cwork, lwork, rwork, lrwork, info )
1537 *
1538  scalem = rwork(1)
1539  numrank = nint(rwork(2))
1540 *
1541 *
1542  ELSE IF ( ( rsvec .AND. ( .NOT. lsvec ) .AND. ( .NOT. jracc ) )
1543  $ .OR.
1544  $ ( jracc .AND. ( .NOT. lsvec ) .AND. ( nr .NE. n ) ) ) THEN
1545 *
1546 * -> Singular Values and Right Singular Vectors <-
1547 *
1548  IF ( almort ) THEN
1549 *
1550 * .. in this case NR equals N
1551  DO 1998 p = 1, nr
1552  CALL zcopy( n-p+1, a(p,p), lda, v(p,p), 1 )
1553  CALL zlacgv( n-p+1, v(p,p), 1 )
1554  1998 CONTINUE
1555  CALL zlaset( 'U', nr-1,nr-1, czero, czero, v(1,2), ldv )
1556 *
1557  CALL zgesvj( 'L','U','N', n, nr, v, ldv, sva, nr, a, lda,
1558  $ cwork, lwork, rwork, lrwork, info )
1559  scalem = rwork(1)
1560  numrank = nint(rwork(2))
1561 
1562  ELSE
1563 *
1564 * .. two more QR factorizations ( one QRF is not enough, two require
1565 * accumulated product of Jacobi rotations, three are perfect )
1566 *
1567  CALL zlaset( 'L', nr-1,nr-1, czero, czero, a(2,1), lda )
1568  CALL zgelqf( nr,n, a, lda, cwork, cwork(n+1), lwork-n, ierr)
1569  CALL zlacpy( 'L', nr, nr, a, lda, v, ldv )
1570  CALL zlaset( 'U', nr-1,nr-1, czero, czero, v(1,2), ldv )
1571  CALL zgeqrf( nr, nr, v, ldv, cwork(n+1), cwork(2*n+1),
1572  $ lwork-2*n, ierr )
1573  DO 8998 p = 1, nr
1574  CALL zcopy( nr-p+1, v(p,p), ldv, v(p,p), 1 )
1575  CALL zlacgv( nr-p+1, v(p,p), 1 )
1576  8998 CONTINUE
1577  CALL zlaset('U', nr-1, nr-1, czero, czero, v(1,2), ldv)
1578 *
1579  CALL zgesvj( 'L', 'U','N', nr, nr, v,ldv, sva, nr, u,
1580  $ ldu, cwork(n+1), lwork-n, rwork, lrwork, info )
1581  scalem = rwork(1)
1582  numrank = nint(rwork(2))
1583  IF ( nr .LT. n ) THEN
1584  CALL zlaset( 'A',n-nr, nr, czero,czero, v(nr+1,1), ldv )
1585  CALL zlaset( 'A',nr, n-nr, czero,czero, v(1,nr+1), ldv )
1586  CALL zlaset( 'A',n-nr,n-nr,czero,cone, v(nr+1,nr+1),ldv )
1587  END IF
1588 *
1589  CALL zunmlq( 'L', 'C', n, n, nr, a, lda, cwork,
1590  $ v, ldv, cwork(n+1), lwork-n, ierr )
1591 *
1592  END IF
1593 * .. permute the rows of V
1594 * DO 8991 p = 1, N
1595 * CALL ZCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA )
1596 * 8991 CONTINUE
1597 * CALL ZLACPY( 'All', N, N, A, LDA, V, LDV )
1598  CALL zlapmr( .false., n, n, v, ldv, iwork )
1599 *
1600  IF ( transp ) THEN
1601  CALL zlacpy( 'A', n, n, v, ldv, u, ldu )
1602  END IF
1603 *
1604  ELSE IF ( jracc .AND. (.NOT. lsvec) .AND. ( nr.EQ. n ) ) THEN
1605 *
1606  CALL zlaset( 'L', n-1,n-1, czero, czero, a(2,1), lda )
1607 *
1608  CALL zgesvj( 'U','N','V', n, n, a, lda, sva, n, v, ldv,
1609  $ cwork, lwork, rwork, lrwork, info )
1610  scalem = rwork(1)
1611  numrank = nint(rwork(2))
1612  CALL zlapmr( .false., n, n, v, ldv, iwork )
1613 *
1614  ELSE IF ( lsvec .AND. ( .NOT. rsvec ) ) THEN
1615 *
1616 * .. Singular Values and Left Singular Vectors ..
1617 *
1618 * .. second preconditioning step to avoid need to accumulate
1619 * Jacobi rotations in the Jacobi iterations.
1620  DO 1965 p = 1, nr
1621  CALL zcopy( n-p+1, a(p,p), lda, u(p,p), 1 )
1622  CALL zlacgv( n-p+1, u(p,p), 1 )
1623  1965 CONTINUE
1624  CALL zlaset( 'U', nr-1, nr-1, czero, czero, u(1,2), ldu )
1625 *
1626  CALL zgeqrf( n, nr, u, ldu, cwork(n+1), cwork(2*n+1),
1627  $ lwork-2*n, ierr )
1628 *
1629  DO 1967 p = 1, nr - 1
1630  CALL zcopy( nr-p, u(p,p+1), ldu, u(p+1,p), 1 )
1631  CALL zlacgv( n-p+1, u(p,p), 1 )
1632  1967 CONTINUE
1633  CALL zlaset( 'U', nr-1, nr-1, czero, czero, u(1,2), ldu )
1634 *
1635  CALL zgesvj( 'L', 'U', 'N', nr,nr, u, ldu, sva, nr, a,
1636  $ lda, cwork(n+1), lwork-n, rwork, lrwork, info )
1637  scalem = rwork(1)
1638  numrank = nint(rwork(2))
1639 *
1640  IF ( nr .LT. m ) THEN
1641  CALL zlaset( 'A', m-nr, nr,czero, czero, u(nr+1,1), ldu )
1642  IF ( nr .LT. n1 ) THEN
1643  CALL zlaset( 'A',nr, n1-nr, czero, czero, u(1,nr+1),ldu )
1644  CALL zlaset( 'A',m-nr,n1-nr,czero,cone,u(nr+1,nr+1),ldu )
1645  END IF
1646  END IF
1647 *
1648  CALL zunmqr( 'L', 'N', m, n1, n, a, lda, cwork, u,
1649  $ ldu, cwork(n+1), lwork-n, ierr )
1650 *
1651  IF ( rowpiv )
1652  $ CALL zlaswp( n1, u, ldu, 1, m-1, iwork(iwoff+1), -1 )
1653 *
1654  DO 1974 p = 1, n1
1655  xsc = one / dznrm2( m, u(1,p), 1 )
1656  CALL zdscal( m, xsc, u(1,p), 1 )
1657  1974 CONTINUE
1658 *
1659  IF ( transp ) THEN
1660  CALL zlacpy( 'A', n, n, u, ldu, v, ldv )
1661  END IF
1662 *
1663  ELSE
1664 *
1665 * .. Full SVD ..
1666 *
1667  IF ( .NOT. jracc ) THEN
1668 *
1669  IF ( .NOT. almort ) THEN
1670 *
1671 * Second Preconditioning Step (QRF [with pivoting])
1672 * Note that the composition of TRANSPOSE, QRF and TRANSPOSE is
1673 * equivalent to an LQF CALL. Since in many libraries the QRF
1674 * seems to be better optimized than the LQF, we do explicit
1675 * transpose and use the QRF. This is subject to changes in an
1676 * optimized implementation of ZGEJSV.
1677 *
1678  DO 1968 p = 1, nr
1679  CALL zcopy( n-p+1, a(p,p), lda, v(p,p), 1 )
1680  CALL zlacgv( n-p+1, v(p,p), 1 )
1681  1968 CONTINUE
1682 *
1683 * .. the following two loops perturb small entries to avoid
1684 * denormals in the second QR factorization, where they are
1685 * as good as zeros. This is done to avoid painfully slow
1686 * computation with denormals. The relative size of the perturbation
1687 * is a parameter that can be changed by the implementer.
1688 * This perturbation device will be obsolete on machines with
1689 * properly implemented arithmetic.
1690 * To switch it off, set L2PERT=.FALSE. To remove it from the
1691 * code, remove the action under L2PERT=.TRUE., leave the ELSE part.
1692 * The following two loops should be blocked and fused with the
1693 * transposed copy above.
1694 *
1695  IF ( l2pert ) THEN
1696  xsc = sqrt(small)
1697  DO 2969 q = 1, nr
1698  ctemp = dcmplx(xsc*abs( v(q,q) ),zero)
1699  DO 2968 p = 1, n
1700  IF ( ( p .GT. q ) .AND. ( abs(v(p,q)) .LE. temp1 )
1701  $ .OR. ( p .LT. q ) )
1702 * $ V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) )
1703  $ v(p,q) = ctemp
1704  IF ( p .LT. q ) v(p,q) = - v(p,q)
1705  2968 CONTINUE
1706  2969 CONTINUE
1707  ELSE
1708  CALL zlaset( 'U', nr-1, nr-1, czero, czero, v(1,2), ldv )
1709  END IF
1710 *
1711 * Estimate the row scaled condition number of R1
1712 * (If R1 is rectangular, N > NR, then the condition number
1713 * of the leading NR x NR submatrix is estimated.)
1714 *
1715  CALL zlacpy( 'L', nr, nr, v, ldv, cwork(2*n+1), nr )
1716  DO 3950 p = 1, nr
1717  temp1 = dznrm2(nr-p+1,cwork(2*n+(p-1)*nr+p),1)
1718  CALL zdscal(nr-p+1,one/temp1,cwork(2*n+(p-1)*nr+p),1)
1719  3950 CONTINUE
1720  CALL zpocon('L',nr,cwork(2*n+1),nr,one,temp1,
1721  $ cwork(2*n+nr*nr+1),rwork,ierr)
1722  condr1 = one / sqrt(temp1)
1723 * .. here need a second opinion on the condition number
1724 * .. then assume worst case scenario
1725 * R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
1726 * more conservative <=> CONDR1 .LT. SQRT(DBLE(N))
1727 *
1728  cond_ok = sqrt(sqrt(dble(nr)))
1729 *[TP] COND_OK is a tuning parameter.
1730 *
1731  IF ( condr1 .LT. cond_ok ) THEN
1732 * .. the second QRF without pivoting. Note: in an optimized
1733 * implementation, this QRF should be implemented as the QRF
1734 * of a lower triangular matrix.
1735 * R1^* = Q2 * R2
1736  CALL zgeqrf( n, nr, v, ldv, cwork(n+1), cwork(2*n+1),
1737  $ lwork-2*n, ierr )
1738 *
1739  IF ( l2pert ) THEN
1740  xsc = sqrt(small)/epsln
1741  DO 3959 p = 2, nr
1742  DO 3958 q = 1, p - 1
1743  ctemp=dcmplx(xsc*min(abs(v(p,p)),abs(v(q,q))),
1744  $ zero)
1745  IF ( abs(v(q,p)) .LE. temp1 )
1746 * $ V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) )
1747  $ v(q,p) = ctemp
1748  3958 CONTINUE
1749  3959 CONTINUE
1750  END IF
1751 *
1752  IF ( nr .NE. n )
1753  $ CALL zlacpy( 'A', n, nr, v, ldv, cwork(2*n+1), n )
1754 * .. save ...
1755 *
1756 * .. this transposed copy should be better than naive
1757  DO 1969 p = 1, nr - 1
1758  CALL zcopy( nr-p, v(p,p+1), ldv, v(p+1,p), 1 )
1759  CALL zlacgv(nr-p+1, v(p,p), 1 )
1760  1969 CONTINUE
1761  v(nr,nr)=conjg(v(nr,nr))
1762 *
1763  condr2 = condr1
1764 *
1765  ELSE
1766 *
1767 * .. ill-conditioned case: second QRF with pivoting
1768 * Note that windowed pivoting would be equally good
1769 * numerically, and more run-time efficient. So, in
1770 * an optimal implementation, the next call to ZGEQP3
1771 * should be replaced with eg. CALL ZGEQPX (ACM TOMS #782)
1772 * with properly (carefully) chosen parameters.
1773 *
1774 * R1^* * P2 = Q2 * R2
1775  DO 3003 p = 1, nr
1776  iwork(n+p) = 0
1777  3003 CONTINUE
1778  CALL zgeqp3( n, nr, v, ldv, iwork(n+1), cwork(n+1),
1779  $ cwork(2*n+1), lwork-2*n, rwork, ierr )
1780 ** CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
1781 ** $ LWORK-2*N, IERR )
1782  IF ( l2pert ) THEN
1783  xsc = sqrt(small)
1784  DO 3969 p = 2, nr
1785  DO 3968 q = 1, p - 1
1786  ctemp=dcmplx(xsc*min(abs(v(p,p)),abs(v(q,q))),
1787  $ zero)
1788  IF ( abs(v(q,p)) .LE. temp1 )
1789 * $ V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) )
1790  $ v(q,p) = ctemp
1791  3968 CONTINUE
1792  3969 CONTINUE
1793  END IF
1794 *
1795  CALL zlacpy( 'A', n, nr, v, ldv, cwork(2*n+1), n )
1796 *
1797  IF ( l2pert ) THEN
1798  xsc = sqrt(small)
1799  DO 8970 p = 2, nr
1800  DO 8971 q = 1, p - 1
1801  ctemp=dcmplx(xsc*min(abs(v(p,p)),abs(v(q,q))),
1802  $ zero)
1803 * V(p,q) = - TEMP1*( V(q,p) / ABS(V(q,p)) )
1804  v(p,q) = - ctemp
1805  8971 CONTINUE
1806  8970 CONTINUE
1807  ELSE
1808  CALL zlaset( 'L',nr-1,nr-1,czero,czero,v(2,1),ldv )
1809  END IF
1810 * Now, compute R2 = L3 * Q3, the LQ factorization.
1811  CALL zgelqf( nr, nr, v, ldv, cwork(2*n+n*nr+1),
1812  $ cwork(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr, ierr )
1813 * .. and estimate the condition number
1814  CALL zlacpy( 'L',nr,nr,v,ldv,cwork(2*n+n*nr+nr+1),nr )
1815  DO 4950 p = 1, nr
1816  temp1 = dznrm2( p, cwork(2*n+n*nr+nr+p), nr )
1817  CALL zdscal( p, one/temp1, cwork(2*n+n*nr+nr+p), nr )
1818  4950 CONTINUE
1819  CALL zpocon( 'L',nr,cwork(2*n+n*nr+nr+1),nr,one,temp1,
1820  $ cwork(2*n+n*nr+nr+nr*nr+1),rwork,ierr )
1821  condr2 = one / sqrt(temp1)
1822 *
1823 *
1824  IF ( condr2 .GE. cond_ok ) THEN
1825 * .. save the Householder vectors used for Q3
1826 * (this overwrites the copy of R2, as it will not be
1827 * needed in this branch, but it does not overwritte the
1828 * Huseholder vectors of Q2.).
1829  CALL zlacpy( 'U', nr, nr, v, ldv, cwork(2*n+1), n )
1830 * .. and the rest of the information on Q3 is in
1831 * WORK(2*N+N*NR+1:2*N+N*NR+N)
1832  END IF
1833 *
1834  END IF
1835 *
1836  IF ( l2pert ) THEN
1837  xsc = sqrt(small)
1838  DO 4968 q = 2, nr
1839  ctemp = xsc * v(q,q)
1840  DO 4969 p = 1, q - 1
1841 * V(p,q) = - TEMP1*( V(p,q) / ABS(V(p,q)) )
1842  v(p,q) = - ctemp
1843  4969 CONTINUE
1844  4968 CONTINUE
1845  ELSE
1846  CALL zlaset( 'U', nr-1,nr-1, czero,czero, v(1,2), ldv )
1847  END IF
1848 *
1849 * Second preconditioning finished; continue with Jacobi SVD
1850 * The input matrix is lower trinagular.
1851 *
1852 * Recover the right singular vectors as solution of a well
1853 * conditioned triangular matrix equation.
1854 *
1855  IF ( condr1 .LT. cond_ok ) THEN
1856 *
1857  CALL zgesvj( 'L','U','N',nr,nr,v,ldv,sva,nr,u, ldu,
1858  $ cwork(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,rwork,
1859  $ lrwork, info )
1860  scalem = rwork(1)
1861  numrank = nint(rwork(2))
1862  DO 3970 p = 1, nr
1863  CALL zcopy( nr, v(1,p), 1, u(1,p), 1 )
1864  CALL zdscal( nr, sva(p), v(1,p), 1 )
1865  3970 CONTINUE
1866 
1867 * .. pick the right matrix equation and solve it
1868 *
1869  IF ( nr .EQ. n ) THEN
1870 * :)) .. best case, R1 is inverted. The solution of this matrix
1871 * equation is Q2*V2 = the product of the Jacobi rotations
1872 * used in ZGESVJ, premultiplied with the orthogonal matrix
1873 * from the second QR factorization.
1874  CALL ztrsm('L','U','N','N', nr,nr,cone, a,lda, v,ldv)
1875  ELSE
1876 * .. R1 is well conditioned, but non-square. Adjoint of R2
1877 * is inverted to get the product of the Jacobi rotations
1878 * used in ZGESVJ. The Q-factor from the second QR
1879 * factorization is then built in explicitly.
1880  CALL ztrsm('L','U','C','N',nr,nr,cone,cwork(2*n+1),
1881  $ n,v,ldv)
1882  IF ( nr .LT. n ) THEN
1883  CALL zlaset('A',n-nr,nr,czero,czero,v(nr+1,1),ldv)
1884  CALL zlaset('A',nr,n-nr,czero,czero,v(1,nr+1),ldv)
1885  CALL zlaset('A',n-nr,n-nr,czero,cone,v(nr+1,nr+1),ldv)
1886  END IF
1887  CALL zunmqr('L','N',n,n,nr,cwork(2*n+1),n,cwork(n+1),
1888  $ v,ldv,cwork(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr)
1889  END IF
1890 *
1891  ELSE IF ( condr2 .LT. cond_ok ) THEN
1892 *
1893 * The matrix R2 is inverted. The solution of the matrix equation
1894 * is Q3^* * V3 = the product of the Jacobi rotations (appplied to
1895 * the lower triangular L3 from the LQ factorization of
1896 * R2=L3*Q3), pre-multiplied with the transposed Q3.
1897  CALL zgesvj( 'L', 'U', 'N', nr, nr, v, ldv, sva, nr, u,
1898  $ ldu, cwork(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr,
1899  $ rwork, lrwork, info )
1900  scalem = rwork(1)
1901  numrank = nint(rwork(2))
1902  DO 3870 p = 1, nr
1903  CALL zcopy( nr, v(1,p), 1, u(1,p), 1 )
1904  CALL zdscal( nr, sva(p), u(1,p), 1 )
1905  3870 CONTINUE
1906  CALL ztrsm('L','U','N','N',nr,nr,cone,cwork(2*n+1),n,
1907  $ u,ldu)
1908 * .. apply the permutation from the second QR factorization
1909  DO 873 q = 1, nr
1910  DO 872 p = 1, nr
1911  cwork(2*n+n*nr+nr+iwork(n+p)) = u(p,q)
1912  872 CONTINUE
1913  DO 874 p = 1, nr
1914  u(p,q) = cwork(2*n+n*nr+nr+p)
1915  874 CONTINUE
1916  873 CONTINUE
1917  IF ( nr .LT. n ) THEN
1918  CALL zlaset( 'A',n-nr,nr,czero,czero,v(nr+1,1),ldv )
1919  CALL zlaset( 'A',nr,n-nr,czero,czero,v(1,nr+1),ldv )
1920  CALL zlaset('A',n-nr,n-nr,czero,cone,v(nr+1,nr+1),ldv)
1921  END IF
1922  CALL zunmqr( 'L','N',n,n,nr,cwork(2*n+1),n,cwork(n+1),
1923  $ v,ldv,cwork(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
1924  ELSE
1925 * Last line of defense.
1926 * #:( This is a rather pathological case: no scaled condition
1927 * improvement after two pivoted QR factorizations. Other
1928 * possibility is that the rank revealing QR factorization
1929 * or the condition estimator has failed, or the COND_OK
1930 * is set very close to ONE (which is unnecessary). Normally,
1931 * this branch should never be executed, but in rare cases of
1932 * failure of the RRQR or condition estimator, the last line of
1933 * defense ensures that ZGEJSV completes the task.
1934 * Compute the full SVD of L3 using ZGESVJ with explicit
1935 * accumulation of Jacobi rotations.
1936  CALL zgesvj( 'L', 'U', 'V', nr, nr, v, ldv, sva, nr, u,
1937  $ ldu, cwork(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr,
1938  $ rwork, lrwork, info )
1939  scalem = rwork(1)
1940  numrank = nint(rwork(2))
1941  IF ( nr .LT. n ) THEN
1942  CALL zlaset( 'A',n-nr,nr,czero,czero,v(nr+1,1),ldv )
1943  CALL zlaset( 'A',nr,n-nr,czero,czero,v(1,nr+1),ldv )
1944  CALL zlaset('A',n-nr,n-nr,czero,cone,v(nr+1,nr+1),ldv)
1945  END IF
1946  CALL zunmqr( 'L','N',n,n,nr,cwork(2*n+1),n,cwork(n+1),
1947  $ v,ldv,cwork(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
1948 *
1949  CALL zunmlq( 'L', 'C', nr, nr, nr, cwork(2*n+1), n,
1950  $ cwork(2*n+n*nr+1), u, ldu, cwork(2*n+n*nr+nr+1),
1951  $ lwork-2*n-n*nr-nr, ierr )
1952  DO 773 q = 1, nr
1953  DO 772 p = 1, nr
1954  cwork(2*n+n*nr+nr+iwork(n+p)) = u(p,q)
1955  772 CONTINUE
1956  DO 774 p = 1, nr
1957  u(p,q) = cwork(2*n+n*nr+nr+p)
1958  774 CONTINUE
1959  773 CONTINUE
1960 *
1961  END IF
1962 *
1963 * Permute the rows of V using the (column) permutation from the
1964 * first QRF. Also, scale the columns to make them unit in
1965 * Euclidean norm. This applies to all cases.
1966 *
1967  temp1 = sqrt(dble(n)) * epsln
1968  DO 1972 q = 1, n
1969  DO 972 p = 1, n
1970  cwork(2*n+n*nr+nr+iwork(p)) = v(p,q)
1971  972 CONTINUE
1972  DO 973 p = 1, n
1973  v(p,q) = cwork(2*n+n*nr+nr+p)
1974  973 CONTINUE
1975  xsc = one / dznrm2( n, v(1,q), 1 )
1976  IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1977  $ CALL zdscal( n, xsc, v(1,q), 1 )
1978  1972 CONTINUE
1979 * At this moment, V contains the right singular vectors of A.
1980 * Next, assemble the left singular vector matrix U (M x N).
1981  IF ( nr .LT. m ) THEN
1982  CALL zlaset('A', m-nr, nr, czero, czero, u(nr+1,1), ldu)
1983  IF ( nr .LT. n1 ) THEN
1984  CALL zlaset('A',nr,n1-nr,czero,czero,u(1,nr+1),ldu)
1985  CALL zlaset('A',m-nr,n1-nr,czero,cone,
1986  $ u(nr+1,nr+1),ldu)
1987  END IF
1988  END IF
1989 *
1990 * The Q matrix from the first QRF is built into the left singular
1991 * matrix U. This applies to all cases.
1992 *
1993  CALL zunmqr( 'L', 'N', m, n1, n, a, lda, cwork, u,
1994  $ ldu, cwork(n+1), lwork-n, ierr )
1995 
1996 * The columns of U are normalized. The cost is O(M*N) flops.
1997  temp1 = sqrt(dble(m)) * epsln
1998  DO 1973 p = 1, nr
1999  xsc = one / dznrm2( m, u(1,p), 1 )
2000  IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
2001  $ CALL zdscal( m, xsc, u(1,p), 1 )
2002  1973 CONTINUE
2003 *
2004 * If the initial QRF is computed with row pivoting, the left
2005 * singular vectors must be adjusted.
2006 *
2007  IF ( rowpiv )
2008  $ CALL zlaswp( n1, u, ldu, 1, m-1, iwork(iwoff+1), -1 )
2009 *
2010  ELSE
2011 *
2012 * .. the initial matrix A has almost orthogonal columns and
2013 * the second QRF is not needed
2014 *
2015  CALL zlacpy( 'U', n, n, a, lda, cwork(n+1), n )
2016  IF ( l2pert ) THEN
2017  xsc = sqrt(small)
2018  DO 5970 p = 2, n
2019  ctemp = xsc * cwork( n + (p-1)*n + p )
2020  DO 5971 q = 1, p - 1
2021 * CWORK(N+(q-1)*N+p)=-TEMP1 * ( CWORK(N+(p-1)*N+q) /
2022 * $ ABS(CWORK(N+(p-1)*N+q)) )
2023  cwork(n+(q-1)*n+p)=-ctemp
2024  5971 CONTINUE
2025  5970 CONTINUE
2026  ELSE
2027  CALL zlaset( 'L',n-1,n-1,czero,czero,cwork(n+2),n )
2028  END IF
2029 *
2030  CALL zgesvj( 'U', 'U', 'N', n, n, cwork(n+1), n, sva,
2031  $ n, u, ldu, cwork(n+n*n+1), lwork-n-n*n, rwork, lrwork,
2032  $ info )
2033 *
2034  scalem = rwork(1)
2035  numrank = nint(rwork(2))
2036  DO 6970 p = 1, n
2037  CALL zcopy( n, cwork(n+(p-1)*n+1), 1, u(1,p), 1 )
2038  CALL zdscal( n, sva(p), cwork(n+(p-1)*n+1), 1 )
2039  6970 CONTINUE
2040 *
2041  CALL ztrsm( 'L', 'U', 'N', 'N', n, n,
2042  $ cone, a, lda, cwork(n+1), n )
2043  DO 6972 p = 1, n
2044  CALL zcopy( n, cwork(n+p), n, v(iwork(p),1), ldv )
2045  6972 CONTINUE
2046  temp1 = sqrt(dble(n))*epsln
2047  DO 6971 p = 1, n
2048  xsc = one / dznrm2( n, v(1,p), 1 )
2049  IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
2050  $ CALL zdscal( n, xsc, v(1,p), 1 )
2051  6971 CONTINUE
2052 *
2053 * Assemble the left singular vector matrix U (M x N).
2054 *
2055  IF ( n .LT. m ) THEN
2056  CALL zlaset( 'A', m-n, n, czero, czero, u(n+1,1), ldu )
2057  IF ( n .LT. n1 ) THEN
2058  CALL zlaset('A',n, n1-n, czero, czero, u(1,n+1),ldu)
2059  CALL zlaset( 'A',m-n,n1-n, czero, cone,u(n+1,n+1),ldu)
2060  END IF
2061  END IF
2062  CALL zunmqr( 'L', 'N', m, n1, n, a, lda, cwork, u,
2063  $ ldu, cwork(n+1), lwork-n, ierr )
2064  temp1 = sqrt(dble(m))*epsln
2065  DO 6973 p = 1, n1
2066  xsc = one / dznrm2( m, u(1,p), 1 )
2067  IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
2068  $ CALL zdscal( m, xsc, u(1,p), 1 )
2069  6973 CONTINUE
2070 *
2071  IF ( rowpiv )
2072  $ CALL zlaswp( n1, u, ldu, 1, m-1, iwork(iwoff+1), -1 )
2073 *
2074  END IF
2075 *
2076 * end of the >> almost orthogonal case << in the full SVD
2077 *
2078  ELSE
2079 *
2080 * This branch deploys a preconditioned Jacobi SVD with explicitly
2081 * accumulated rotations. It is included as optional, mainly for
2082 * experimental purposes. It does perform well, and can also be used.
2083 * In this implementation, this branch will be automatically activated
2084 * if the condition number sigma_max(A) / sigma_min(A) is predicted
2085 * to be greater than the overflow threshold. This is because the
2086 * a posteriori computation of the singular vectors assumes robust
2087 * implementation of BLAS and some LAPACK procedures, capable of working
2088 * in presence of extreme values, e.g. when the singular values spread from
2089 * the underflow to the overflow threshold.
2090 *
2091  DO 7968 p = 1, nr
2092  CALL zcopy( n-p+1, a(p,p), lda, v(p,p), 1 )
2093  CALL zlacgv( n-p+1, v(p,p), 1 )
2094  7968 CONTINUE
2095 *
2096  IF ( l2pert ) THEN
2097  xsc = sqrt(small/epsln)
2098  DO 5969 q = 1, nr
2099  ctemp = dcmplx(xsc*abs( v(q,q) ),zero)
2100  DO 5968 p = 1, n
2101  IF ( ( p .GT. q ) .AND. ( abs(v(p,q)) .LE. temp1 )
2102  $ .OR. ( p .LT. q ) )
2103 * $ V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) )
2104  $ v(p,q) = ctemp
2105  IF ( p .LT. q ) v(p,q) = - v(p,q)
2106  5968 CONTINUE
2107  5969 CONTINUE
2108  ELSE
2109  CALL zlaset( 'U', nr-1, nr-1, czero, czero, v(1,2), ldv )
2110  END IF
2111 
2112  CALL zgeqrf( n, nr, v, ldv, cwork(n+1), cwork(2*n+1),
2113  $ lwork-2*n, ierr )
2114  CALL zlacpy( 'L', n, nr, v, ldv, cwork(2*n+1), n )
2115 *
2116  DO 7969 p = 1, nr
2117  CALL zcopy( nr-p+1, v(p,p), ldv, u(p,p), 1 )
2118  CALL zlacgv( nr-p+1, u(p,p), 1 )
2119  7969 CONTINUE
2120 
2121  IF ( l2pert ) THEN
2122  xsc = sqrt(small/epsln)
2123  DO 9970 q = 2, nr
2124  DO 9971 p = 1, q - 1
2125  ctemp = dcmplx(xsc * min(abs(u(p,p)),abs(u(q,q))),
2126  $ zero)
2127 * U(p,q) = - TEMP1 * ( U(q,p) / ABS(U(q,p)) )
2128  u(p,q) = - ctemp
2129  9971 CONTINUE
2130  9970 CONTINUE
2131  ELSE
2132  CALL zlaset('U', nr-1, nr-1, czero, czero, u(1,2), ldu )
2133  END IF
2134 
2135  CALL zgesvj( 'L', 'U', 'V', nr, nr, u, ldu, sva,
2136  $ n, v, ldv, cwork(2*n+n*nr+1), lwork-2*n-n*nr,
2137  $ rwork, lrwork, info )
2138  scalem = rwork(1)
2139  numrank = nint(rwork(2))
2140 
2141  IF ( nr .LT. n ) THEN
2142  CALL zlaset( 'A',n-nr,nr,czero,czero,v(nr+1,1),ldv )
2143  CALL zlaset( 'A',nr,n-nr,czero,czero,v(1,nr+1),ldv )
2144  CALL zlaset( 'A',n-nr,n-nr,czero,cone,v(nr+1,nr+1),ldv )
2145  END IF
2146 
2147  CALL zunmqr( 'L','N',n,n,nr,cwork(2*n+1),n,cwork(n+1),
2148  $ v,ldv,cwork(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
2149 *
2150 * Permute the rows of V using the (column) permutation from the
2151 * first QRF. Also, scale the columns to make them unit in
2152 * Euclidean norm. This applies to all cases.
2153 *
2154  temp1 = sqrt(dble(n)) * epsln
2155  DO 7972 q = 1, n
2156  DO 8972 p = 1, n
2157  cwork(2*n+n*nr+nr+iwork(p)) = v(p,q)
2158  8972 CONTINUE
2159  DO 8973 p = 1, n
2160  v(p,q) = cwork(2*n+n*nr+nr+p)
2161  8973 CONTINUE
2162  xsc = one / dznrm2( n, v(1,q), 1 )
2163  IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
2164  $ CALL zdscal( n, xsc, v(1,q), 1 )
2165  7972 CONTINUE
2166 *
2167 * At this moment, V contains the right singular vectors of A.
2168 * Next, assemble the left singular vector matrix U (M x N).
2169 *
2170  IF ( nr .LT. m ) THEN
2171  CALL zlaset( 'A', m-nr, nr, czero, czero, u(nr+1,1), ldu )
2172  IF ( nr .LT. n1 ) THEN
2173  CALL zlaset('A',nr, n1-nr, czero, czero, u(1,nr+1),ldu)
2174  CALL zlaset('A',m-nr,n1-nr, czero, cone,u(nr+1,nr+1),ldu)
2175  END IF
2176  END IF
2177 *
2178  CALL zunmqr( 'L', 'N', m, n1, n, a, lda, cwork, u,
2179  $ ldu, cwork(n+1), lwork-n, ierr )
2180 *
2181  IF ( rowpiv )
2182  $ CALL zlaswp( n1, u, ldu, 1, m-1, iwork(iwoff+1), -1 )
2183 *
2184 *
2185  END IF
2186  IF ( transp ) THEN
2187 * .. swap U and V because the procedure worked on A^*
2188  DO 6974 p = 1, n
2189  CALL zswap( n, u(1,p), 1, v(1,p), 1 )
2190  6974 CONTINUE
2191  END IF
2192 *
2193  END IF
2194 * end of the full SVD
2195 *
2196 * Undo scaling, if necessary (and possible)
2197 *
2198  IF ( uscal2 .LE. (big/sva(1))*uscal1 ) THEN
2199  CALL dlascl( 'G', 0, 0, uscal1, uscal2, nr, 1, sva, n, ierr )
2200  uscal1 = one
2201  uscal2 = one
2202  END IF
2203 *
2204  IF ( nr .LT. n ) THEN
2205  DO 3004 p = nr+1, n
2206  sva(p) = zero
2207  3004 CONTINUE
2208  END IF
2209 *
2210  rwork(1) = uscal2 * scalem
2211  rwork(2) = uscal1
2212  IF ( errest ) rwork(3) = sconda
2213  IF ( lsvec .AND. rsvec ) THEN
2214  rwork(4) = condr1
2215  rwork(5) = condr2
2216  END IF
2217  IF ( l2tran ) THEN
2218  rwork(6) = entra
2219  rwork(7) = entrat
2220  END IF
2221 *
2222  iwork(1) = nr
2223  iwork(2) = numrank
2224  iwork(3) = warning
2225  IF ( transp ) THEN
2226  iwork(4) = 1
2227  ELSE
2228  iwork(4) = -1
2229  END IF
2230 
2231 *
2232  RETURN
2233 * ..
2234 * .. END OF ZGEJSV
2235 * ..
2236  END
2237 *
zlacgv
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:76
zlapmr
subroutine zlapmr(FORWRD, M, N, X, LDX, K)
ZLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: zlapmr.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
ztrsm
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:182
zgelqf
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
Definition: zgelqf.f:145
zlassq
subroutine zlassq(N, X, INCX, SCALE, SUMSQ)
ZLASSQ updates a sum of squares represented in scaled form.
Definition: zlassq.f:108
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
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:83
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
zungqr
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
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
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:83
zgejsv
subroutine zgejsv(JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, M, N, A, LDA, SVA, U, LDU, V, LDV, CWORK, LWORK, RWORK, LRWORK, IWORK, INFO)
ZGEJSV
Definition: zgejsv.f:571
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
dlassq
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f:105
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
zgesvj
subroutine zgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, CWORK, LWORK, RWORK, LRWORK, INFO)
ZGESVJ
Definition: zgesvj.f:353
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81