Annotation of rpl/lapack/lapack/zgejsv.f, revision 1.9
1.6 bertrand 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.
1.8 bertrand 83: *> = 'G': Computation as with 'F' with an additional estimate of the
1.6 bertrand 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
1.8 bertrand 89: *> numerically rank deficient. The error in the computed
1.6 bertrand 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
1.8 bertrand 120: *> computed as the product of Jacobi rotations, if JOBT = 'N'.
1.6 bertrand 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
1.8 bertrand 134: *> SQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN,
1.6 bertrand 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).
1.8 bertrand 232: *> If JOBU = 'W' .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N),
1.6 bertrand 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;
1.8 bertrand 254: *> If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N),
1.6 bertrand 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: *>
1.8 bertrand 285: *> 1. If only SIGMA is needed ( JOBU = 'N', JOBV = 'N' ) and
1.6 bertrand 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)).
1.8 bertrand 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'):
1.6 bertrand 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
1.8 bertrand 321: *> 3.1 .. no scaled condition estimate requested (JOBE = 'N'):
1.6 bertrand 322: *> -> the minimal requirement is LWORK >= 3*N.
323: *> -> For optimal performance:
1.8 bertrand 324: *> if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
1.6 bertrand 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:
1.8 bertrand 332: *> if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
1.6 bertrand 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)).
1.8 bertrand 337: *> 4. If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and
338: *> 4.1. if JOBV = 'V'
1.6 bertrand 339: *> the minimal requirement is LWORK >= 5*N+2*N*N.
1.8 bertrand 340: *> 4.2. if JOBV = 'J' the minimal requirement is
1.6 bertrand 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
1.8 bertrand 359: *> column equilibrated A. (If JOBA = 'E' or 'G')
1.6 bertrand 360: *> SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1).
1.9 ! bertrand 361: *> It is computed using ZPOCON. It holds
1.6 bertrand 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
1.9 ! bertrand 370: *> useful for the analysis of the algorithm. They are provided for
1.6 bertrand 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.
1.8 bertrand 378: *> The following two parameters are computed if JOBT = 'T'.
1.6 bertrand 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:
1.8 bertrand 459: *> If IWORK(3) = 1 then some of the column norms of A
1.6 bertrand 460: *> were denormalized floats. The requested high accuracy
461: *> is not warranted by the data.
1.8 bertrand 462: *> IWORK(4) = 1 or -1. If IWORK(4) = 1, then the procedure used A^* to
1.6 bertrand 463: *> do the job as specified by the JOB parameters.
1.8 bertrand 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
1.6 bertrand 466: *> IWORK for the job parameters used in the call.
467: *> \endverbatim
468: *>
469: *> \param[out] INFO
470: *> \verbatim
471: *> INFO is INTEGER
1.8 bertrand 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.
1.6 bertrand 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: *> \ingroup complex16GEsing
487: *
488: *> \par Further Details:
489: * =====================
490: *>
491: *> \verbatim
492: *>
493: *> ZGEJSV implements a preconditioned Jacobi SVD algorithm. It uses ZGEQP3,
494: *> ZGEQRF, and ZGELQF as preprocessors and preconditioners. Optionally, an
495: *> additional row pivoting can be used as a preprocessor, which in some
496: *> cases results in much higher accuracy. An example is matrix A with the
497: *> structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
498: *> diagonal matrices and C is well-conditioned matrix. In that case, complete
499: *> pivoting in the first QR factorizations provides accuracy dependent on the
500: *> condition number of C, and independent of D1, D2. Such higher accuracy is
501: *> not completely understood theoretically, but it works well in practice.
502: *> Further, if A can be written as A = B*D, with well-conditioned B and some
503: *> diagonal D, then the high accuracy is guaranteed, both theoretically and
504: *> in software, independent of D. For more details see [1], [2].
505: *> The computational range for the singular values can be the full range
506: *> ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
507: *> & LAPACK routines called by ZGEJSV are implemented to work in that range.
508: *> If that is not the case, then the restriction for safe computation with
509: *> the singular values in the range of normalized IEEE numbers is that the
510: *> spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
511: *> overflow. This code (ZGEJSV) is best used in this restricted range,
512: *> meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are
513: *> returned as zeros. See JOBR for details on this.
514: *> Further, this implementation is somewhat slower than the one described
515: *> in [1,2] due to replacement of some non-LAPACK components, and because
516: *> the choice of some tuning parameters in the iterative part (ZGESVJ) is
517: *> left to the implementer on a particular machine.
518: *> The rank revealing QR factorization (in this code: ZGEQP3) should be
519: *> implemented as in [3]. We have a new version of ZGEQP3 under development
520: *> that is more robust than the current one in LAPACK, with a cleaner cut in
521: *> rank deficient cases. It will be available in the SIGMA library [4].
522: *> If M is much larger than N, it is obvious that the initial QRF with
523: *> column pivoting can be preprocessed by the QRF without pivoting. That
524: *> well known trick is not used in ZGEJSV because in some cases heavy row
525: *> weighting can be treated with complete pivoting. The overhead in cases
526: *> M much larger than N is then only due to pivoting, but the benefits in
527: *> terms of accuracy have prevailed. The implementer/user can incorporate
528: *> this extra QRF step easily. The implementer can also improve data movement
529: *> (matrix transpose, matrix copy, matrix transposed copy) - this
530: *> implementation of ZGEJSV uses only the simplest, naive data movement.
531: *> \endverbatim
532: *
533: *> \par Contributor:
534: * ==================
535: *>
536: *> Zlatko Drmac, Department of Mathematics, Faculty of Science,
537: *> University of Zagreb (Zagreb, Croatia); drmac@math.hr
538: *
539: *> \par References:
540: * ================
541: *>
542: *> \verbatim
543: *>
544: *> [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
545: *> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
546: *> LAPACK Working note 169.
547: *> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
548: *> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
549: *> LAPACK Working note 170.
550: *> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
551: *> factorization software - a case study.
552: *> ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
553: *> LAPACK Working note 176.
554: *> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
555: *> QSVD, (H,K)-SVD computations.
556: *> Department of Mathematics, University of Zagreb, 2008, 2016.
557: *> \endverbatim
558: *
559: *> \par Bugs, examples and comments:
560: * =================================
561: *>
562: *> Please report all bugs and send interesting examples and/or comments to
563: *> drmac@math.hr. Thank you.
564: *>
565: * =====================================================================
566: SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
567: $ M, N, A, LDA, SVA, U, LDU, V, LDV,
568: $ CWORK, LWORK, RWORK, LRWORK, IWORK, INFO )
569: *
1.9 ! bertrand 570: * -- LAPACK computational routine --
1.6 bertrand 571: * -- LAPACK is a software package provided by Univ. of Tennessee, --
572: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
573: *
574: * .. Scalar Arguments ..
575: IMPLICIT NONE
576: INTEGER INFO, LDA, LDU, LDV, LWORK, LRWORK, M, N
577: * ..
578: * .. Array Arguments ..
579: COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ),
580: $ CWORK( LWORK )
581: DOUBLE PRECISION SVA( N ), RWORK( LRWORK )
582: INTEGER IWORK( * )
583: CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
584: * ..
585: *
586: * ===========================================================================
587: *
588: * .. Local Parameters ..
589: DOUBLE PRECISION ZERO, ONE
590: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
591: COMPLEX*16 CZERO, CONE
592: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), CONE = ( 1.0D0, 0.0D0 ) )
593: * ..
594: * .. Local Scalars ..
595: COMPLEX*16 CTEMP
596: DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1,
597: $ COND_OK, CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN,
598: $ MAXPRJ, SCALEM, SCONDA, SFMIN, SMALL, TEMP1,
599: $ USCAL1, USCAL2, XSC
600: INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING
601: LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LQUERY,
602: $ LSVEC, L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN, NOSCAL,
603: $ ROWPIV, RSVEC, TRANSP
604: *
605: INTEGER OPTWRK, MINWRK, MINRWRK, MINIWRK
606: INTEGER LWCON, LWLQF, LWQP3, LWQRF, LWUNMLQ, LWUNMQR, LWUNMQRM,
607: $ LWSVDJ, LWSVDJV, LRWQP3, LRWCON, LRWSVDJ, IWOFF
608: INTEGER LWRK_ZGELQF, LWRK_ZGEQP3, LWRK_ZGEQP3N, LWRK_ZGEQRF,
609: $ LWRK_ZGESVJ, LWRK_ZGESVJV, LWRK_ZGESVJU, LWRK_ZUNMLQ,
610: $ LWRK_ZUNMQR, LWRK_ZUNMQRM
611: * ..
612: * .. Local Arrays
613: COMPLEX*16 CDUMMY(1)
614: DOUBLE PRECISION RDUMMY(1)
615: *
616: * .. Intrinsic Functions ..
617: INTRINSIC ABS, DCMPLX, CONJG, DLOG, MAX, MIN, DBLE, NINT, SQRT
618: * ..
619: * .. External Functions ..
620: DOUBLE PRECISION DLAMCH, DZNRM2
621: INTEGER IDAMAX, IZAMAX
622: LOGICAL LSAME
623: EXTERNAL IDAMAX, IZAMAX, LSAME, DLAMCH, DZNRM2
624: * ..
625: * .. External Subroutines ..
626: EXTERNAL DLASSQ, ZCOPY, ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLAPMR,
627: $ ZLASCL, DLASCL, ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ,
628: $ ZUNMQR, ZPOCON, DSCAL, ZDSCAL, ZSWAP, ZTRSM, ZLACGV,
629: $ XERBLA
630: *
631: EXTERNAL ZGESVJ
632: * ..
633: *
634: * Test the input arguments
635: *
636: LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
637: JRACC = LSAME( JOBV, 'J' )
638: RSVEC = LSAME( JOBV, 'V' ) .OR. JRACC
639: ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' )
640: L2RANK = LSAME( JOBA, 'R' )
641: L2ABER = LSAME( JOBA, 'A' )
642: ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' )
643: L2TRAN = LSAME( JOBT, 'T' ) .AND. ( M .EQ. N )
644: L2KILL = LSAME( JOBR, 'R' )
645: DEFR = LSAME( JOBR, 'N' )
646: L2PERT = LSAME( JOBP, 'P' )
647: *
648: LQUERY = ( LWORK .EQ. -1 ) .OR. ( LRWORK .EQ. -1 )
649: *
650: IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR.
651: $ ERREST .OR. LSAME( JOBA, 'C' ) )) THEN
652: INFO = - 1
653: ELSE IF ( .NOT.( LSVEC .OR. LSAME( JOBU, 'N' ) .OR.
654: $ ( LSAME( JOBU, 'W' ) .AND. RSVEC .AND. L2TRAN ) ) ) THEN
655: INFO = - 2
656: ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR.
657: $ ( LSAME( JOBV, 'W' ) .AND. LSVEC .AND. L2TRAN ) ) ) THEN
658: INFO = - 3
659: ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) ) THEN
660: INFO = - 4
661: ELSE IF ( .NOT. ( LSAME(JOBT,'T') .OR. LSAME(JOBT,'N') ) ) THEN
662: INFO = - 5
663: ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN
664: INFO = - 6
665: ELSE IF ( M .LT. 0 ) THEN
666: INFO = - 7
667: ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN
668: INFO = - 8
669: ELSE IF ( LDA .LT. M ) THEN
670: INFO = - 10
671: ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN
672: INFO = - 13
673: ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN
674: INFO = - 15
675: ELSE
676: * #:)
677: INFO = 0
678: END IF
679: *
680: IF ( INFO .EQ. 0 ) THEN
681: * .. compute the minimal and the optimal workspace lengths
682: * [[The expressions for computing the minimal and the optimal
683: * values of LCWORK, LRWORK are written with a lot of redundancy and
684: * can be simplified. However, this verbose form is useful for
685: * maintenance and modifications of the code.]]
686: *
687: * .. minimal workspace length for ZGEQP3 of an M x N matrix,
688: * ZGEQRF of an N x N matrix, ZGELQF of an N x N matrix,
689: * ZUNMLQ for computing N x N matrix, ZUNMQR for computing N x N
690: * matrix, ZUNMQR for computing M x N matrix, respectively.
691: LWQP3 = N+1
692: LWQRF = MAX( 1, N )
693: LWLQF = MAX( 1, N )
694: LWUNMLQ = MAX( 1, N )
695: LWUNMQR = MAX( 1, N )
696: LWUNMQRM = MAX( 1, M )
697: * .. minimal workspace length for ZPOCON of an N x N matrix
698: LWCON = 2 * N
699: * .. minimal workspace length for ZGESVJ of an N x N matrix,
700: * without and with explicit accumulation of Jacobi rotations
701: LWSVDJ = MAX( 2 * N, 1 )
702: LWSVDJV = MAX( 2 * N, 1 )
703: * .. minimal REAL workspace length for ZGEQP3, ZPOCON, ZGESVJ
1.8 bertrand 704: LRWQP3 = 2 * N
1.6 bertrand 705: LRWCON = N
706: LRWSVDJ = N
707: IF ( LQUERY ) THEN
708: CALL ZGEQP3( M, N, A, LDA, IWORK, CDUMMY, CDUMMY, -1,
709: $ RDUMMY, IERR )
1.9 ! bertrand 710: LWRK_ZGEQP3 = INT( CDUMMY(1) )
1.6 bertrand 711: CALL ZGEQRF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )
1.9 ! bertrand 712: LWRK_ZGEQRF = INT( CDUMMY(1) )
1.6 bertrand 713: CALL ZGELQF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )
1.9 ! bertrand 714: LWRK_ZGELQF = INT( CDUMMY(1) )
1.6 bertrand 715: END IF
716: MINWRK = 2
717: OPTWRK = 2
718: MINIWRK = N
719: IF ( .NOT. (LSVEC .OR. RSVEC ) ) THEN
720: * .. minimal and optimal sizes of the complex workspace if
721: * only the singular values are requested
722: IF ( ERREST ) THEN
723: MINWRK = MAX( N+LWQP3, N**2+LWCON, N+LWQRF, LWSVDJ )
724: ELSE
725: MINWRK = MAX( N+LWQP3, N+LWQRF, LWSVDJ )
726: END IF
727: IF ( LQUERY ) THEN
728: CALL ZGESVJ( 'L', 'N', 'N', N, N, A, LDA, SVA, N, V,
729: $ LDV, CDUMMY, -1, RDUMMY, -1, IERR )
1.9 ! bertrand 730: LWRK_ZGESVJ = INT( CDUMMY(1) )
1.6 bertrand 731: IF ( ERREST ) THEN
732: OPTWRK = MAX( N+LWRK_ZGEQP3, N**2+LWCON,
733: $ N+LWRK_ZGEQRF, LWRK_ZGESVJ )
734: ELSE
735: OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWRK_ZGEQRF,
736: $ LWRK_ZGESVJ )
737: END IF
738: END IF
739: IF ( L2TRAN .OR. ROWPIV ) THEN
740: IF ( ERREST ) THEN
741: MINRWRK = MAX( 7, 2*M, LRWQP3, LRWCON, LRWSVDJ )
742: ELSE
743: MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ )
744: END IF
745: ELSE
746: IF ( ERREST ) THEN
747: MINRWRK = MAX( 7, LRWQP3, LRWCON, LRWSVDJ )
748: ELSE
749: MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
750: END IF
751: END IF
752: IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
753: ELSE IF ( RSVEC .AND. (.NOT.LSVEC) ) THEN
754: * .. minimal and optimal sizes of the complex workspace if the
755: * singular values and the right singular vectors are requested
756: IF ( ERREST ) THEN
757: MINWRK = MAX( N+LWQP3, LWCON, LWSVDJ, N+LWLQF,
758: $ 2*N+LWQRF, N+LWSVDJ, N+LWUNMLQ )
759: ELSE
760: MINWRK = MAX( N+LWQP3, LWSVDJ, N+LWLQF, 2*N+LWQRF,
761: $ N+LWSVDJ, N+LWUNMLQ )
762: END IF
763: IF ( LQUERY ) THEN
764: CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,
765: $ LDA, CDUMMY, -1, RDUMMY, -1, IERR )
1.9 ! bertrand 766: LWRK_ZGESVJ = INT( CDUMMY(1) )
1.6 bertrand 767: CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,
768: $ V, LDV, CDUMMY, -1, IERR )
1.9 ! bertrand 769: LWRK_ZUNMLQ = INT( CDUMMY(1) )
1.6 bertrand 770: IF ( ERREST ) THEN
771: OPTWRK = MAX( N+LWRK_ZGEQP3, LWCON, LWRK_ZGESVJ,
772: $ N+LWRK_ZGELQF, 2*N+LWRK_ZGEQRF,
773: $ N+LWRK_ZGESVJ, N+LWRK_ZUNMLQ )
774: ELSE
775: OPTWRK = MAX( N+LWRK_ZGEQP3, LWRK_ZGESVJ,N+LWRK_ZGELQF,
776: $ 2*N+LWRK_ZGEQRF, N+LWRK_ZGESVJ,
777: $ N+LWRK_ZUNMLQ )
778: END IF
779: END IF
780: IF ( L2TRAN .OR. ROWPIV ) THEN
781: IF ( ERREST ) THEN
782: MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ, LRWCON )
783: ELSE
784: MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ )
785: END IF
786: ELSE
787: IF ( ERREST ) THEN
788: MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
789: ELSE
790: MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
791: END IF
792: END IF
793: IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
794: ELSE IF ( LSVEC .AND. (.NOT.RSVEC) ) THEN
795: * .. minimal and optimal sizes of the complex workspace if the
796: * singular values and the left singular vectors are requested
797: IF ( ERREST ) THEN
798: MINWRK = N + MAX( LWQP3,LWCON,N+LWQRF,LWSVDJ,LWUNMQRM )
799: ELSE
800: MINWRK = N + MAX( LWQP3, N+LWQRF, LWSVDJ, LWUNMQRM )
801: END IF
802: IF ( LQUERY ) THEN
803: CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,
804: $ LDA, CDUMMY, -1, RDUMMY, -1, IERR )
1.9 ! bertrand 805: LWRK_ZGESVJ = INT( CDUMMY(1) )
1.6 bertrand 806: CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
807: $ LDU, CDUMMY, -1, IERR )
1.9 ! bertrand 808: LWRK_ZUNMQRM = INT( CDUMMY(1) )
1.6 bertrand 809: IF ( ERREST ) THEN
810: OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, N+LWRK_ZGEQRF,
811: $ LWRK_ZGESVJ, LWRK_ZUNMQRM )
812: ELSE
813: OPTWRK = N + MAX( LWRK_ZGEQP3, N+LWRK_ZGEQRF,
814: $ LWRK_ZGESVJ, LWRK_ZUNMQRM )
815: END IF
816: END IF
817: IF ( L2TRAN .OR. ROWPIV ) THEN
818: IF ( ERREST ) THEN
819: MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ, LRWCON )
820: ELSE
821: MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ )
822: END IF
823: ELSE
824: IF ( ERREST ) THEN
825: MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
826: ELSE
827: MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
828: END IF
829: END IF
830: IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
831: ELSE
832: * .. minimal and optimal sizes of the complex workspace if the
833: * full SVD is requested
834: IF ( .NOT. JRACC ) THEN
835: IF ( ERREST ) THEN
836: MINWRK = MAX( N+LWQP3, N+LWCON, 2*N+N**2+LWCON,
837: $ 2*N+LWQRF, 2*N+LWQP3,
838: $ 2*N+N**2+N+LWLQF, 2*N+N**2+N+N**2+LWCON,
839: $ 2*N+N**2+N+LWSVDJ, 2*N+N**2+N+LWSVDJV,
840: $ 2*N+N**2+N+LWUNMQR,2*N+N**2+N+LWUNMLQ,
841: $ N+N**2+LWSVDJ, N+LWUNMQRM )
842: ELSE
843: MINWRK = MAX( N+LWQP3, 2*N+N**2+LWCON,
844: $ 2*N+LWQRF, 2*N+LWQP3,
845: $ 2*N+N**2+N+LWLQF, 2*N+N**2+N+N**2+LWCON,
846: $ 2*N+N**2+N+LWSVDJ, 2*N+N**2+N+LWSVDJV,
847: $ 2*N+N**2+N+LWUNMQR,2*N+N**2+N+LWUNMLQ,
848: $ N+N**2+LWSVDJ, N+LWUNMQRM )
849: END IF
850: MINIWRK = MINIWRK + N
851: IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
852: ELSE
853: IF ( ERREST ) THEN
854: MINWRK = MAX( N+LWQP3, N+LWCON, 2*N+LWQRF,
855: $ 2*N+N**2+LWSVDJV, 2*N+N**2+N+LWUNMQR,
856: $ N+LWUNMQRM )
857: ELSE
858: MINWRK = MAX( N+LWQP3, 2*N+LWQRF,
859: $ 2*N+N**2+LWSVDJV, 2*N+N**2+N+LWUNMQR,
860: $ N+LWUNMQRM )
861: END IF
862: IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
863: END IF
864: IF ( LQUERY ) THEN
865: CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
866: $ LDU, CDUMMY, -1, IERR )
1.9 ! bertrand 867: LWRK_ZUNMQRM = INT( CDUMMY(1) )
1.6 bertrand 868: CALL ZUNMQR( 'L', 'N', N, N, N, A, LDA, CDUMMY, U,
869: $ LDU, CDUMMY, -1, IERR )
1.9 ! bertrand 870: LWRK_ZUNMQR = INT( CDUMMY(1) )
1.6 bertrand 871: IF ( .NOT. JRACC ) THEN
872: CALL ZGEQP3( N,N, A, LDA, IWORK, CDUMMY,CDUMMY, -1,
873: $ RDUMMY, IERR )
1.9 ! bertrand 874: LWRK_ZGEQP3N = INT( CDUMMY(1) )
1.6 bertrand 875: CALL ZGESVJ( 'L', 'U', 'N', N, N, U, LDU, SVA,
876: $ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
1.9 ! bertrand 877: LWRK_ZGESVJ = INT( CDUMMY(1) )
1.6 bertrand 878: CALL ZGESVJ( 'U', 'U', 'N', N, N, U, LDU, SVA,
879: $ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
1.9 ! bertrand 880: LWRK_ZGESVJU = INT( CDUMMY(1) )
1.6 bertrand 881: CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,
882: $ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
1.9 ! bertrand 883: LWRK_ZGESVJV = INT( CDUMMY(1) )
1.6 bertrand 884: CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,
885: $ V, LDV, CDUMMY, -1, IERR )
1.9 ! bertrand 886: LWRK_ZUNMLQ = INT( CDUMMY(1) )
1.6 bertrand 887: IF ( ERREST ) THEN
888: OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON,
889: $ 2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF,
890: $ 2*N+LWRK_ZGEQP3N,
891: $ 2*N+N**2+N+LWRK_ZGELQF,
892: $ 2*N+N**2+N+N**2+LWCON,
893: $ 2*N+N**2+N+LWRK_ZGESVJ,
894: $ 2*N+N**2+N+LWRK_ZGESVJV,
895: $ 2*N+N**2+N+LWRK_ZUNMQR,
896: $ 2*N+N**2+N+LWRK_ZUNMLQ,
897: $ N+N**2+LWRK_ZGESVJU,
898: $ N+LWRK_ZUNMQRM )
899: ELSE
900: OPTWRK = MAX( N+LWRK_ZGEQP3,
901: $ 2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF,
902: $ 2*N+LWRK_ZGEQP3N,
903: $ 2*N+N**2+N+LWRK_ZGELQF,
904: $ 2*N+N**2+N+N**2+LWCON,
905: $ 2*N+N**2+N+LWRK_ZGESVJ,
906: $ 2*N+N**2+N+LWRK_ZGESVJV,
907: $ 2*N+N**2+N+LWRK_ZUNMQR,
908: $ 2*N+N**2+N+LWRK_ZUNMLQ,
909: $ N+N**2+LWRK_ZGESVJU,
910: $ N+LWRK_ZUNMQRM )
911: END IF
912: ELSE
913: CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,
914: $ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
1.9 ! bertrand 915: LWRK_ZGESVJV = INT( CDUMMY(1) )
1.6 bertrand 916: CALL ZUNMQR( 'L', 'N', N, N, N, CDUMMY, N, CDUMMY,
917: $ V, LDV, CDUMMY, -1, IERR )
1.9 ! bertrand 918: LWRK_ZUNMQR = INT( CDUMMY(1) )
1.6 bertrand 919: CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
920: $ LDU, CDUMMY, -1, IERR )
1.9 ! bertrand 921: LWRK_ZUNMQRM = INT( CDUMMY(1) )
1.6 bertrand 922: IF ( ERREST ) THEN
923: OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON,
924: $ 2*N+LWRK_ZGEQRF, 2*N+N**2,
925: $ 2*N+N**2+LWRK_ZGESVJV,
926: $ 2*N+N**2+N+LWRK_ZUNMQR,N+LWRK_ZUNMQRM )
927: ELSE
928: OPTWRK = MAX( N+LWRK_ZGEQP3, 2*N+LWRK_ZGEQRF,
929: $ 2*N+N**2, 2*N+N**2+LWRK_ZGESVJV,
930: $ 2*N+N**2+N+LWRK_ZUNMQR,
931: $ N+LWRK_ZUNMQRM )
932: END IF
933: END IF
934: END IF
935: IF ( L2TRAN .OR. ROWPIV ) THEN
936: MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ, LRWCON )
937: ELSE
938: MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
939: END IF
940: END IF
941: MINWRK = MAX( 2, MINWRK )
1.8 bertrand 942: OPTWRK = MAX( MINWRK, OPTWRK )
1.6 bertrand 943: IF ( LWORK .LT. MINWRK .AND. (.NOT.LQUERY) ) INFO = - 17
944: IF ( LRWORK .LT. MINRWRK .AND. (.NOT.LQUERY) ) INFO = - 19
945: END IF
946: *
947: IF ( INFO .NE. 0 ) THEN
948: * #:(
949: CALL XERBLA( 'ZGEJSV', - INFO )
950: RETURN
951: ELSE IF ( LQUERY ) THEN
952: CWORK(1) = OPTWRK
953: CWORK(2) = MINWRK
954: RWORK(1) = MINRWRK
955: IWORK(1) = MAX( 4, MINIWRK )
956: RETURN
957: END IF
958: *
959: * Quick return for void matrix (Y3K safe)
960: * #:)
961: IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
962: IWORK(1:4) = 0
963: RWORK(1:7) = 0
964: RETURN
965: ENDIF
966: *
967: * Determine whether the matrix U should be M x N or M x M
968: *
969: IF ( LSVEC ) THEN
970: N1 = N
971: IF ( LSAME( JOBU, 'F' ) ) N1 = M
972: END IF
973: *
974: * Set numerical parameters
975: *
976: *! NOTE: Make sure DLAMCH() does not fail on the target architecture.
977: *
978: EPSLN = DLAMCH('Epsilon')
979: SFMIN = DLAMCH('SafeMinimum')
980: SMALL = SFMIN / EPSLN
981: BIG = DLAMCH('O')
982: * BIG = ONE / SFMIN
983: *
984: * Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
985: *
986: *(!) If necessary, scale SVA() to protect the largest norm from
987: * overflow. It is possible that this scaling pushes the smallest
988: * column norm left from the underflow threshold (extreme case).
989: *
990: SCALEM = ONE / SQRT(DBLE(M)*DBLE(N))
991: NOSCAL = .TRUE.
992: GOSCAL = .TRUE.
993: DO 1874 p = 1, N
994: AAPP = ZERO
995: AAQQ = ONE
996: CALL ZLASSQ( M, A(1,p), 1, AAPP, AAQQ )
997: IF ( AAPP .GT. BIG ) THEN
998: INFO = - 9
999: CALL XERBLA( 'ZGEJSV', -INFO )
1000: RETURN
1001: END IF
1002: AAQQ = SQRT(AAQQ)
1003: IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL ) THEN
1004: SVA(p) = AAPP * AAQQ
1005: ELSE
1006: NOSCAL = .FALSE.
1007: SVA(p) = AAPP * ( AAQQ * SCALEM )
1008: IF ( GOSCAL ) THEN
1009: GOSCAL = .FALSE.
1010: CALL DSCAL( p-1, SCALEM, SVA, 1 )
1011: END IF
1012: END IF
1013: 1874 CONTINUE
1014: *
1015: IF ( NOSCAL ) SCALEM = ONE
1016: *
1017: AAPP = ZERO
1018: AAQQ = BIG
1019: DO 4781 p = 1, N
1020: AAPP = MAX( AAPP, SVA(p) )
1021: IF ( SVA(p) .NE. ZERO ) AAQQ = MIN( AAQQ, SVA(p) )
1022: 4781 CONTINUE
1023: *
1024: * Quick return for zero M x N matrix
1025: * #:)
1026: IF ( AAPP .EQ. ZERO ) THEN
1027: IF ( LSVEC ) CALL ZLASET( 'G', M, N1, CZERO, CONE, U, LDU )
1028: IF ( RSVEC ) CALL ZLASET( 'G', N, N, CZERO, CONE, V, LDV )
1029: RWORK(1) = ONE
1030: RWORK(2) = ONE
1031: IF ( ERREST ) RWORK(3) = ONE
1032: IF ( LSVEC .AND. RSVEC ) THEN
1033: RWORK(4) = ONE
1034: RWORK(5) = ONE
1035: END IF
1036: IF ( L2TRAN ) THEN
1037: RWORK(6) = ZERO
1038: RWORK(7) = ZERO
1039: END IF
1040: IWORK(1) = 0
1041: IWORK(2) = 0
1042: IWORK(3) = 0
1043: IWORK(4) = -1
1044: RETURN
1045: END IF
1046: *
1047: * Issue warning if denormalized column norms detected. Override the
1048: * high relative accuracy request. Issue licence to kill nonzero columns
1049: * (set them to zero) whose norm is less than sigma_max / BIG (roughly).
1050: * #:(
1051: WARNING = 0
1052: IF ( AAQQ .LE. SFMIN ) THEN
1053: L2RANK = .TRUE.
1054: L2KILL = .TRUE.
1055: WARNING = 1
1056: END IF
1057: *
1058: * Quick return for one-column matrix
1059: * #:)
1060: IF ( N .EQ. 1 ) THEN
1061: *
1062: IF ( LSVEC ) THEN
1063: CALL ZLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR )
1064: CALL ZLACPY( 'A', M, 1, A, LDA, U, LDU )
1065: * computing all M left singular vectors of the M x 1 matrix
1066: IF ( N1 .NE. N ) THEN
1067: CALL ZGEQRF( M, N, U,LDU, CWORK, CWORK(N+1),LWORK-N,IERR )
1068: CALL ZUNGQR( M,N1,1, U,LDU,CWORK,CWORK(N+1),LWORK-N,IERR )
1069: CALL ZCOPY( M, A(1,1), 1, U(1,1), 1 )
1070: END IF
1071: END IF
1072: IF ( RSVEC ) THEN
1073: V(1,1) = CONE
1074: END IF
1075: IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN
1076: SVA(1) = SVA(1) / SCALEM
1077: SCALEM = ONE
1078: END IF
1079: RWORK(1) = ONE / SCALEM
1080: RWORK(2) = ONE
1081: IF ( SVA(1) .NE. ZERO ) THEN
1082: IWORK(1) = 1
1083: IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN
1084: IWORK(2) = 1
1085: ELSE
1086: IWORK(2) = 0
1087: END IF
1088: ELSE
1089: IWORK(1) = 0
1090: IWORK(2) = 0
1091: END IF
1092: IWORK(3) = 0
1093: IWORK(4) = -1
1094: IF ( ERREST ) RWORK(3) = ONE
1095: IF ( LSVEC .AND. RSVEC ) THEN
1096: RWORK(4) = ONE
1097: RWORK(5) = ONE
1098: END IF
1099: IF ( L2TRAN ) THEN
1100: RWORK(6) = ZERO
1101: RWORK(7) = ZERO
1102: END IF
1103: RETURN
1104: *
1105: END IF
1106: *
1107: TRANSP = .FALSE.
1108: *
1109: AATMAX = -ONE
1110: AATMIN = BIG
1111: IF ( ROWPIV .OR. L2TRAN ) THEN
1112: *
1113: * Compute the row norms, needed to determine row pivoting sequence
1114: * (in the case of heavily row weighted A, row pivoting is strongly
1115: * advised) and to collect information needed to compare the
1116: * structures of A * A^* and A^* * A (in the case L2TRAN.EQ..TRUE.).
1117: *
1118: IF ( L2TRAN ) THEN
1119: DO 1950 p = 1, M
1120: XSC = ZERO
1121: TEMP1 = ONE
1122: CALL ZLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
1123: * ZLASSQ gets both the ell_2 and the ell_infinity norm
1124: * in one pass through the vector
1125: RWORK(M+p) = XSC * SCALEM
1126: RWORK(p) = XSC * (SCALEM*SQRT(TEMP1))
1127: AATMAX = MAX( AATMAX, RWORK(p) )
1128: IF (RWORK(p) .NE. ZERO)
1129: $ AATMIN = MIN(AATMIN,RWORK(p))
1130: 1950 CONTINUE
1131: ELSE
1132: DO 1904 p = 1, M
1133: RWORK(M+p) = SCALEM*ABS( A(p,IZAMAX(N,A(p,1),LDA)) )
1134: AATMAX = MAX( AATMAX, RWORK(M+p) )
1135: AATMIN = MIN( AATMIN, RWORK(M+p) )
1136: 1904 CONTINUE
1137: END IF
1138: *
1139: END IF
1140: *
1141: * For square matrix A try to determine whether A^* would be better
1142: * input for the preconditioned Jacobi SVD, with faster convergence.
1143: * The decision is based on an O(N) function of the vector of column
1144: * and row norms of A, based on the Shannon entropy. This should give
1145: * the right choice in most cases when the difference actually matters.
1146: * It may fail and pick the slower converging side.
1147: *
1148: ENTRA = ZERO
1149: ENTRAT = ZERO
1150: IF ( L2TRAN ) THEN
1151: *
1152: XSC = ZERO
1153: TEMP1 = ONE
1154: CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
1155: TEMP1 = ONE / TEMP1
1156: *
1157: ENTRA = ZERO
1158: DO 1113 p = 1, N
1159: BIG1 = ( ( SVA(p) / XSC )**2 ) * TEMP1
1160: IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1)
1161: 1113 CONTINUE
1162: ENTRA = - ENTRA / DLOG(DBLE(N))
1163: *
1164: * Now, SVA().^2/Trace(A^* * A) is a point in the probability simplex.
1165: * It is derived from the diagonal of A^* * A. Do the same with the
1166: * diagonal of A * A^*, compute the entropy of the corresponding
1167: * probability distribution. Note that A * A^* and A^* * A have the
1168: * same trace.
1169: *
1170: ENTRAT = ZERO
1171: DO 1114 p = 1, M
1172: BIG1 = ( ( RWORK(p) / XSC )**2 ) * TEMP1
1173: IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1)
1174: 1114 CONTINUE
1175: ENTRAT = - ENTRAT / DLOG(DBLE(M))
1176: *
1177: * Analyze the entropies and decide A or A^*. Smaller entropy
1178: * usually means better input for the algorithm.
1179: *
1180: TRANSP = ( ENTRAT .LT. ENTRA )
1181: *
1182: * If A^* is better than A, take the adjoint of A. This is allowed
1183: * only for square matrices, M=N.
1184: IF ( TRANSP ) THEN
1185: * In an optimal implementation, this trivial transpose
1186: * should be replaced with faster transpose.
1187: DO 1115 p = 1, N - 1
1188: A(p,p) = CONJG(A(p,p))
1189: DO 1116 q = p + 1, N
1190: CTEMP = CONJG(A(q,p))
1191: A(q,p) = CONJG(A(p,q))
1192: A(p,q) = CTEMP
1193: 1116 CONTINUE
1194: 1115 CONTINUE
1195: A(N,N) = CONJG(A(N,N))
1196: DO 1117 p = 1, N
1197: RWORK(M+p) = SVA(p)
1198: SVA(p) = RWORK(p)
1199: * previously computed row 2-norms are now column 2-norms
1200: * of the transposed matrix
1201: 1117 CONTINUE
1202: TEMP1 = AAPP
1203: AAPP = AATMAX
1204: AATMAX = TEMP1
1205: TEMP1 = AAQQ
1206: AAQQ = AATMIN
1207: AATMIN = TEMP1
1208: KILL = LSVEC
1209: LSVEC = RSVEC
1210: RSVEC = KILL
1211: IF ( LSVEC ) N1 = N
1212: *
1213: ROWPIV = .TRUE.
1214: END IF
1215: *
1216: END IF
1217: * END IF L2TRAN
1218: *
1219: * Scale the matrix so that its maximal singular value remains less
1220: * than SQRT(BIG) -- the matrix is scaled so that its maximal column
1221: * has Euclidean norm equal to SQRT(BIG/N). The only reason to keep
1222: * SQRT(BIG) instead of BIG is the fact that ZGEJSV uses LAPACK and
1223: * BLAS routines that, in some implementations, are not capable of
1224: * working in the full interval [SFMIN,BIG] and that they may provoke
1225: * overflows in the intermediate results. If the singular values spread
1226: * from SFMIN to BIG, then ZGESVJ will compute them. So, in that case,
1227: * one should use ZGESVJ instead of ZGEJSV.
1228: * >> change in the April 2016 update: allow bigger range, i.e. the
1229: * largest column is allowed up to BIG/N and ZGESVJ will do the rest.
1230: BIG1 = SQRT( BIG )
1231: TEMP1 = SQRT( BIG / DBLE(N) )
1232: * TEMP1 = BIG/DBLE(N)
1233: *
1234: CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
1235: IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN
1236: AAQQ = ( AAQQ / AAPP ) * TEMP1
1237: ELSE
1238: AAQQ = ( AAQQ * TEMP1 ) / AAPP
1239: END IF
1240: TEMP1 = TEMP1 * SCALEM
1241: CALL ZLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR )
1242: *
1243: * To undo scaling at the end of this procedure, multiply the
1244: * computed singular values with USCAL2 / USCAL1.
1245: *
1246: USCAL1 = TEMP1
1247: USCAL2 = AAPP
1248: *
1249: IF ( L2KILL ) THEN
1250: * L2KILL enforces computation of nonzero singular values in
1251: * the restricted range of condition number of the initial A,
1252: * sigma_max(A) / sigma_min(A) approx. SQRT(BIG)/SQRT(SFMIN).
1253: XSC = SQRT( SFMIN )
1254: ELSE
1255: XSC = SMALL
1256: *
1257: * Now, if the condition number of A is too big,
1258: * sigma_max(A) / sigma_min(A) .GT. SQRT(BIG/N) * EPSLN / SFMIN,
1259: * as a precaution measure, the full SVD is computed using ZGESVJ
1260: * with accumulated Jacobi rotations. This provides numerically
1261: * more robust computation, at the cost of slightly increased run
1262: * time. Depending on the concrete implementation of BLAS and LAPACK
1263: * (i.e. how they behave in presence of extreme ill-conditioning) the
1264: * implementor may decide to remove this switch.
1265: IF ( ( AAQQ.LT.SQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN
1266: JRACC = .TRUE.
1267: END IF
1268: *
1269: END IF
1270: IF ( AAQQ .LT. XSC ) THEN
1271: DO 700 p = 1, N
1272: IF ( SVA(p) .LT. XSC ) THEN
1273: CALL ZLASET( 'A', M, 1, CZERO, CZERO, A(1,p), LDA )
1274: SVA(p) = ZERO
1275: END IF
1276: 700 CONTINUE
1277: END IF
1278: *
1279: * Preconditioning using QR factorization with pivoting
1280: *
1281: IF ( ROWPIV ) THEN
1282: * Optional row permutation (Bjoerck row pivoting):
1283: * A result by Cox and Higham shows that the Bjoerck's
1284: * row pivoting combined with standard column pivoting
1285: * has similar effect as Powell-Reid complete pivoting.
1286: * The ell-infinity norms of A are made nonincreasing.
1287: IF ( ( LSVEC .AND. RSVEC ) .AND. .NOT.( JRACC ) ) THEN
1288: IWOFF = 2*N
1289: ELSE
1290: IWOFF = N
1291: END IF
1292: DO 1952 p = 1, M - 1
1293: q = IDAMAX( M-p+1, RWORK(M+p), 1 ) + p - 1
1294: IWORK(IWOFF+p) = q
1295: IF ( p .NE. q ) THEN
1296: TEMP1 = RWORK(M+p)
1297: RWORK(M+p) = RWORK(M+q)
1298: RWORK(M+q) = TEMP1
1299: END IF
1300: 1952 CONTINUE
1301: CALL ZLASWP( N, A, LDA, 1, M-1, IWORK(IWOFF+1), 1 )
1302: END IF
1303: *
1304: * End of the preparation phase (scaling, optional sorting and
1305: * transposing, optional flushing of small columns).
1306: *
1307: * Preconditioning
1308: *
1309: * If the full SVD is needed, the right singular vectors are computed
1310: * from a matrix equation, and for that we need theoretical analysis
1311: * of the Businger-Golub pivoting. So we use ZGEQP3 as the first RR QRF.
1312: * In all other cases the first RR QRF can be chosen by other criteria
1313: * (eg speed by replacing global with restricted window pivoting, such
1314: * as in xGEQPX from TOMS # 782). Good results will be obtained using
1315: * xGEQPX with properly (!) chosen numerical parameters.
1.9 ! bertrand 1316: * Any improvement of ZGEQP3 improves overall performance of ZGEJSV.
1.6 bertrand 1317: *
1318: * A * P1 = Q1 * [ R1^* 0]^*:
1319: DO 1963 p = 1, N
1320: * .. all columns are free columns
1321: IWORK(p) = 0
1322: 1963 CONTINUE
1323: CALL ZGEQP3( M, N, A, LDA, IWORK, CWORK, CWORK(N+1), LWORK-N,
1324: $ RWORK, IERR )
1325: *
1326: * The upper triangular matrix R1 from the first QRF is inspected for
1327: * rank deficiency and possibilities for deflation, or possible
1328: * ill-conditioning. Depending on the user specified flag L2RANK,
1329: * the procedure explores possibilities to reduce the numerical
1330: * rank by inspecting the computed upper triangular factor. If
1331: * L2RANK or L2ABER are up, then ZGEJSV will compute the SVD of
1332: * A + dA, where ||dA|| <= f(M,N)*EPSLN.
1333: *
1334: NR = 1
1335: IF ( L2ABER ) THEN
1336: * Standard absolute error bound suffices. All sigma_i with
1337: * sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
1.8 bertrand 1338: * aggressive enforcement of lower numerical rank by introducing a
1.6 bertrand 1339: * backward error of the order of N*EPSLN*||A||.
1340: TEMP1 = SQRT(DBLE(N))*EPSLN
1341: DO 3001 p = 2, N
1342: IF ( ABS(A(p,p)) .GE. (TEMP1*ABS(A(1,1))) ) THEN
1343: NR = NR + 1
1344: ELSE
1345: GO TO 3002
1346: END IF
1347: 3001 CONTINUE
1348: 3002 CONTINUE
1349: ELSE IF ( L2RANK ) THEN
1.8 bertrand 1350: * .. similarly as above, only slightly more gentle (less aggressive).
1.6 bertrand 1351: * Sudden drop on the diagonal of R1 is used as the criterion for
1352: * close-to-rank-deficient.
1353: TEMP1 = SQRT(SFMIN)
1354: DO 3401 p = 2, N
1355: IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR.
1356: $ ( ABS(A(p,p)) .LT. SMALL ) .OR.
1357: $ ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402
1358: NR = NR + 1
1359: 3401 CONTINUE
1360: 3402 CONTINUE
1361: *
1362: ELSE
1363: * The goal is high relative accuracy. However, if the matrix
1364: * has high scaled condition number the relative accuracy is in
1365: * general not feasible. Later on, a condition number estimator
1366: * will be deployed to estimate the scaled condition number.
1367: * Here we just remove the underflowed part of the triangular
1368: * factor. This prevents the situation in which the code is
1369: * working hard to get the accuracy not warranted by the data.
1370: TEMP1 = SQRT(SFMIN)
1371: DO 3301 p = 2, N
1372: IF ( ( ABS(A(p,p)) .LT. SMALL ) .OR.
1373: $ ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302
1374: NR = NR + 1
1375: 3301 CONTINUE
1376: 3302 CONTINUE
1377: *
1378: END IF
1379: *
1380: ALMORT = .FALSE.
1381: IF ( NR .EQ. N ) THEN
1382: MAXPRJ = ONE
1383: DO 3051 p = 2, N
1384: TEMP1 = ABS(A(p,p)) / SVA(IWORK(p))
1385: MAXPRJ = MIN( MAXPRJ, TEMP1 )
1386: 3051 CONTINUE
1387: IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE.
1388: END IF
1389: *
1390: *
1391: SCONDA = - ONE
1392: CONDR1 = - ONE
1393: CONDR2 = - ONE
1394: *
1395: IF ( ERREST ) THEN
1396: IF ( N .EQ. NR ) THEN
1397: IF ( RSVEC ) THEN
1398: * .. V is available as workspace
1399: CALL ZLACPY( 'U', N, N, A, LDA, V, LDV )
1400: DO 3053 p = 1, N
1401: TEMP1 = SVA(IWORK(p))
1402: CALL ZDSCAL( p, ONE/TEMP1, V(1,p), 1 )
1403: 3053 CONTINUE
1404: IF ( LSVEC )THEN
1405: CALL ZPOCON( 'U', N, V, LDV, ONE, TEMP1,
1406: $ CWORK(N+1), RWORK, IERR )
1407: ELSE
1408: CALL ZPOCON( 'U', N, V, LDV, ONE, TEMP1,
1409: $ CWORK, RWORK, IERR )
1410: END IF
1411: *
1412: ELSE IF ( LSVEC ) THEN
1413: * .. U is available as workspace
1414: CALL ZLACPY( 'U', N, N, A, LDA, U, LDU )
1415: DO 3054 p = 1, N
1416: TEMP1 = SVA(IWORK(p))
1417: CALL ZDSCAL( p, ONE/TEMP1, U(1,p), 1 )
1418: 3054 CONTINUE
1419: CALL ZPOCON( 'U', N, U, LDU, ONE, TEMP1,
1420: $ CWORK(N+1), RWORK, IERR )
1421: ELSE
1422: CALL ZLACPY( 'U', N, N, A, LDA, CWORK, N )
1423: *[] CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N )
1424: * Change: here index shifted by N to the left, CWORK(1:N)
1425: * not needed for SIGMA only computation
1426: DO 3052 p = 1, N
1427: TEMP1 = SVA(IWORK(p))
1428: *[] CALL ZDSCAL( p, ONE/TEMP1, CWORK(N+(p-1)*N+1), 1 )
1429: CALL ZDSCAL( p, ONE/TEMP1, CWORK((p-1)*N+1), 1 )
1430: 3052 CONTINUE
1431: * .. the columns of R are scaled to have unit Euclidean lengths.
1432: *[] CALL ZPOCON( 'U', N, CWORK(N+1), N, ONE, TEMP1,
1433: *[] $ CWORK(N+N*N+1), RWORK, IERR )
1434: CALL ZPOCON( 'U', N, CWORK, N, ONE, TEMP1,
1435: $ CWORK(N*N+1), RWORK, IERR )
1436: *
1437: END IF
1438: IF ( TEMP1 .NE. ZERO ) THEN
1439: SCONDA = ONE / SQRT(TEMP1)
1440: ELSE
1441: SCONDA = - ONE
1442: END IF
1443: * SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1).
1444: * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
1445: ELSE
1446: SCONDA = - ONE
1447: END IF
1448: END IF
1449: *
1450: L2PERT = L2PERT .AND. ( ABS( A(1,1)/A(NR,NR) ) .GT. SQRT(BIG1) )
1451: * If there is no violent scaling, artificial perturbation is not needed.
1452: *
1453: * Phase 3:
1454: *
1455: IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
1456: *
1457: * Singular Values only
1458: *
1459: * .. transpose A(1:NR,1:N)
1460: DO 1946 p = 1, MIN( N-1, NR )
1461: CALL ZCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 )
1462: CALL ZLACGV( N-p+1, A(p,p), 1 )
1463: 1946 CONTINUE
1464: IF ( NR .EQ. N ) A(N,N) = CONJG(A(N,N))
1465: *
1466: * The following two DO-loops introduce small relative perturbation
1467: * into the strict upper triangle of the lower triangular matrix.
1468: * Small entries below the main diagonal are also changed.
1469: * This modification is useful if the computing environment does not
1470: * provide/allow FLUSH TO ZERO underflow, for it prevents many
1471: * annoying denormalized numbers in case of strongly scaled matrices.
1472: * The perturbation is structured so that it does not introduce any
1473: * new perturbation of the singular values, and it does not destroy
1474: * the job done by the preconditioner.
1475: * The licence for this perturbation is in the variable L2PERT, which
1476: * should be .FALSE. if FLUSH TO ZERO underflow is active.
1477: *
1478: IF ( .NOT. ALMORT ) THEN
1479: *
1480: IF ( L2PERT ) THEN
1481: * XSC = SQRT(SMALL)
1482: XSC = EPSLN / DBLE(N)
1483: DO 4947 q = 1, NR
1484: CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
1485: DO 4949 p = 1, N
1486: IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) )
1487: $ .OR. ( p .LT. q ) )
1488: * $ A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) )
1489: $ A(p,q) = CTEMP
1490: 4949 CONTINUE
1491: 4947 CONTINUE
1492: ELSE
1493: CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, A(1,2),LDA )
1494: END IF
1495: *
1496: * .. second preconditioning using the QR factorization
1497: *
1498: CALL ZGEQRF( N,NR, A,LDA, CWORK, CWORK(N+1),LWORK-N, IERR )
1499: *
1500: * .. and transpose upper to lower triangular
1501: DO 1948 p = 1, NR - 1
1502: CALL ZCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 )
1503: CALL ZLACGV( NR-p+1, A(p,p), 1 )
1504: 1948 CONTINUE
1505: *
1506: END IF
1507: *
1508: * Row-cyclic Jacobi SVD algorithm with column pivoting
1509: *
1510: * .. again some perturbation (a "background noise") is added
1511: * to drown denormals
1512: IF ( L2PERT ) THEN
1513: * XSC = SQRT(SMALL)
1514: XSC = EPSLN / DBLE(N)
1515: DO 1947 q = 1, NR
1516: CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
1517: DO 1949 p = 1, NR
1518: IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) )
1519: $ .OR. ( p .LT. q ) )
1520: * $ A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) )
1521: $ A(p,q) = CTEMP
1522: 1949 CONTINUE
1523: 1947 CONTINUE
1524: ELSE
1525: CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, A(1,2), LDA )
1526: END IF
1527: *
1528: * .. and one-sided Jacobi rotations are started on a lower
1529: * triangular matrix (plus perturbation which is ignored in
1530: * the part which destroys triangular form (confusing?!))
1531: *
1532: CALL ZGESVJ( 'L', 'N', 'N', NR, NR, A, LDA, SVA,
1533: $ N, V, LDV, CWORK, LWORK, RWORK, LRWORK, INFO )
1534: *
1535: SCALEM = RWORK(1)
1536: NUMRANK = NINT(RWORK(2))
1537: *
1538: *
1539: ELSE IF ( ( RSVEC .AND. ( .NOT. LSVEC ) .AND. ( .NOT. JRACC ) )
1540: $ .OR.
1541: $ ( JRACC .AND. ( .NOT. LSVEC ) .AND. ( NR .NE. N ) ) ) THEN
1542: *
1543: * -> Singular Values and Right Singular Vectors <-
1544: *
1545: IF ( ALMORT ) THEN
1546: *
1547: * .. in this case NR equals N
1548: DO 1998 p = 1, NR
1549: CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1550: CALL ZLACGV( N-p+1, V(p,p), 1 )
1551: 1998 CONTINUE
1552: CALL ZLASET( 'U', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV )
1553: *
1554: CALL ZGESVJ( 'L','U','N', N, NR, V, LDV, SVA, NR, A, LDA,
1555: $ CWORK, LWORK, RWORK, LRWORK, INFO )
1556: SCALEM = RWORK(1)
1557: NUMRANK = NINT(RWORK(2))
1558:
1559: ELSE
1560: *
1561: * .. two more QR factorizations ( one QRF is not enough, two require
1562: * accumulated product of Jacobi rotations, three are perfect )
1563: *
1564: CALL ZLASET( 'L', NR-1,NR-1, CZERO, CZERO, A(2,1), LDA )
1565: CALL ZGELQF( NR,N, A, LDA, CWORK, CWORK(N+1), LWORK-N, IERR)
1566: CALL ZLACPY( 'L', NR, NR, A, LDA, V, LDV )
1567: CALL ZLASET( 'U', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV )
1568: CALL ZGEQRF( NR, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
1569: $ LWORK-2*N, IERR )
1570: DO 8998 p = 1, NR
1571: CALL ZCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
1572: CALL ZLACGV( NR-p+1, V(p,p), 1 )
1573: 8998 CONTINUE
1574: CALL ZLASET('U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV)
1575: *
1576: CALL ZGESVJ( 'L', 'U','N', NR, NR, V,LDV, SVA, NR, U,
1577: $ LDU, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )
1578: SCALEM = RWORK(1)
1579: NUMRANK = NINT(RWORK(2))
1580: IF ( NR .LT. N ) THEN
1581: CALL ZLASET( 'A',N-NR, NR, CZERO,CZERO, V(NR+1,1), LDV )
1582: CALL ZLASET( 'A',NR, N-NR, CZERO,CZERO, V(1,NR+1), LDV )
1583: CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE, V(NR+1,NR+1),LDV )
1584: END IF
1585: *
1586: CALL ZUNMLQ( 'L', 'C', N, N, NR, A, LDA, CWORK,
1587: $ V, LDV, CWORK(N+1), LWORK-N, IERR )
1588: *
1589: END IF
1590: * .. permute the rows of V
1591: * DO 8991 p = 1, N
1592: * CALL ZCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA )
1593: * 8991 CONTINUE
1594: * CALL ZLACPY( 'All', N, N, A, LDA, V, LDV )
1595: CALL ZLAPMR( .FALSE., N, N, V, LDV, IWORK )
1596: *
1597: IF ( TRANSP ) THEN
1598: CALL ZLACPY( 'A', N, N, V, LDV, U, LDU )
1599: END IF
1600: *
1601: ELSE IF ( JRACC .AND. (.NOT. LSVEC) .AND. ( NR.EQ. N ) ) THEN
1602: *
1603: CALL ZLASET( 'L', N-1,N-1, CZERO, CZERO, A(2,1), LDA )
1604: *
1605: CALL ZGESVJ( 'U','N','V', N, N, A, LDA, SVA, N, V, LDV,
1606: $ CWORK, LWORK, RWORK, LRWORK, INFO )
1607: SCALEM = RWORK(1)
1608: NUMRANK = NINT(RWORK(2))
1609: CALL ZLAPMR( .FALSE., N, N, V, LDV, IWORK )
1610: *
1611: ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN
1612: *
1613: * .. Singular Values and Left Singular Vectors ..
1614: *
1615: * .. second preconditioning step to avoid need to accumulate
1616: * Jacobi rotations in the Jacobi iterations.
1617: DO 1965 p = 1, NR
1618: CALL ZCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 )
1619: CALL ZLACGV( N-p+1, U(p,p), 1 )
1620: 1965 CONTINUE
1621: CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
1622: *
1623: CALL ZGEQRF( N, NR, U, LDU, CWORK(N+1), CWORK(2*N+1),
1624: $ LWORK-2*N, IERR )
1625: *
1626: DO 1967 p = 1, NR - 1
1627: CALL ZCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 )
1628: CALL ZLACGV( N-p+1, U(p,p), 1 )
1629: 1967 CONTINUE
1630: CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
1631: *
1632: CALL ZGESVJ( 'L', 'U', 'N', NR,NR, U, LDU, SVA, NR, A,
1633: $ LDA, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )
1634: SCALEM = RWORK(1)
1635: NUMRANK = NINT(RWORK(2))
1636: *
1637: IF ( NR .LT. M ) THEN
1638: CALL ZLASET( 'A', M-NR, NR,CZERO, CZERO, U(NR+1,1), LDU )
1639: IF ( NR .LT. N1 ) THEN
1640: CALL ZLASET( 'A',NR, N1-NR, CZERO, CZERO, U(1,NR+1),LDU )
1641: CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,U(NR+1,NR+1),LDU )
1642: END IF
1643: END IF
1644: *
1645: CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
1646: $ LDU, CWORK(N+1), LWORK-N, IERR )
1647: *
1648: IF ( ROWPIV )
1649: $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
1650: *
1651: DO 1974 p = 1, N1
1652: XSC = ONE / DZNRM2( M, U(1,p), 1 )
1653: CALL ZDSCAL( M, XSC, U(1,p), 1 )
1654: 1974 CONTINUE
1655: *
1656: IF ( TRANSP ) THEN
1657: CALL ZLACPY( 'A', N, N, U, LDU, V, LDV )
1658: END IF
1659: *
1660: ELSE
1661: *
1662: * .. Full SVD ..
1663: *
1664: IF ( .NOT. JRACC ) THEN
1665: *
1666: IF ( .NOT. ALMORT ) THEN
1667: *
1668: * Second Preconditioning Step (QRF [with pivoting])
1669: * Note that the composition of TRANSPOSE, QRF and TRANSPOSE is
1670: * equivalent to an LQF CALL. Since in many libraries the QRF
1671: * seems to be better optimized than the LQF, we do explicit
1672: * transpose and use the QRF. This is subject to changes in an
1673: * optimized implementation of ZGEJSV.
1674: *
1675: DO 1968 p = 1, NR
1676: CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1677: CALL ZLACGV( N-p+1, V(p,p), 1 )
1678: 1968 CONTINUE
1679: *
1680: * .. the following two loops perturb small entries to avoid
1681: * denormals in the second QR factorization, where they are
1682: * as good as zeros. This is done to avoid painfully slow
1683: * computation with denormals. The relative size of the perturbation
1684: * is a parameter that can be changed by the implementer.
1685: * This perturbation device will be obsolete on machines with
1686: * properly implemented arithmetic.
1687: * To switch it off, set L2PERT=.FALSE. To remove it from the
1688: * code, remove the action under L2PERT=.TRUE., leave the ELSE part.
1689: * The following two loops should be blocked and fused with the
1690: * transposed copy above.
1691: *
1692: IF ( L2PERT ) THEN
1693: XSC = SQRT(SMALL)
1694: DO 2969 q = 1, NR
1695: CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO)
1696: DO 2968 p = 1, N
1697: IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 )
1698: $ .OR. ( p .LT. q ) )
1699: * $ V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) )
1700: $ V(p,q) = CTEMP
1701: IF ( p .LT. q ) V(p,q) = - V(p,q)
1702: 2968 CONTINUE
1703: 2969 CONTINUE
1704: ELSE
1705: CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV )
1706: END IF
1707: *
1708: * Estimate the row scaled condition number of R1
1709: * (If R1 is rectangular, N > NR, then the condition number
1710: * of the leading NR x NR submatrix is estimated.)
1711: *
1712: CALL ZLACPY( 'L', NR, NR, V, LDV, CWORK(2*N+1), NR )
1713: DO 3950 p = 1, NR
1714: TEMP1 = DZNRM2(NR-p+1,CWORK(2*N+(p-1)*NR+p),1)
1715: CALL ZDSCAL(NR-p+1,ONE/TEMP1,CWORK(2*N+(p-1)*NR+p),1)
1716: 3950 CONTINUE
1717: CALL ZPOCON('L',NR,CWORK(2*N+1),NR,ONE,TEMP1,
1718: $ CWORK(2*N+NR*NR+1),RWORK,IERR)
1719: CONDR1 = ONE / SQRT(TEMP1)
1.8 bertrand 1720: * .. here need a second opinion on the condition number
1.6 bertrand 1721: * .. then assume worst case scenario
1722: * R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
1723: * more conservative <=> CONDR1 .LT. SQRT(DBLE(N))
1724: *
1725: COND_OK = SQRT(SQRT(DBLE(NR)))
1726: *[TP] COND_OK is a tuning parameter.
1727: *
1728: IF ( CONDR1 .LT. COND_OK ) THEN
1729: * .. the second QRF without pivoting. Note: in an optimized
1730: * implementation, this QRF should be implemented as the QRF
1731: * of a lower triangular matrix.
1732: * R1^* = Q2 * R2
1733: CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
1734: $ LWORK-2*N, IERR )
1735: *
1736: IF ( L2PERT ) THEN
1737: XSC = SQRT(SMALL)/EPSLN
1738: DO 3959 p = 2, NR
1739: DO 3958 q = 1, p - 1
1740: CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
1741: $ ZERO)
1742: IF ( ABS(V(q,p)) .LE. TEMP1 )
1743: * $ V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) )
1744: $ V(q,p) = CTEMP
1745: 3958 CONTINUE
1746: 3959 CONTINUE
1747: END IF
1748: *
1749: IF ( NR .NE. N )
1750: $ CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N )
1751: * .. save ...
1752: *
1753: * .. this transposed copy should be better than naive
1754: DO 1969 p = 1, NR - 1
1755: CALL ZCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 )
1756: CALL ZLACGV(NR-p+1, V(p,p), 1 )
1757: 1969 CONTINUE
1758: V(NR,NR)=CONJG(V(NR,NR))
1759: *
1760: CONDR2 = CONDR1
1761: *
1762: ELSE
1763: *
1764: * .. ill-conditioned case: second QRF with pivoting
1.8 bertrand 1765: * Note that windowed pivoting would be equally good
1.6 bertrand 1766: * numerically, and more run-time efficient. So, in
1767: * an optimal implementation, the next call to ZGEQP3
1768: * should be replaced with eg. CALL ZGEQPX (ACM TOMS #782)
1769: * with properly (carefully) chosen parameters.
1770: *
1771: * R1^* * P2 = Q2 * R2
1772: DO 3003 p = 1, NR
1773: IWORK(N+p) = 0
1774: 3003 CONTINUE
1775: CALL ZGEQP3( N, NR, V, LDV, IWORK(N+1), CWORK(N+1),
1776: $ CWORK(2*N+1), LWORK-2*N, RWORK, IERR )
1777: ** CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
1778: ** $ LWORK-2*N, IERR )
1779: IF ( L2PERT ) THEN
1780: XSC = SQRT(SMALL)
1781: DO 3969 p = 2, NR
1782: DO 3968 q = 1, p - 1
1783: CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
1784: $ ZERO)
1785: IF ( ABS(V(q,p)) .LE. TEMP1 )
1786: * $ V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) )
1787: $ V(q,p) = CTEMP
1788: 3968 CONTINUE
1789: 3969 CONTINUE
1790: END IF
1791: *
1792: CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N )
1793: *
1794: IF ( L2PERT ) THEN
1795: XSC = SQRT(SMALL)
1796: DO 8970 p = 2, NR
1797: DO 8971 q = 1, p - 1
1798: CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
1799: $ ZERO)
1800: * V(p,q) = - TEMP1*( V(q,p) / ABS(V(q,p)) )
1801: V(p,q) = - CTEMP
1802: 8971 CONTINUE
1803: 8970 CONTINUE
1804: ELSE
1805: CALL ZLASET( 'L',NR-1,NR-1,CZERO,CZERO,V(2,1),LDV )
1806: END IF
1807: * Now, compute R2 = L3 * Q3, the LQ factorization.
1808: CALL ZGELQF( NR, NR, V, LDV, CWORK(2*N+N*NR+1),
1809: $ CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR )
1810: * .. and estimate the condition number
1811: CALL ZLACPY( 'L',NR,NR,V,LDV,CWORK(2*N+N*NR+NR+1),NR )
1812: DO 4950 p = 1, NR
1813: TEMP1 = DZNRM2( p, CWORK(2*N+N*NR+NR+p), NR )
1814: CALL ZDSCAL( p, ONE/TEMP1, CWORK(2*N+N*NR+NR+p), NR )
1815: 4950 CONTINUE
1816: CALL ZPOCON( 'L',NR,CWORK(2*N+N*NR+NR+1),NR,ONE,TEMP1,
1817: $ CWORK(2*N+N*NR+NR+NR*NR+1),RWORK,IERR )
1818: CONDR2 = ONE / SQRT(TEMP1)
1819: *
1820: *
1821: IF ( CONDR2 .GE. COND_OK ) THEN
1822: * .. save the Householder vectors used for Q3
1.8 bertrand 1823: * (this overwrites the copy of R2, as it will not be
1.6 bertrand 1824: * needed in this branch, but it does not overwritte the
1825: * Huseholder vectors of Q2.).
1826: CALL ZLACPY( 'U', NR, NR, V, LDV, CWORK(2*N+1), N )
1827: * .. and the rest of the information on Q3 is in
1828: * WORK(2*N+N*NR+1:2*N+N*NR+N)
1829: END IF
1830: *
1831: END IF
1832: *
1833: IF ( L2PERT ) THEN
1834: XSC = SQRT(SMALL)
1835: DO 4968 q = 2, NR
1836: CTEMP = XSC * V(q,q)
1837: DO 4969 p = 1, q - 1
1838: * V(p,q) = - TEMP1*( V(p,q) / ABS(V(p,q)) )
1839: V(p,q) = - CTEMP
1840: 4969 CONTINUE
1841: 4968 CONTINUE
1842: ELSE
1843: CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV )
1844: END IF
1845: *
1846: * Second preconditioning finished; continue with Jacobi SVD
1847: * The input matrix is lower trinagular.
1848: *
1849: * Recover the right singular vectors as solution of a well
1850: * conditioned triangular matrix equation.
1851: *
1852: IF ( CONDR1 .LT. COND_OK ) THEN
1853: *
1854: CALL ZGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U, LDU,
1855: $ CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,RWORK,
1856: $ LRWORK, INFO )
1857: SCALEM = RWORK(1)
1858: NUMRANK = NINT(RWORK(2))
1859: DO 3970 p = 1, NR
1860: CALL ZCOPY( NR, V(1,p), 1, U(1,p), 1 )
1861: CALL ZDSCAL( NR, SVA(p), V(1,p), 1 )
1862: 3970 CONTINUE
1863:
1864: * .. pick the right matrix equation and solve it
1865: *
1866: IF ( NR .EQ. N ) THEN
1867: * :)) .. best case, R1 is inverted. The solution of this matrix
1868: * equation is Q2*V2 = the product of the Jacobi rotations
1869: * used in ZGESVJ, premultiplied with the orthogonal matrix
1870: * from the second QR factorization.
1871: CALL ZTRSM('L','U','N','N', NR,NR,CONE, A,LDA, V,LDV)
1872: ELSE
1873: * .. R1 is well conditioned, but non-square. Adjoint of R2
1874: * is inverted to get the product of the Jacobi rotations
1875: * used in ZGESVJ. The Q-factor from the second QR
1876: * factorization is then built in explicitly.
1877: CALL ZTRSM('L','U','C','N',NR,NR,CONE,CWORK(2*N+1),
1878: $ N,V,LDV)
1879: IF ( NR .LT. N ) THEN
1880: CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV)
1881: CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
1882: CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
1883: END IF
1884: CALL ZUNMQR('L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
1885: $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
1886: END IF
1887: *
1888: ELSE IF ( CONDR2 .LT. COND_OK ) THEN
1889: *
1890: * The matrix R2 is inverted. The solution of the matrix equation
1891: * is Q3^* * V3 = the product of the Jacobi rotations (appplied to
1892: * the lower triangular L3 from the LQ factorization of
1893: * R2=L3*Q3), pre-multiplied with the transposed Q3.
1894: CALL ZGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
1895: $ LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR,
1896: $ RWORK, LRWORK, INFO )
1897: SCALEM = RWORK(1)
1898: NUMRANK = NINT(RWORK(2))
1899: DO 3870 p = 1, NR
1900: CALL ZCOPY( NR, V(1,p), 1, U(1,p), 1 )
1901: CALL ZDSCAL( NR, SVA(p), U(1,p), 1 )
1902: 3870 CONTINUE
1903: CALL ZTRSM('L','U','N','N',NR,NR,CONE,CWORK(2*N+1),N,
1904: $ U,LDU)
1905: * .. apply the permutation from the second QR factorization
1906: DO 873 q = 1, NR
1907: DO 872 p = 1, NR
1908: CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
1909: 872 CONTINUE
1910: DO 874 p = 1, NR
1911: U(p,q) = CWORK(2*N+N*NR+NR+p)
1912: 874 CONTINUE
1913: 873 CONTINUE
1914: IF ( NR .LT. N ) THEN
1915: CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
1916: CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
1917: CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
1918: END IF
1919: CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
1920: $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
1921: ELSE
1922: * Last line of defense.
1923: * #:( This is a rather pathological case: no scaled condition
1924: * improvement after two pivoted QR factorizations. Other
1925: * possibility is that the rank revealing QR factorization
1926: * or the condition estimator has failed, or the COND_OK
1927: * is set very close to ONE (which is unnecessary). Normally,
1928: * this branch should never be executed, but in rare cases of
1929: * failure of the RRQR or condition estimator, the last line of
1930: * defense ensures that ZGEJSV completes the task.
1931: * Compute the full SVD of L3 using ZGESVJ with explicit
1932: * accumulation of Jacobi rotations.
1933: CALL ZGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
1934: $ LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR,
1935: $ RWORK, LRWORK, INFO )
1936: SCALEM = RWORK(1)
1937: NUMRANK = NINT(RWORK(2))
1938: IF ( NR .LT. N ) THEN
1939: CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
1940: CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
1941: CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
1942: END IF
1943: CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
1944: $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
1945: *
1946: CALL ZUNMLQ( 'L', 'C', NR, NR, NR, CWORK(2*N+1), N,
1947: $ CWORK(2*N+N*NR+1), U, LDU, CWORK(2*N+N*NR+NR+1),
1948: $ LWORK-2*N-N*NR-NR, IERR )
1949: DO 773 q = 1, NR
1950: DO 772 p = 1, NR
1951: CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
1952: 772 CONTINUE
1953: DO 774 p = 1, NR
1954: U(p,q) = CWORK(2*N+N*NR+NR+p)
1955: 774 CONTINUE
1956: 773 CONTINUE
1957: *
1958: END IF
1959: *
1960: * Permute the rows of V using the (column) permutation from the
1961: * first QRF. Also, scale the columns to make them unit in
1962: * Euclidean norm. This applies to all cases.
1963: *
1964: TEMP1 = SQRT(DBLE(N)) * EPSLN
1965: DO 1972 q = 1, N
1966: DO 972 p = 1, N
1967: CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
1968: 972 CONTINUE
1969: DO 973 p = 1, N
1970: V(p,q) = CWORK(2*N+N*NR+NR+p)
1971: 973 CONTINUE
1972: XSC = ONE / DZNRM2( N, V(1,q), 1 )
1973: IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1974: $ CALL ZDSCAL( N, XSC, V(1,q), 1 )
1975: 1972 CONTINUE
1976: * At this moment, V contains the right singular vectors of A.
1977: * Next, assemble the left singular vector matrix U (M x N).
1978: IF ( NR .LT. M ) THEN
1979: CALL ZLASET('A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU)
1980: IF ( NR .LT. N1 ) THEN
1981: CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
1982: CALL ZLASET('A',M-NR,N1-NR,CZERO,CONE,
1983: $ U(NR+1,NR+1),LDU)
1984: END IF
1985: END IF
1986: *
1987: * The Q matrix from the first QRF is built into the left singular
1988: * matrix U. This applies to all cases.
1989: *
1990: CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
1991: $ LDU, CWORK(N+1), LWORK-N, IERR )
1992:
1993: * The columns of U are normalized. The cost is O(M*N) flops.
1994: TEMP1 = SQRT(DBLE(M)) * EPSLN
1995: DO 1973 p = 1, NR
1996: XSC = ONE / DZNRM2( M, U(1,p), 1 )
1997: IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1998: $ CALL ZDSCAL( M, XSC, U(1,p), 1 )
1999: 1973 CONTINUE
2000: *
2001: * If the initial QRF is computed with row pivoting, the left
2002: * singular vectors must be adjusted.
2003: *
2004: IF ( ROWPIV )
2005: $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
2006: *
2007: ELSE
2008: *
2009: * .. the initial matrix A has almost orthogonal columns and
2010: * the second QRF is not needed
2011: *
2012: CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N )
2013: IF ( L2PERT ) THEN
2014: XSC = SQRT(SMALL)
2015: DO 5970 p = 2, N
2016: CTEMP = XSC * CWORK( N + (p-1)*N + p )
2017: DO 5971 q = 1, p - 1
2018: * CWORK(N+(q-1)*N+p)=-TEMP1 * ( CWORK(N+(p-1)*N+q) /
2019: * $ ABS(CWORK(N+(p-1)*N+q)) )
2020: CWORK(N+(q-1)*N+p)=-CTEMP
2021: 5971 CONTINUE
2022: 5970 CONTINUE
2023: ELSE
2024: CALL ZLASET( 'L',N-1,N-1,CZERO,CZERO,CWORK(N+2),N )
2025: END IF
2026: *
2027: CALL ZGESVJ( 'U', 'U', 'N', N, N, CWORK(N+1), N, SVA,
2028: $ N, U, LDU, CWORK(N+N*N+1), LWORK-N-N*N, RWORK, LRWORK,
2029: $ INFO )
2030: *
2031: SCALEM = RWORK(1)
2032: NUMRANK = NINT(RWORK(2))
2033: DO 6970 p = 1, N
2034: CALL ZCOPY( N, CWORK(N+(p-1)*N+1), 1, U(1,p), 1 )
2035: CALL ZDSCAL( N, SVA(p), CWORK(N+(p-1)*N+1), 1 )
2036: 6970 CONTINUE
2037: *
2038: CALL ZTRSM( 'L', 'U', 'N', 'N', N, N,
2039: $ CONE, A, LDA, CWORK(N+1), N )
2040: DO 6972 p = 1, N
2041: CALL ZCOPY( N, CWORK(N+p), N, V(IWORK(p),1), LDV )
2042: 6972 CONTINUE
2043: TEMP1 = SQRT(DBLE(N))*EPSLN
2044: DO 6971 p = 1, N
2045: XSC = ONE / DZNRM2( N, V(1,p), 1 )
2046: IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
2047: $ CALL ZDSCAL( N, XSC, V(1,p), 1 )
2048: 6971 CONTINUE
2049: *
2050: * Assemble the left singular vector matrix U (M x N).
2051: *
2052: IF ( N .LT. M ) THEN
2053: CALL ZLASET( 'A', M-N, N, CZERO, CZERO, U(N+1,1), LDU )
2054: IF ( N .LT. N1 ) THEN
2055: CALL ZLASET('A',N, N1-N, CZERO, CZERO, U(1,N+1),LDU)
2056: CALL ZLASET( 'A',M-N,N1-N, CZERO, CONE,U(N+1,N+1),LDU)
2057: END IF
2058: END IF
2059: CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
2060: $ LDU, CWORK(N+1), LWORK-N, IERR )
2061: TEMP1 = SQRT(DBLE(M))*EPSLN
2062: DO 6973 p = 1, N1
2063: XSC = ONE / DZNRM2( M, U(1,p), 1 )
2064: IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
2065: $ CALL ZDSCAL( M, XSC, U(1,p), 1 )
2066: 6973 CONTINUE
2067: *
2068: IF ( ROWPIV )
2069: $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
2070: *
2071: END IF
2072: *
2073: * end of the >> almost orthogonal case << in the full SVD
2074: *
2075: ELSE
2076: *
2077: * This branch deploys a preconditioned Jacobi SVD with explicitly
2078: * accumulated rotations. It is included as optional, mainly for
1.8 bertrand 2079: * experimental purposes. It does perform well, and can also be used.
1.6 bertrand 2080: * In this implementation, this branch will be automatically activated
2081: * if the condition number sigma_max(A) / sigma_min(A) is predicted
2082: * to be greater than the overflow threshold. This is because the
2083: * a posteriori computation of the singular vectors assumes robust
2084: * implementation of BLAS and some LAPACK procedures, capable of working
2085: * in presence of extreme values, e.g. when the singular values spread from
2086: * the underflow to the overflow threshold.
2087: *
2088: DO 7968 p = 1, NR
2089: CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
2090: CALL ZLACGV( N-p+1, V(p,p), 1 )
2091: 7968 CONTINUE
2092: *
2093: IF ( L2PERT ) THEN
2094: XSC = SQRT(SMALL/EPSLN)
2095: DO 5969 q = 1, NR
2096: CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO)
2097: DO 5968 p = 1, N
2098: IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 )
2099: $ .OR. ( p .LT. q ) )
2100: * $ V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) )
2101: $ V(p,q) = CTEMP
2102: IF ( p .LT. q ) V(p,q) = - V(p,q)
2103: 5968 CONTINUE
2104: 5969 CONTINUE
2105: ELSE
2106: CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV )
2107: END IF
2108:
2109: CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
2110: $ LWORK-2*N, IERR )
2111: CALL ZLACPY( 'L', N, NR, V, LDV, CWORK(2*N+1), N )
2112: *
2113: DO 7969 p = 1, NR
2114: CALL ZCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 )
2115: CALL ZLACGV( NR-p+1, U(p,p), 1 )
2116: 7969 CONTINUE
2117:
2118: IF ( L2PERT ) THEN
2119: XSC = SQRT(SMALL/EPSLN)
2120: DO 9970 q = 2, NR
2121: DO 9971 p = 1, q - 1
2122: CTEMP = DCMPLX(XSC * MIN(ABS(U(p,p)),ABS(U(q,q))),
2123: $ ZERO)
2124: * U(p,q) = - TEMP1 * ( U(q,p) / ABS(U(q,p)) )
2125: U(p,q) = - CTEMP
2126: 9971 CONTINUE
2127: 9970 CONTINUE
2128: ELSE
2129: CALL ZLASET('U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
2130: END IF
2131:
2132: CALL ZGESVJ( 'L', 'U', 'V', NR, NR, U, LDU, SVA,
2133: $ N, V, LDV, CWORK(2*N+N*NR+1), LWORK-2*N-N*NR,
2134: $ RWORK, LRWORK, INFO )
2135: SCALEM = RWORK(1)
2136: NUMRANK = NINT(RWORK(2))
2137:
2138: IF ( NR .LT. N ) THEN
2139: CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
2140: CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
2141: CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV )
2142: END IF
2143:
2144: CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
2145: $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
2146: *
2147: * Permute the rows of V using the (column) permutation from the
2148: * first QRF. Also, scale the columns to make them unit in
2149: * Euclidean norm. This applies to all cases.
2150: *
2151: TEMP1 = SQRT(DBLE(N)) * EPSLN
2152: DO 7972 q = 1, N
2153: DO 8972 p = 1, N
2154: CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
2155: 8972 CONTINUE
2156: DO 8973 p = 1, N
2157: V(p,q) = CWORK(2*N+N*NR+NR+p)
2158: 8973 CONTINUE
2159: XSC = ONE / DZNRM2( N, V(1,q), 1 )
2160: IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
2161: $ CALL ZDSCAL( N, XSC, V(1,q), 1 )
2162: 7972 CONTINUE
2163: *
2164: * At this moment, V contains the right singular vectors of A.
2165: * Next, assemble the left singular vector matrix U (M x N).
2166: *
2167: IF ( NR .LT. M ) THEN
2168: CALL ZLASET( 'A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU )
2169: IF ( NR .LT. N1 ) THEN
2170: CALL ZLASET('A',NR, N1-NR, CZERO, CZERO, U(1,NR+1),LDU)
2171: CALL ZLASET('A',M-NR,N1-NR, CZERO, CONE,U(NR+1,NR+1),LDU)
2172: END IF
2173: END IF
2174: *
2175: CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
2176: $ LDU, CWORK(N+1), LWORK-N, IERR )
2177: *
2178: IF ( ROWPIV )
2179: $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
2180: *
2181: *
2182: END IF
2183: IF ( TRANSP ) THEN
2184: * .. swap U and V because the procedure worked on A^*
2185: DO 6974 p = 1, N
2186: CALL ZSWAP( N, U(1,p), 1, V(1,p), 1 )
2187: 6974 CONTINUE
2188: END IF
2189: *
2190: END IF
2191: * end of the full SVD
2192: *
2193: * Undo scaling, if necessary (and possible)
2194: *
2195: IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
2196: CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
2197: USCAL1 = ONE
2198: USCAL2 = ONE
2199: END IF
2200: *
2201: IF ( NR .LT. N ) THEN
2202: DO 3004 p = NR+1, N
2203: SVA(p) = ZERO
2204: 3004 CONTINUE
2205: END IF
2206: *
2207: RWORK(1) = USCAL2 * SCALEM
2208: RWORK(2) = USCAL1
2209: IF ( ERREST ) RWORK(3) = SCONDA
2210: IF ( LSVEC .AND. RSVEC ) THEN
2211: RWORK(4) = CONDR1
2212: RWORK(5) = CONDR2
2213: END IF
2214: IF ( L2TRAN ) THEN
2215: RWORK(6) = ENTRA
2216: RWORK(7) = ENTRAT
2217: END IF
2218: *
2219: IWORK(1) = NR
2220: IWORK(2) = NUMRANK
2221: IWORK(3) = WARNING
2222: IF ( TRANSP ) THEN
2223: IWORK(4) = 1
2224: ELSE
2225: IWORK(4) = -1
2226: END IF
2227:
2228: *
2229: RETURN
2230: * ..
2231: * .. END OF ZGEJSV
2232: * ..
2233: END
2234: *
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgejsv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack