Annotation of rpl/lapack/lapack/dgejsv.f, revision 1.21
1.7 bertrand 1: *> \brief \b DGEJSV
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.16 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.7 bertrand 7: *
8: *> \htmlonly
1.16 bertrand 9: *> Download DGEJSV + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgejsv.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgejsv.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgejsv.f">
1.7 bertrand 15: *> [TXT]</a>
1.16 bertrand 16: *> \endhtmlonly
1.7 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
22: * M, N, A, LDA, SVA, U, LDU, V, LDV,
23: * WORK, LWORK, IWORK, INFO )
1.16 bertrand 24: *
1.7 bertrand 25: * .. Scalar Arguments ..
26: * IMPLICIT NONE
27: * INTEGER INFO, LDA, LDU, LDV, LWORK, M, N
28: * ..
29: * .. Array Arguments ..
30: * DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ),
31: * $ WORK( LWORK )
32: * INTEGER IWORK( * )
33: * CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
34: * ..
1.16 bertrand 35: *
1.7 bertrand 36: *
37: *> \par Purpose:
38: * =============
39: *>
40: *> \verbatim
41: *>
42: *> DGEJSV computes the singular value decomposition (SVD) of a real M-by-N
43: *> matrix [A], where M >= N. The SVD of [A] is written as
44: *>
45: *> [A] = [U] * [SIGMA] * [V]^t,
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) orthonormal matrix, and
49: *> [V] is an N-by-N orthogonal 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.
1.13 bertrand 54: *> DGEJSV can sometimes compute tiny singular values and their singular vectors much
1.14 bertrand 55: *> more accurately than other SVD routines, see below under Further Details.
56: *> \endverbatim
1.7 bertrand 57: *
58: * Arguments:
59: * ==========
60: *
61: *> \param[in] JOBA
62: *> \verbatim
63: *> JOBA is CHARACTER*1
64: *> Specifies the level of accuracy:
65: *> = 'C': This option works well (high relative accuracy) if A = B * D,
66: *> with well-conditioned B and arbitrary diagonal matrix D.
67: *> The accuracy cannot be spoiled by COLUMN scaling. The
68: *> accuracy of the computed output depends on the condition of
69: *> B, and the procedure aims at the best theoretical accuracy.
70: *> The relative error max_{i=1:N}|d sigma_i| / sigma_i is
71: *> bounded by f(M,N)*epsilon* cond(B), independent of D.
72: *> The input matrix is preprocessed with the QRF with column
73: *> pivoting. This initial preprocessing and preconditioning by
74: *> a rank revealing QR factorization is common for all values of
75: *> JOBA. Additional actions are specified as follows:
76: *> = 'E': Computation as with 'C' with an additional estimate of the
77: *> condition number of B. It provides a realistic error bound.
78: *> = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings
79: *> D1, D2, and well-conditioned matrix C, this option gives
80: *> higher accuracy than the 'C' option. If the structure of the
81: *> input matrix is not known, and relative accuracy is
82: *> desirable, then this option is advisable. The input matrix A
83: *> is preprocessed with QR factorization with FULL (row and
84: *> column) pivoting.
1.20 bertrand 85: *> = 'G': Computation as with 'F' with an additional estimate of the
1.7 bertrand 86: *> condition number of B, where A=D*B. If A has heavily weighted
87: *> rows, then using this condition number gives too pessimistic
88: *> error bound.
89: *> = 'A': Small singular values are the noise and the matrix is treated
1.16 bertrand 90: *> as numerically rank deficient. The error in the computed
1.7 bertrand 91: *> singular values is bounded by f(m,n)*epsilon*||A||.
92: *> The computed SVD A = U * S * V^t restores A up to
93: *> f(m,n)*epsilon*||A||.
94: *> This gives the procedure the licence to discard (set to zero)
95: *> all singular values below N*epsilon*||A||.
96: *> = 'R': Similar as in 'A'. Rank revealing property of the initial
97: *> QR factorization is used do reveal (using triangular factor)
98: *> a gap sigma_{r+1} < epsilon * sigma_r in which case the
99: *> numerical RANK is declared to be r. The SVD is computed with
100: *> absolute error bounds, but more accurately than with 'A'.
101: *> \endverbatim
102: *>
103: *> \param[in] JOBU
104: *> \verbatim
105: *> JOBU is CHARACTER*1
106: *> Specifies whether to compute the columns of U:
107: *> = 'U': N columns of U are returned in the array U.
108: *> = 'F': full set of M left sing. vectors is returned in the array U.
109: *> = 'W': U may be used as workspace of length M*N. See the description
110: *> of U.
111: *> = 'N': U is not computed.
112: *> \endverbatim
113: *>
114: *> \param[in] JOBV
115: *> \verbatim
116: *> JOBV is CHARACTER*1
117: *> Specifies whether to compute the matrix V:
118: *> = 'V': N columns of V are returned in the array V; Jacobi rotations
119: *> are not explicitly accumulated.
120: *> = 'J': N columns of V are returned in the array V, but they are
121: *> computed as the product of Jacobi rotations. This option is
122: *> allowed only if JOBU .NE. 'N', i.e. in computing the full SVD.
123: *> = 'W': V may be used as workspace of length N*N. See the description
124: *> of V.
125: *> = 'N': V is not computed.
126: *> \endverbatim
127: *>
128: *> \param[in] JOBR
129: *> \verbatim
130: *> JOBR is CHARACTER*1
131: *> Specifies the RANGE for the singular values. Issues the licence to
132: *> set to zero small positive singular values if they are outside
133: *> specified range. If A .NE. 0 is scaled so that the largest singular
134: *> value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues
135: *> the licence to kill columns of A whose norm in c*A is less than
1.20 bertrand 136: *> DSQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN,
1.7 bertrand 137: *> where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E').
138: *> = 'N': Do not kill small columns of c*A. This option assumes that
139: *> BLAS and QR factorizations and triangular solvers are
140: *> implemented to work in that range. If the condition of A
141: *> is greater than BIG, use DGESVJ.
142: *> = 'R': RESTRICTED range for sigma(c*A) is [DSQRT(SFMIN), DSQRT(BIG)]
143: *> (roughly, as described above). This option is recommended.
144: *> ~~~~~~~~~~~~~~~~~~~~~~~~~~~
145: *> For computing the singular values in the FULL range [SFMIN,BIG]
146: *> use DGESVJ.
147: *> \endverbatim
148: *>
149: *> \param[in] JOBT
150: *> \verbatim
151: *> JOBT is CHARACTER*1
152: *> If the matrix is square then the procedure may determine to use
153: *> transposed A if A^t seems to be better with respect to convergence.
154: *> If the matrix is not square, JOBT is ignored. This is subject to
155: *> changes in the future.
156: *> The decision is based on two values of entropy over the adjoint
157: *> orbit of A^t * A. See the descriptions of WORK(6) and WORK(7).
158: *> = 'T': transpose if entropy test indicates possibly faster
159: *> convergence of Jacobi process if A^t is taken as input. If A is
160: *> replaced with A^t, then the row pivoting is included automatically.
161: *> = 'N': do not speculate.
162: *> This option can be used to compute only the singular values, or the
163: *> full SVD (U, SIGMA and V). For only one set of singular vectors
164: *> (U or V), the caller should provide both U and V, as one of the
165: *> matrices is used as workspace if the matrix A is transposed.
166: *> The implementer can easily remove this constraint and make the
167: *> code more complicated. See the descriptions of U and V.
168: *> \endverbatim
169: *>
170: *> \param[in] JOBP
171: *> \verbatim
172: *> JOBP is CHARACTER*1
173: *> Issues the licence to introduce structured perturbations to drown
174: *> denormalized numbers. This licence should be active if the
175: *> denormals are poorly implemented, causing slow computation,
176: *> especially in cases of fast convergence (!). For details see [1,2].
177: *> For the sake of simplicity, this perturbations are included only
178: *> when the full SVD or only the singular values are requested. The
179: *> implementer/user can easily add the perturbation for the cases of
180: *> computing one set of singular vectors.
181: *> = 'P': introduce perturbation
182: *> = 'N': do not perturb
183: *> \endverbatim
184: *>
185: *> \param[in] M
186: *> \verbatim
187: *> M is INTEGER
188: *> The number of rows of the input matrix A. M >= 0.
189: *> \endverbatim
190: *>
191: *> \param[in] N
192: *> \verbatim
193: *> N is INTEGER
194: *> The number of columns of the input matrix A. M >= N >= 0.
195: *> \endverbatim
196: *>
197: *> \param[in,out] A
198: *> \verbatim
199: *> A is DOUBLE PRECISION array, dimension (LDA,N)
200: *> On entry, the M-by-N matrix A.
201: *> \endverbatim
202: *>
203: *> \param[in] LDA
204: *> \verbatim
205: *> LDA is INTEGER
206: *> The leading dimension of the array A. LDA >= max(1,M).
207: *> \endverbatim
208: *>
209: *> \param[out] SVA
210: *> \verbatim
211: *> SVA is DOUBLE PRECISION array, dimension (N)
212: *> On exit,
213: *> - For WORK(1)/WORK(2) = ONE: The singular values of A. During the
214: *> computation SVA contains Euclidean column norms of the
215: *> iterated matrices in the array A.
216: *> - For WORK(1) .NE. WORK(2): The singular values of A are
217: *> (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if
218: *> sigma_max(A) overflows or if small singular values have been
219: *> saved from underflow by scaling the input matrix A.
220: *> - If JOBR='R' then some of the singular values may be returned
221: *> as exact zeros obtained by "set to zero" because they are
222: *> below the numerical rank threshold or are denormalized numbers.
223: *> \endverbatim
224: *>
225: *> \param[out] U
226: *> \verbatim
227: *> U is DOUBLE PRECISION array, dimension ( LDU, N )
228: *> If JOBU = 'U', then U contains on exit the M-by-N matrix of
229: *> the left singular vectors.
230: *> If JOBU = 'F', then U contains on exit the M-by-M matrix of
231: *> the left singular vectors, including an ONB
232: *> of the orthogonal complement of the Range(A).
1.20 bertrand 233: *> If JOBU = 'W' .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N),
1.7 bertrand 234: *> then U is used as workspace if the procedure
235: *> replaces A with A^t. In that case, [V] is computed
236: *> in U as left singular vectors of A^t and then
237: *> copied back to the V array. This 'W' option is just
238: *> a reminder to the caller that in this case U is
239: *> reserved as workspace of length N*N.
1.14 bertrand 240: *> If JOBU = 'N' U is not referenced, unless JOBT='T'.
1.7 bertrand 241: *> \endverbatim
242: *>
243: *> \param[in] LDU
244: *> \verbatim
245: *> LDU is INTEGER
246: *> The leading dimension of the array U, LDU >= 1.
247: *> IF JOBU = 'U' or 'F' or 'W', then LDU >= M.
248: *> \endverbatim
249: *>
250: *> \param[out] V
251: *> \verbatim
252: *> V is DOUBLE PRECISION array, dimension ( LDV, N )
253: *> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
254: *> the right singular vectors;
1.20 bertrand 255: *> If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N),
1.7 bertrand 256: *> then V is used as workspace if the pprocedure
257: *> replaces A with A^t. In that case, [U] is computed
258: *> in V as right singular vectors of A^t and then
259: *> copied back to the U array. This 'W' option is just
260: *> a reminder to the caller that in this case V is
261: *> reserved as workspace of length N*N.
1.14 bertrand 262: *> If JOBV = 'N' V is not referenced, unless JOBT='T'.
1.7 bertrand 263: *> \endverbatim
264: *>
265: *> \param[in] LDV
266: *> \verbatim
267: *> LDV is INTEGER
268: *> The leading dimension of the array V, LDV >= 1.
269: *> If JOBV = 'V' or 'J' or 'W', then LDV >= N.
270: *> \endverbatim
271: *>
272: *> \param[out] WORK
273: *> \verbatim
1.18 bertrand 274: *> WORK is DOUBLE PRECISION array, dimension (LWORK)
1.20 bertrand 275: *> On exit, if N > 0 .AND. M > 0 (else not referenced),
1.7 bertrand 276: *> WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such
277: *> that SCALE*SVA(1:N) are the computed singular values
278: *> of A. (See the description of SVA().)
279: *> WORK(2) = See the description of WORK(1).
280: *> WORK(3) = SCONDA is an estimate for the condition number of
1.20 bertrand 281: *> column equilibrated A. (If JOBA = 'E' or 'G')
1.7 bertrand 282: *> SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
283: *> It is computed using DPOCON. It holds
284: *> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
285: *> where R is the triangular factor from the QRF of A.
286: *> However, if R is truncated and the numerical rank is
287: *> determined to be strictly smaller than N, SCONDA is
288: *> returned as -1, thus indicating that the smallest
289: *> singular values might be lost.
290: *>
291: *> If full SVD is needed, the following two condition numbers are
1.21 ! bertrand 292: *> useful for the analysis of the algorithm. They are provided for
1.7 bertrand 293: *> a developer/implementer who is familiar with the details of
294: *> the method.
295: *>
296: *> WORK(4) = an estimate of the scaled condition number of the
297: *> triangular factor in the first QR factorization.
298: *> WORK(5) = an estimate of the scaled condition number of the
299: *> triangular factor in the second QR factorization.
1.20 bertrand 300: *> The following two parameters are computed if JOBT = 'T'.
1.7 bertrand 301: *> They are provided for a developer/implementer who is familiar
302: *> with the details of the method.
303: *>
304: *> WORK(6) = the entropy of A^t*A :: this is the Shannon entropy
305: *> of diag(A^t*A) / Trace(A^t*A) taken as point in the
306: *> probability simplex.
307: *> WORK(7) = the entropy of A*A^t.
308: *> \endverbatim
309: *>
310: *> \param[in] LWORK
311: *> \verbatim
312: *> LWORK is INTEGER
313: *> Length of WORK to confirm proper allocation of work space.
314: *> LWORK depends on the job:
315: *>
1.20 bertrand 316: *> If only SIGMA is needed (JOBU = 'N', JOBV = 'N') and
317: *> -> .. no scaled condition estimate required (JOBE = 'N'):
1.7 bertrand 318: *> LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.
319: *> ->> For optimal performance (blocked code) the optimal value
320: *> is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal
321: *> block size for DGEQP3 and DGEQRF.
1.16 bertrand 322: *> In general, optimal LWORK is computed as
323: *> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 7).
1.7 bertrand 324: *> -> .. an estimate of the scaled condition number of A is
325: *> required (JOBA='E', 'G'). In this case, LWORK is the maximum
326: *> of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4*N,7).
1.16 bertrand 327: *> ->> For optimal performance (blocked code) the optimal value
1.7 bertrand 328: *> is LWORK >= max(2*M+N,3*N+(N+1)*NB, N*N+4*N, 7).
329: *> In general, the optimal length LWORK is computed as
1.16 bertrand 330: *> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF),
1.7 bertrand 331: *> N+N*N+LWORK(DPOCON),7).
332: *>
1.20 bertrand 333: *> If SIGMA and the right singular vectors are needed (JOBV = 'V'),
1.7 bertrand 334: *> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
335: *> -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
1.14 bertrand 336: *> where NB is the optimal block size for DGEQP3, DGEQRF, DGELQF,
1.7 bertrand 337: *> DORMLQ. In general, the optimal length LWORK is computed as
1.16 bertrand 338: *> LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON),
1.14 bertrand 339: *> N+LWORK(DGELQF), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)).
1.7 bertrand 340: *>
341: *> If SIGMA and the left singular vectors are needed
342: *> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
343: *> -> For optimal performance:
1.20 bertrand 344: *> if JOBU = 'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
345: *> if JOBU = 'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
1.7 bertrand 346: *> where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR.
347: *> In general, the optimal length LWORK is computed as
348: *> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON),
1.16 bertrand 349: *> 2*N+LWORK(DGEQRF), N+LWORK(DORMQR)).
1.20 bertrand 350: *> Here LWORK(DORMQR) equals N*NB (for JOBU = 'U') or
351: *> M*NB (for JOBU = 'F').
1.16 bertrand 352: *>
1.20 bertrand 353: *> If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and
354: *> -> if JOBV = 'V'
1.16 bertrand 355: *> the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N).
1.20 bertrand 356: *> -> if JOBV = 'J' the minimal requirement is
1.7 bertrand 357: *> LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6).
358: *> -> For optimal performance, LWORK should be additionally
359: *> larger than N+M*NB, where NB is the optimal block size
360: *> for DORMQR.
361: *> \endverbatim
362: *>
363: *> \param[out] IWORK
364: *> \verbatim
1.18 bertrand 365: *> IWORK is INTEGER array, dimension (M+3*N).
1.7 bertrand 366: *> On exit,
367: *> IWORK(1) = the numerical rank determined after the initial
368: *> QR factorization with pivoting. See the descriptions
369: *> of JOBA and JOBR.
370: *> IWORK(2) = the number of the computed nonzero singular values
371: *> IWORK(3) = if nonzero, a warning message:
1.20 bertrand 372: *> If IWORK(3) = 1 then some of the column norms of A
1.7 bertrand 373: *> were denormalized floats. The requested high accuracy
374: *> is not warranted by the data.
375: *> \endverbatim
376: *>
377: *> \param[out] INFO
378: *> \verbatim
379: *> INFO is INTEGER
1.20 bertrand 380: *> < 0: if INFO = -i, then the i-th argument had an illegal value.
381: *> = 0: successful exit;
382: *> > 0: DGEJSV did not converge in the maximal allowed number
383: *> of sweeps. The computed values may be inaccurate.
1.7 bertrand 384: *> \endverbatim
385: *
386: * Authors:
387: * ========
388: *
1.16 bertrand 389: *> \author Univ. of Tennessee
390: *> \author Univ. of California Berkeley
391: *> \author Univ. of Colorado Denver
392: *> \author NAG Ltd.
1.7 bertrand 393: *
1.10 bertrand 394: *> \ingroup doubleGEsing
1.7 bertrand 395: *
396: *> \par Further Details:
397: * =====================
398: *>
399: *> \verbatim
400: *>
401: *> DGEJSV implements a preconditioned Jacobi SVD algorithm. It uses DGEQP3,
402: *> DGEQRF, and DGELQF as preprocessors and preconditioners. Optionally, an
403: *> additional row pivoting can be used as a preprocessor, which in some
404: *> cases results in much higher accuracy. An example is matrix A with the
405: *> structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
406: *> diagonal matrices and C is well-conditioned matrix. In that case, complete
407: *> pivoting in the first QR factorizations provides accuracy dependent on the
408: *> condition number of C, and independent of D1, D2. Such higher accuracy is
409: *> not completely understood theoretically, but it works well in practice.
410: *> Further, if A can be written as A = B*D, with well-conditioned B and some
411: *> diagonal D, then the high accuracy is guaranteed, both theoretically and
412: *> in software, independent of D. For more details see [1], [2].
413: *> The computational range for the singular values can be the full range
414: *> ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
415: *> & LAPACK routines called by DGEJSV are implemented to work in that range.
416: *> If that is not the case, then the restriction for safe computation with
417: *> the singular values in the range of normalized IEEE numbers is that the
418: *> spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
419: *> overflow. This code (DGEJSV) is best used in this restricted range,
420: *> meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are
421: *> returned as zeros. See JOBR for details on this.
422: *> Further, this implementation is somewhat slower than the one described
423: *> in [1,2] due to replacement of some non-LAPACK components, and because
424: *> the choice of some tuning parameters in the iterative part (DGESVJ) is
425: *> left to the implementer on a particular machine.
426: *> The rank revealing QR factorization (in this code: DGEQP3) should be
427: *> implemented as in [3]. We have a new version of DGEQP3 under development
428: *> that is more robust than the current one in LAPACK, with a cleaner cut in
1.16 bertrand 429: *> rank deficient cases. It will be available in the SIGMA library [4].
430: *> If M is much larger than N, it is obvious that the initial QRF with
1.7 bertrand 431: *> column pivoting can be preprocessed by the QRF without pivoting. That
432: *> well known trick is not used in DGEJSV because in some cases heavy row
433: *> weighting can be treated with complete pivoting. The overhead in cases
434: *> M much larger than N is then only due to pivoting, but the benefits in
435: *> terms of accuracy have prevailed. The implementer/user can incorporate
436: *> this extra QRF step easily. The implementer can also improve data movement
437: *> (matrix transpose, matrix copy, matrix transposed copy) - this
438: *> implementation of DGEJSV uses only the simplest, naive data movement.
439: *> \endverbatim
440: *
441: *> \par Contributors:
442: * ==================
443: *>
444: *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
445: *
446: *> \par References:
447: * ================
448: *>
449: *> \verbatim
450: *>
451: *> [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
452: *> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
453: *> LAPACK Working note 169.
454: *> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
455: *> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
456: *> LAPACK Working note 170.
457: *> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
458: *> factorization software - a case study.
459: *> ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
460: *> LAPACK Working note 176.
461: *> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
462: *> QSVD, (H,K)-SVD computations.
463: *> Department of Mathematics, University of Zagreb, 2008.
464: *> \endverbatim
465: *
466: *> \par Bugs, examples and comments:
467: * =================================
468: *>
469: *> Please report all bugs and send interesting examples and/or comments to
470: *> drmac@math.hr. Thank you.
471: *>
472: * =====================================================================
1.1 bertrand 473: SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
1.6 bertrand 474: $ M, N, A, LDA, SVA, U, LDU, V, LDV,
475: $ WORK, LWORK, IWORK, INFO )
1.1 bertrand 476: *
1.21 ! bertrand 477: * -- LAPACK computational routine --
1.1 bertrand 478: * -- LAPACK is a software package provided by Univ. of Tennessee, --
479: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
480: *
481: * .. Scalar Arguments ..
482: IMPLICIT NONE
483: INTEGER INFO, LDA, LDU, LDV, LWORK, M, N
484: * ..
485: * .. Array Arguments ..
486: DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ),
1.6 bertrand 487: $ WORK( LWORK )
1.1 bertrand 488: INTEGER IWORK( * )
489: CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
490: * ..
491: *
1.6 bertrand 492: * ===========================================================================
1.1 bertrand 493: *
494: * .. Local Parameters ..
495: DOUBLE PRECISION ZERO, ONE
496: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
497: * ..
498: * .. Local Scalars ..
499: DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK,
1.6 bertrand 500: $ CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN, MAXPRJ, SCALEM,
501: $ SCONDA, SFMIN, SMALL, TEMP1, USCAL1, USCAL2, XSC
1.1 bertrand 502: INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING
503: LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LSVEC,
1.6 bertrand 504: $ L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN,
505: $ NOSCAL, ROWPIV, RSVEC, TRANSP
1.1 bertrand 506: * ..
507: * .. Intrinsic Functions ..
1.13 bertrand 508: INTRINSIC DABS, DLOG, MAX, MIN, DBLE, IDNINT, DSIGN, DSQRT
1.1 bertrand 509: * ..
510: * .. External Functions ..
511: DOUBLE PRECISION DLAMCH, DNRM2
512: INTEGER IDAMAX
513: LOGICAL LSAME
514: EXTERNAL IDAMAX, LSAME, DLAMCH, DNRM2
515: * ..
516: * .. External Subroutines ..
517: EXTERNAL DCOPY, DGELQF, DGEQP3, DGEQRF, DLACPY, DLASCL,
1.6 bertrand 518: $ DLASET, DLASSQ, DLASWP, DORGQR, DORMLQ,
519: $ DORMQR, DPOCON, DSCAL, DSWAP, DTRSM, XERBLA
1.1 bertrand 520: *
521: EXTERNAL DGESVJ
522: * ..
523: *
524: * Test the input arguments
525: *
526: LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
527: JRACC = LSAME( JOBV, 'J' )
528: RSVEC = LSAME( JOBV, 'V' ) .OR. JRACC
529: ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' )
530: L2RANK = LSAME( JOBA, 'R' )
531: L2ABER = LSAME( JOBA, 'A' )
532: ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' )
533: L2TRAN = LSAME( JOBT, 'T' )
534: L2KILL = LSAME( JOBR, 'R' )
535: DEFR = LSAME( JOBR, 'N' )
536: L2PERT = LSAME( JOBP, 'P' )
537: *
538: IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR.
1.6 bertrand 539: $ ERREST .OR. LSAME( JOBA, 'C' ) )) THEN
1.1 bertrand 540: INFO = - 1
541: ELSE IF ( .NOT.( LSVEC .OR. LSAME( JOBU, 'N' ) .OR.
1.6 bertrand 542: $ LSAME( JOBU, 'W' )) ) THEN
1.1 bertrand 543: INFO = - 2
544: ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR.
1.6 bertrand 545: $ LSAME( JOBV, 'W' )) .OR. ( JRACC .AND. (.NOT.LSVEC) ) ) THEN
1.1 bertrand 546: INFO = - 3
547: ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) ) THEN
548: INFO = - 4
549: ELSE IF ( .NOT. ( L2TRAN .OR. LSAME( JOBT, 'N' ) ) ) THEN
550: INFO = - 5
551: ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN
552: INFO = - 6
553: ELSE IF ( M .LT. 0 ) THEN
554: INFO = - 7
555: ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN
556: INFO = - 8
557: ELSE IF ( LDA .LT. M ) THEN
558: INFO = - 10
559: ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN
560: INFO = - 13
561: ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN
1.16 bertrand 562: INFO = - 15
1.1 bertrand 563: ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND.
1.13 bertrand 564: & (LWORK .LT. MAX(7,4*N+1,2*M+N))) .OR.
1.6 bertrand 565: & (.NOT.(LSVEC .OR. RSVEC) .AND. ERREST .AND.
1.13 bertrand 566: & (LWORK .LT. MAX(7,4*N+N*N,2*M+N))) .OR.
567: & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX(7,2*M+N,4*N+1)))
1.6 bertrand 568: & .OR.
1.13 bertrand 569: & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX(7,2*M+N,4*N+1)))
1.6 bertrand 570: & .OR.
1.16 bertrand 571: & (LSVEC .AND. RSVEC .AND. (.NOT.JRACC) .AND.
1.13 bertrand 572: & (LWORK.LT.MAX(2*M+N,6*N+2*N*N)))
1.6 bertrand 573: & .OR. (LSVEC .AND. RSVEC .AND. JRACC .AND.
1.13 bertrand 574: & LWORK.LT.MAX(2*M+N,4*N+N*N,2*N+N*N+6)))
1.1 bertrand 575: & THEN
576: INFO = - 17
577: ELSE
578: * #:)
579: INFO = 0
580: END IF
581: *
582: IF ( INFO .NE. 0 ) THEN
583: * #:(
584: CALL XERBLA( 'DGEJSV', - INFO )
1.6 bertrand 585: RETURN
1.1 bertrand 586: END IF
587: *
588: * Quick return for void matrix (Y3K safe)
589: * #:)
1.14 bertrand 590: IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
591: IWORK(1:3) = 0
592: WORK(1:7) = 0
593: RETURN
594: ENDIF
1.1 bertrand 595: *
596: * Determine whether the matrix U should be M x N or M x M
597: *
598: IF ( LSVEC ) THEN
599: N1 = N
600: IF ( LSAME( JOBU, 'F' ) ) N1 = M
601: END IF
602: *
603: * Set numerical parameters
604: *
605: *! NOTE: Make sure DLAMCH() does not fail on the target architecture.
606: *
607: EPSLN = DLAMCH('Epsilon')
608: SFMIN = DLAMCH('SafeMinimum')
609: SMALL = SFMIN / EPSLN
610: BIG = DLAMCH('O')
611: * BIG = ONE / SFMIN
612: *
613: * Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
614: *
615: *(!) If necessary, scale SVA() to protect the largest norm from
616: * overflow. It is possible that this scaling pushes the smallest
617: * column norm left from the underflow threshold (extreme case).
618: *
619: SCALEM = ONE / DSQRT(DBLE(M)*DBLE(N))
620: NOSCAL = .TRUE.
621: GOSCAL = .TRUE.
622: DO 1874 p = 1, N
623: AAPP = ZERO
1.4 bertrand 624: AAQQ = ONE
1.1 bertrand 625: CALL DLASSQ( M, A(1,p), 1, AAPP, AAQQ )
626: IF ( AAPP .GT. BIG ) THEN
627: INFO = - 9
628: CALL XERBLA( 'DGEJSV', -INFO )
629: RETURN
630: END IF
631: AAQQ = DSQRT(AAQQ)
632: IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL ) THEN
633: SVA(p) = AAPP * AAQQ
634: ELSE
635: NOSCAL = .FALSE.
636: SVA(p) = AAPP * ( AAQQ * SCALEM )
637: IF ( GOSCAL ) THEN
638: GOSCAL = .FALSE.
639: CALL DSCAL( p-1, SCALEM, SVA, 1 )
640: END IF
641: END IF
642: 1874 CONTINUE
643: *
644: IF ( NOSCAL ) SCALEM = ONE
645: *
646: AAPP = ZERO
647: AAQQ = BIG
648: DO 4781 p = 1, N
1.13 bertrand 649: AAPP = MAX( AAPP, SVA(p) )
650: IF ( SVA(p) .NE. ZERO ) AAQQ = MIN( AAQQ, SVA(p) )
1.1 bertrand 651: 4781 CONTINUE
652: *
653: * Quick return for zero M x N matrix
654: * #:)
655: IF ( AAPP .EQ. ZERO ) THEN
656: IF ( LSVEC ) CALL DLASET( 'G', M, N1, ZERO, ONE, U, LDU )
657: IF ( RSVEC ) CALL DLASET( 'G', N, N, ZERO, ONE, V, LDV )
658: WORK(1) = ONE
659: WORK(2) = ONE
660: IF ( ERREST ) WORK(3) = ONE
661: IF ( LSVEC .AND. RSVEC ) THEN
662: WORK(4) = ONE
663: WORK(5) = ONE
664: END IF
665: IF ( L2TRAN ) THEN
666: WORK(6) = ZERO
667: WORK(7) = ZERO
668: END IF
669: IWORK(1) = 0
670: IWORK(2) = 0
1.6 bertrand 671: IWORK(3) = 0
1.1 bertrand 672: RETURN
673: END IF
674: *
675: * Issue warning if denormalized column norms detected. Override the
676: * high relative accuracy request. Issue licence to kill columns
677: * (set them to zero) whose norm is less than sigma_max / BIG (roughly).
678: * #:(
679: WARNING = 0
680: IF ( AAQQ .LE. SFMIN ) THEN
681: L2RANK = .TRUE.
682: L2KILL = .TRUE.
683: WARNING = 1
684: END IF
685: *
686: * Quick return for one-column matrix
687: * #:)
688: IF ( N .EQ. 1 ) THEN
689: *
690: IF ( LSVEC ) THEN
691: CALL DLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR )
692: CALL DLACPY( 'A', M, 1, A, LDA, U, LDU )
693: * computing all M left singular vectors of the M x 1 matrix
694: IF ( N1 .NE. N ) THEN
695: CALL DGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,IERR )
696: CALL DORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,IERR )
697: CALL DCOPY( M, A(1,1), 1, U(1,1), 1 )
698: END IF
699: END IF
700: IF ( RSVEC ) THEN
701: V(1,1) = ONE
702: END IF
703: IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN
704: SVA(1) = SVA(1) / SCALEM
705: SCALEM = ONE
706: END IF
707: WORK(1) = ONE / SCALEM
708: WORK(2) = ONE
709: IF ( SVA(1) .NE. ZERO ) THEN
710: IWORK(1) = 1
711: IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN
712: IWORK(2) = 1
713: ELSE
714: IWORK(2) = 0
715: END IF
716: ELSE
717: IWORK(1) = 0
718: IWORK(2) = 0
719: END IF
1.14 bertrand 720: IWORK(3) = 0
1.1 bertrand 721: IF ( ERREST ) WORK(3) = ONE
722: IF ( LSVEC .AND. RSVEC ) THEN
723: WORK(4) = ONE
724: WORK(5) = ONE
725: END IF
726: IF ( L2TRAN ) THEN
727: WORK(6) = ZERO
728: WORK(7) = ZERO
729: END IF
730: RETURN
731: *
732: END IF
733: *
734: TRANSP = .FALSE.
735: L2TRAN = L2TRAN .AND. ( M .EQ. N )
736: *
737: AATMAX = -ONE
738: AATMIN = BIG
739: IF ( ROWPIV .OR. L2TRAN ) THEN
740: *
741: * Compute the row norms, needed to determine row pivoting sequence
742: * (in the case of heavily row weighted A, row pivoting is strongly
743: * advised) and to collect information needed to compare the
744: * structures of A * A^t and A^t * A (in the case L2TRAN.EQ..TRUE.).
745: *
746: IF ( L2TRAN ) THEN
747: DO 1950 p = 1, M
748: XSC = ZERO
1.4 bertrand 749: TEMP1 = ONE
1.1 bertrand 750: CALL DLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
751: * DLASSQ gets both the ell_2 and the ell_infinity norm
752: * in one pass through the vector
753: WORK(M+N+p) = XSC * SCALEM
754: WORK(N+p) = XSC * (SCALEM*DSQRT(TEMP1))
1.13 bertrand 755: AATMAX = MAX( AATMAX, WORK(N+p) )
756: IF (WORK(N+p) .NE. ZERO) AATMIN = MIN(AATMIN,WORK(N+p))
1.1 bertrand 757: 1950 CONTINUE
758: ELSE
759: DO 1904 p = 1, M
760: WORK(M+N+p) = SCALEM*DABS( A(p,IDAMAX(N,A(p,1),LDA)) )
1.13 bertrand 761: AATMAX = MAX( AATMAX, WORK(M+N+p) )
762: AATMIN = MIN( AATMIN, WORK(M+N+p) )
1.1 bertrand 763: 1904 CONTINUE
764: END IF
765: *
766: END IF
767: *
768: * For square matrix A try to determine whether A^t would be better
769: * input for the preconditioned Jacobi SVD, with faster convergence.
770: * The decision is based on an O(N) function of the vector of column
771: * and row norms of A, based on the Shannon entropy. This should give
772: * the right choice in most cases when the difference actually matters.
773: * It may fail and pick the slower converging side.
774: *
775: ENTRA = ZERO
776: ENTRAT = ZERO
777: IF ( L2TRAN ) THEN
778: *
779: XSC = ZERO
1.4 bertrand 780: TEMP1 = ONE
1.1 bertrand 781: CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
782: TEMP1 = ONE / TEMP1
783: *
784: ENTRA = ZERO
785: DO 1113 p = 1, N
786: BIG1 = ( ( SVA(p) / XSC )**2 ) * TEMP1
787: IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1)
788: 1113 CONTINUE
789: ENTRA = - ENTRA / DLOG(DBLE(N))
790: *
791: * Now, SVA().^2/Trace(A^t * A) is a point in the probability simplex.
792: * It is derived from the diagonal of A^t * A. Do the same with the
793: * diagonal of A * A^t, compute the entropy of the corresponding
794: * probability distribution. Note that A * A^t and A^t * A have the
795: * same trace.
796: *
797: ENTRAT = ZERO
798: DO 1114 p = N+1, N+M
799: BIG1 = ( ( WORK(p) / XSC )**2 ) * TEMP1
800: IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1)
801: 1114 CONTINUE
802: ENTRAT = - ENTRAT / DLOG(DBLE(M))
803: *
804: * Analyze the entropies and decide A or A^t. Smaller entropy
805: * usually means better input for the algorithm.
806: *
807: TRANSP = ( ENTRAT .LT. ENTRA )
808: *
809: * If A^t is better than A, transpose A.
810: *
811: IF ( TRANSP ) THEN
812: * In an optimal implementation, this trivial transpose
813: * should be replaced with faster transpose.
814: DO 1115 p = 1, N - 1
815: DO 1116 q = p + 1, N
816: TEMP1 = A(q,p)
817: A(q,p) = A(p,q)
818: A(p,q) = TEMP1
819: 1116 CONTINUE
820: 1115 CONTINUE
821: DO 1117 p = 1, N
822: WORK(M+N+p) = SVA(p)
823: SVA(p) = WORK(N+p)
824: 1117 CONTINUE
825: TEMP1 = AAPP
826: AAPP = AATMAX
827: AATMAX = TEMP1
828: TEMP1 = AAQQ
829: AAQQ = AATMIN
830: AATMIN = TEMP1
831: KILL = LSVEC
832: LSVEC = RSVEC
833: RSVEC = KILL
1.4 bertrand 834: IF ( LSVEC ) N1 = N
1.1 bertrand 835: *
836: ROWPIV = .TRUE.
837: END IF
838: *
839: END IF
840: * END IF L2TRAN
841: *
842: * Scale the matrix so that its maximal singular value remains less
843: * than DSQRT(BIG) -- the matrix is scaled so that its maximal column
844: * has Euclidean norm equal to DSQRT(BIG/N). The only reason to keep
845: * DSQRT(BIG) instead of BIG is the fact that DGEJSV uses LAPACK and
846: * BLAS routines that, in some implementations, are not capable of
847: * working in the full interval [SFMIN,BIG] and that they may provoke
848: * overflows in the intermediate results. If the singular values spread
849: * from SFMIN to BIG, then DGESVJ will compute them. So, in that case,
850: * one should use DGESVJ instead of DGEJSV.
851: *
852: BIG1 = DSQRT( BIG )
853: TEMP1 = DSQRT( BIG / DBLE(N) )
854: *
855: CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
856: IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN
857: AAQQ = ( AAQQ / AAPP ) * TEMP1
858: ELSE
859: AAQQ = ( AAQQ * TEMP1 ) / AAPP
860: END IF
861: TEMP1 = TEMP1 * SCALEM
862: CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR )
863: *
864: * To undo scaling at the end of this procedure, multiply the
865: * computed singular values with USCAL2 / USCAL1.
866: *
867: USCAL1 = TEMP1
868: USCAL2 = AAPP
869: *
870: IF ( L2KILL ) THEN
871: * L2KILL enforces computation of nonzero singular values in
872: * the restricted range of condition number of the initial A,
873: * sigma_max(A) / sigma_min(A) approx. DSQRT(BIG)/DSQRT(SFMIN).
874: XSC = DSQRT( SFMIN )
875: ELSE
876: XSC = SMALL
877: *
878: * Now, if the condition number of A is too big,
879: * sigma_max(A) / sigma_min(A) .GT. DSQRT(BIG/N) * EPSLN / SFMIN,
880: * as a precaution measure, the full SVD is computed using DGESVJ
881: * with accumulated Jacobi rotations. This provides numerically
882: * more robust computation, at the cost of slightly increased run
883: * time. Depending on the concrete implementation of BLAS and LAPACK
884: * (i.e. how they behave in presence of extreme ill-conditioning) the
885: * implementor may decide to remove this switch.
886: IF ( ( AAQQ.LT.DSQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN
887: JRACC = .TRUE.
888: END IF
889: *
890: END IF
891: IF ( AAQQ .LT. XSC ) THEN
892: DO 700 p = 1, N
893: IF ( SVA(p) .LT. XSC ) THEN
894: CALL DLASET( 'A', M, 1, ZERO, ZERO, A(1,p), LDA )
895: SVA(p) = ZERO
896: END IF
897: 700 CONTINUE
898: END IF
899: *
900: * Preconditioning using QR factorization with pivoting
901: *
902: IF ( ROWPIV ) THEN
903: * Optional row permutation (Bjoerck row pivoting):
904: * A result by Cox and Higham shows that the Bjoerck's
905: * row pivoting combined with standard column pivoting
906: * has similar effect as Powell-Reid complete pivoting.
907: * The ell-infinity norms of A are made nonincreasing.
908: DO 1952 p = 1, M - 1
909: q = IDAMAX( M-p+1, WORK(M+N+p), 1 ) + p - 1
910: IWORK(2*N+p) = q
911: IF ( p .NE. q ) THEN
912: TEMP1 = WORK(M+N+p)
913: WORK(M+N+p) = WORK(M+N+q)
914: WORK(M+N+q) = TEMP1
915: END IF
916: 1952 CONTINUE
917: CALL DLASWP( N, A, LDA, 1, M-1, IWORK(2*N+1), 1 )
918: END IF
919: *
920: * End of the preparation phase (scaling, optional sorting and
921: * transposing, optional flushing of small columns).
922: *
923: * Preconditioning
924: *
925: * If the full SVD is needed, the right singular vectors are computed
926: * from a matrix equation, and for that we need theoretical analysis
927: * of the Businger-Golub pivoting. So we use DGEQP3 as the first RR QRF.
928: * In all other cases the first RR QRF can be chosen by other criteria
929: * (eg speed by replacing global with restricted window pivoting, such
930: * as in SGEQPX from TOMS # 782). Good results will be obtained using
931: * SGEQPX with properly (!) chosen numerical parameters.
1.21 ! bertrand 932: * Any improvement of DGEQP3 improves overall performance of DGEJSV.
1.1 bertrand 933: *
934: * A * P1 = Q1 * [ R1^t 0]^t:
935: DO 1963 p = 1, N
936: * .. all columns are free columns
937: IWORK(p) = 0
938: 1963 CONTINUE
939: CALL DGEQP3( M,N,A,LDA, IWORK,WORK, WORK(N+1),LWORK-N, IERR )
940: *
941: * The upper triangular matrix R1 from the first QRF is inspected for
942: * rank deficiency and possibilities for deflation, or possible
943: * ill-conditioning. Depending on the user specified flag L2RANK,
944: * the procedure explores possibilities to reduce the numerical
945: * rank by inspecting the computed upper triangular factor. If
946: * L2RANK or L2ABER are up, then DGEJSV will compute the SVD of
947: * A + dA, where ||dA|| <= f(M,N)*EPSLN.
948: *
949: NR = 1
950: IF ( L2ABER ) THEN
951: * Standard absolute error bound suffices. All sigma_i with
952: * sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
1.20 bertrand 953: * aggressive enforcement of lower numerical rank by introducing a
1.1 bertrand 954: * backward error of the order of N*EPSLN*||A||.
955: TEMP1 = DSQRT(DBLE(N))*EPSLN
956: DO 3001 p = 2, N
957: IF ( DABS(A(p,p)) .GE. (TEMP1*DABS(A(1,1))) ) THEN
958: NR = NR + 1
959: ELSE
960: GO TO 3002
961: END IF
962: 3001 CONTINUE
963: 3002 CONTINUE
964: ELSE IF ( L2RANK ) THEN
1.20 bertrand 965: * .. similarly as above, only slightly more gentle (less aggressive).
1.1 bertrand 966: * Sudden drop on the diagonal of R1 is used as the criterion for
1.16 bertrand 967: * close-to-rank-deficient.
1.1 bertrand 968: TEMP1 = DSQRT(SFMIN)
969: DO 3401 p = 2, N
970: IF ( ( DABS(A(p,p)) .LT. (EPSLN*DABS(A(p-1,p-1))) ) .OR.
1.6 bertrand 971: $ ( DABS(A(p,p)) .LT. SMALL ) .OR.
972: $ ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402
1.1 bertrand 973: NR = NR + 1
974: 3401 CONTINUE
975: 3402 CONTINUE
976: *
977: ELSE
978: * The goal is high relative accuracy. However, if the matrix
979: * has high scaled condition number the relative accuracy is in
980: * general not feasible. Later on, a condition number estimator
981: * will be deployed to estimate the scaled condition number.
982: * Here we just remove the underflowed part of the triangular
983: * factor. This prevents the situation in which the code is
984: * working hard to get the accuracy not warranted by the data.
985: TEMP1 = DSQRT(SFMIN)
986: DO 3301 p = 2, N
987: IF ( ( DABS(A(p,p)) .LT. SMALL ) .OR.
1.6 bertrand 988: $ ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302
1.1 bertrand 989: NR = NR + 1
990: 3301 CONTINUE
991: 3302 CONTINUE
992: *
993: END IF
994: *
995: ALMORT = .FALSE.
996: IF ( NR .EQ. N ) THEN
997: MAXPRJ = ONE
998: DO 3051 p = 2, N
999: TEMP1 = DABS(A(p,p)) / SVA(IWORK(p))
1.13 bertrand 1000: MAXPRJ = MIN( MAXPRJ, TEMP1 )
1.1 bertrand 1001: 3051 CONTINUE
1002: IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE.
1003: END IF
1004: *
1005: *
1006: SCONDA = - ONE
1007: CONDR1 = - ONE
1008: CONDR2 = - ONE
1009: *
1010: IF ( ERREST ) THEN
1011: IF ( N .EQ. NR ) THEN
1012: IF ( RSVEC ) THEN
1013: * .. V is available as workspace
1014: CALL DLACPY( 'U', N, N, A, LDA, V, LDV )
1015: DO 3053 p = 1, N
1016: TEMP1 = SVA(IWORK(p))
1017: CALL DSCAL( p, ONE/TEMP1, V(1,p), 1 )
1018: 3053 CONTINUE
1019: CALL DPOCON( 'U', N, V, LDV, ONE, TEMP1,
1.6 bertrand 1020: $ WORK(N+1), IWORK(2*N+M+1), IERR )
1.1 bertrand 1021: ELSE IF ( LSVEC ) THEN
1022: * .. U is available as workspace
1023: CALL DLACPY( 'U', N, N, A, LDA, U, LDU )
1024: DO 3054 p = 1, N
1025: TEMP1 = SVA(IWORK(p))
1026: CALL DSCAL( p, ONE/TEMP1, U(1,p), 1 )
1027: 3054 CONTINUE
1028: CALL DPOCON( 'U', N, U, LDU, ONE, TEMP1,
1.6 bertrand 1029: $ WORK(N+1), IWORK(2*N+M+1), IERR )
1.1 bertrand 1030: ELSE
1031: CALL DLACPY( 'U', N, N, A, LDA, WORK(N+1), N )
1032: DO 3052 p = 1, N
1033: TEMP1 = SVA(IWORK(p))
1034: CALL DSCAL( p, ONE/TEMP1, WORK(N+(p-1)*N+1), 1 )
1035: 3052 CONTINUE
1036: * .. the columns of R are scaled to have unit Euclidean lengths.
1037: CALL DPOCON( 'U', N, WORK(N+1), N, ONE, TEMP1,
1.6 bertrand 1038: $ WORK(N+N*N+1), IWORK(2*N+M+1), IERR )
1.1 bertrand 1039: END IF
1040: SCONDA = ONE / DSQRT(TEMP1)
1041: * SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
1042: * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
1043: ELSE
1044: SCONDA = - ONE
1045: END IF
1046: END IF
1047: *
1048: L2PERT = L2PERT .AND. ( DABS( A(1,1)/A(NR,NR) ) .GT. DSQRT(BIG1) )
1049: * If there is no violent scaling, artificial perturbation is not needed.
1050: *
1051: * Phase 3:
1052: *
1053: IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
1054: *
1055: * Singular Values only
1056: *
1057: * .. transpose A(1:NR,1:N)
1.13 bertrand 1058: DO 1946 p = 1, MIN( N-1, NR )
1.1 bertrand 1059: CALL DCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 )
1060: 1946 CONTINUE
1061: *
1062: * The following two DO-loops introduce small relative perturbation
1063: * into the strict upper triangle of the lower triangular matrix.
1064: * Small entries below the main diagonal are also changed.
1065: * This modification is useful if the computing environment does not
1066: * provide/allow FLUSH TO ZERO underflow, for it prevents many
1067: * annoying denormalized numbers in case of strongly scaled matrices.
1068: * The perturbation is structured so that it does not introduce any
1069: * new perturbation of the singular values, and it does not destroy
1070: * the job done by the preconditioner.
1071: * The licence for this perturbation is in the variable L2PERT, which
1072: * should be .FALSE. if FLUSH TO ZERO underflow is active.
1073: *
1074: IF ( .NOT. ALMORT ) THEN
1075: *
1076: IF ( L2PERT ) THEN
1077: * XSC = DSQRT(SMALL)
1078: XSC = EPSLN / DBLE(N)
1079: DO 4947 q = 1, NR
1080: TEMP1 = XSC*DABS(A(q,q))
1081: DO 4949 p = 1, N
1082: IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) )
1.6 bertrand 1083: $ .OR. ( p .LT. q ) )
1084: $ A(p,q) = DSIGN( TEMP1, A(p,q) )
1.1 bertrand 1085: 4949 CONTINUE
1086: 4947 CONTINUE
1087: ELSE
1088: CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, A(1,2),LDA )
1089: END IF
1090: *
1091: * .. second preconditioning using the QR factorization
1092: *
1093: CALL DGEQRF( N,NR, A,LDA, WORK, WORK(N+1),LWORK-N, IERR )
1094: *
1095: * .. and transpose upper to lower triangular
1096: DO 1948 p = 1, NR - 1
1097: CALL DCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 )
1098: 1948 CONTINUE
1099: *
1100: END IF
1101: *
1102: * Row-cyclic Jacobi SVD algorithm with column pivoting
1103: *
1104: * .. again some perturbation (a "background noise") is added
1105: * to drown denormals
1106: IF ( L2PERT ) THEN
1107: * XSC = DSQRT(SMALL)
1108: XSC = EPSLN / DBLE(N)
1109: DO 1947 q = 1, NR
1110: TEMP1 = XSC*DABS(A(q,q))
1111: DO 1949 p = 1, NR
1112: IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) )
1.6 bertrand 1113: $ .OR. ( p .LT. q ) )
1114: $ A(p,q) = DSIGN( TEMP1, A(p,q) )
1.1 bertrand 1115: 1949 CONTINUE
1116: 1947 CONTINUE
1117: ELSE
1118: CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2), LDA )
1119: END IF
1120: *
1121: * .. and one-sided Jacobi rotations are started on a lower
1122: * triangular matrix (plus perturbation which is ignored in
1123: * the part which destroys triangular form (confusing?!))
1124: *
1125: CALL DGESVJ( 'L', 'NoU', 'NoV', NR, NR, A, LDA, SVA,
1.6 bertrand 1126: $ N, V, LDV, WORK, LWORK, INFO )
1.1 bertrand 1127: *
1128: SCALEM = WORK(1)
1129: NUMRANK = IDNINT(WORK(2))
1130: *
1131: *
1132: ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN
1133: *
1134: * -> Singular Values and Right Singular Vectors <-
1135: *
1136: IF ( ALMORT ) THEN
1137: *
1138: * .. in this case NR equals N
1139: DO 1998 p = 1, NR
1140: CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1141: 1998 CONTINUE
1142: CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1143: *
1144: CALL DGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA,
1.6 bertrand 1145: $ WORK, LWORK, INFO )
1.1 bertrand 1146: SCALEM = WORK(1)
1147: NUMRANK = IDNINT(WORK(2))
1148:
1149: ELSE
1150: *
1151: * .. two more QR factorizations ( one QRF is not enough, two require
1152: * accumulated product of Jacobi rotations, three are perfect )
1153: *
1154: CALL DLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1), LDA )
1155: CALL DGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N, IERR)
1156: CALL DLACPY( 'Lower', NR, NR, A, LDA, V, LDV )
1157: CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1158: CALL DGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1.6 bertrand 1159: $ LWORK-2*N, IERR )
1.1 bertrand 1160: DO 8998 p = 1, NR
1161: CALL DCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
1162: 8998 CONTINUE
1163: CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1164: *
1165: CALL DGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,
1.6 bertrand 1166: $ LDU, WORK(N+1), LWORK, INFO )
1.1 bertrand 1167: SCALEM = WORK(N+1)
1168: NUMRANK = IDNINT(WORK(N+2))
1169: IF ( NR .LT. N ) THEN
1170: CALL DLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1), LDV )
1171: CALL DLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1), LDV )
1172: CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1), LDV )
1173: END IF
1174: *
1175: CALL DORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK,
1.6 bertrand 1176: $ V, LDV, WORK(N+1), LWORK-N, IERR )
1.1 bertrand 1177: *
1178: END IF
1179: *
1180: DO 8991 p = 1, N
1181: CALL DCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA )
1182: 8991 CONTINUE
1183: CALL DLACPY( 'All', N, N, A, LDA, V, LDV )
1184: *
1185: IF ( TRANSP ) THEN
1186: CALL DLACPY( 'All', N, N, V, LDV, U, LDU )
1187: END IF
1188: *
1189: ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN
1190: *
1191: * .. Singular Values and Left Singular Vectors ..
1192: *
1193: * .. second preconditioning step to avoid need to accumulate
1194: * Jacobi rotations in the Jacobi iterations.
1195: DO 1965 p = 1, NR
1196: CALL DCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 )
1197: 1965 CONTINUE
1198: CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
1199: *
1200: CALL DGEQRF( N, NR, U, LDU, WORK(N+1), WORK(2*N+1),
1.6 bertrand 1201: $ LWORK-2*N, IERR )
1.1 bertrand 1202: *
1203: DO 1967 p = 1, NR - 1
1204: CALL DCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 )
1205: 1967 CONTINUE
1206: CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
1207: *
1208: CALL DGESVJ( 'Lower', 'U', 'N', NR,NR, U, LDU, SVA, NR, A,
1.6 bertrand 1209: $ LDA, WORK(N+1), LWORK-N, INFO )
1.1 bertrand 1210: SCALEM = WORK(N+1)
1211: NUMRANK = IDNINT(WORK(N+2))
1212: *
1213: IF ( NR .LT. M ) THEN
1214: CALL DLASET( 'A', M-NR, NR,ZERO, ZERO, U(NR+1,1), LDU )
1215: IF ( NR .LT. N1 ) THEN
1216: CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1), LDU )
1217: CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1), LDU )
1218: END IF
1219: END IF
1220: *
1221: CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
1.6 bertrand 1222: $ LDU, WORK(N+1), LWORK-N, IERR )
1.1 bertrand 1223: *
1224: IF ( ROWPIV )
1.6 bertrand 1225: $ CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1.1 bertrand 1226: *
1227: DO 1974 p = 1, N1
1228: XSC = ONE / DNRM2( M, U(1,p), 1 )
1229: CALL DSCAL( M, XSC, U(1,p), 1 )
1230: 1974 CONTINUE
1231: *
1232: IF ( TRANSP ) THEN
1233: CALL DLACPY( 'All', N, N, U, LDU, V, LDV )
1234: END IF
1235: *
1236: ELSE
1237: *
1238: * .. Full SVD ..
1239: *
1240: IF ( .NOT. JRACC ) THEN
1241: *
1242: IF ( .NOT. ALMORT ) THEN
1243: *
1244: * Second Preconditioning Step (QRF [with pivoting])
1245: * Note that the composition of TRANSPOSE, QRF and TRANSPOSE is
1246: * equivalent to an LQF CALL. Since in many libraries the QRF
1247: * seems to be better optimized than the LQF, we do explicit
1248: * transpose and use the QRF. This is subject to changes in an
1249: * optimized implementation of DGEJSV.
1250: *
1251: DO 1968 p = 1, NR
1252: CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1253: 1968 CONTINUE
1254: *
1255: * .. the following two loops perturb small entries to avoid
1256: * denormals in the second QR factorization, where they are
1257: * as good as zeros. This is done to avoid painfully slow
1258: * computation with denormals. The relative size of the perturbation
1259: * is a parameter that can be changed by the implementer.
1260: * This perturbation device will be obsolete on machines with
1261: * properly implemented arithmetic.
1262: * To switch it off, set L2PERT=.FALSE. To remove it from the
1263: * code, remove the action under L2PERT=.TRUE., leave the ELSE part.
1264: * The following two loops should be blocked and fused with the
1265: * transposed copy above.
1266: *
1267: IF ( L2PERT ) THEN
1268: XSC = DSQRT(SMALL)
1269: DO 2969 q = 1, NR
1270: TEMP1 = XSC*DABS( V(q,q) )
1271: DO 2968 p = 1, N
1272: IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 )
1.6 bertrand 1273: $ .OR. ( p .LT. q ) )
1274: $ V(p,q) = DSIGN( TEMP1, V(p,q) )
1275: IF ( p .LT. q ) V(p,q) = - V(p,q)
1.1 bertrand 1276: 2968 CONTINUE
1277: 2969 CONTINUE
1278: ELSE
1279: CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1280: END IF
1281: *
1282: * Estimate the row scaled condition number of R1
1283: * (If R1 is rectangular, N > NR, then the condition number
1284: * of the leading NR x NR submatrix is estimated.)
1285: *
1286: CALL DLACPY( 'L', NR, NR, V, LDV, WORK(2*N+1), NR )
1287: DO 3950 p = 1, NR
1288: TEMP1 = DNRM2(NR-p+1,WORK(2*N+(p-1)*NR+p),1)
1289: CALL DSCAL(NR-p+1,ONE/TEMP1,WORK(2*N+(p-1)*NR+p),1)
1290: 3950 CONTINUE
1291: CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1,
1.6 bertrand 1292: $ WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR)
1.1 bertrand 1293: CONDR1 = ONE / DSQRT(TEMP1)
1.20 bertrand 1294: * .. here need a second opinion on the condition number
1.1 bertrand 1295: * .. then assume worst case scenario
1296: * R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
1297: * more conservative <=> CONDR1 .LT. DSQRT(DBLE(N))
1298: *
1299: COND_OK = DSQRT(DBLE(NR))
1300: *[TP] COND_OK is a tuning parameter.
1301:
1302: IF ( CONDR1 .LT. COND_OK ) THEN
1303: * .. the second QRF without pivoting. Note: in an optimized
1304: * implementation, this QRF should be implemented as the QRF
1305: * of a lower triangular matrix.
1306: * R1^t = Q2 * R2
1307: CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1.6 bertrand 1308: $ LWORK-2*N, IERR )
1.1 bertrand 1309: *
1310: IF ( L2PERT ) THEN
1311: XSC = DSQRT(SMALL)/EPSLN
1312: DO 3959 p = 2, NR
1313: DO 3958 q = 1, p - 1
1.13 bertrand 1314: TEMP1 = XSC * MIN(DABS(V(p,p)),DABS(V(q,q)))
1.1 bertrand 1315: IF ( DABS(V(q,p)) .LE. TEMP1 )
1.6 bertrand 1316: $ V(q,p) = DSIGN( TEMP1, V(q,p) )
1.1 bertrand 1317: 3958 CONTINUE
1318: 3959 CONTINUE
1319: END IF
1320: *
1321: IF ( NR .NE. N )
1.6 bertrand 1322: $ CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N )
1.1 bertrand 1323: * .. save ...
1324: *
1325: * .. this transposed copy should be better than naive
1326: DO 1969 p = 1, NR - 1
1327: CALL DCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 )
1328: 1969 CONTINUE
1329: *
1330: CONDR2 = CONDR1
1331: *
1332: ELSE
1333: *
1334: * .. ill-conditioned case: second QRF with pivoting
1.20 bertrand 1335: * Note that windowed pivoting would be equally good
1.1 bertrand 1336: * numerically, and more run-time efficient. So, in
1337: * an optimal implementation, the next call to DGEQP3
1338: * should be replaced with eg. CALL SGEQPX (ACM TOMS #782)
1339: * with properly (carefully) chosen parameters.
1340: *
1341: * R1^t * P2 = Q2 * R2
1342: DO 3003 p = 1, NR
1343: IWORK(N+p) = 0
1344: 3003 CONTINUE
1345: CALL DGEQP3( N, NR, V, LDV, IWORK(N+1), WORK(N+1),
1.6 bertrand 1346: $ WORK(2*N+1), LWORK-2*N, IERR )
1.1 bertrand 1347: ** CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1.6 bertrand 1348: ** $ LWORK-2*N, IERR )
1.1 bertrand 1349: IF ( L2PERT ) THEN
1350: XSC = DSQRT(SMALL)
1351: DO 3969 p = 2, NR
1352: DO 3968 q = 1, p - 1
1.13 bertrand 1353: TEMP1 = XSC * MIN(DABS(V(p,p)),DABS(V(q,q)))
1.1 bertrand 1354: IF ( DABS(V(q,p)) .LE. TEMP1 )
1.6 bertrand 1355: $ V(q,p) = DSIGN( TEMP1, V(q,p) )
1.1 bertrand 1356: 3968 CONTINUE
1357: 3969 CONTINUE
1358: END IF
1359: *
1360: CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N )
1361: *
1362: IF ( L2PERT ) THEN
1363: XSC = DSQRT(SMALL)
1364: DO 8970 p = 2, NR
1365: DO 8971 q = 1, p - 1
1.13 bertrand 1366: TEMP1 = XSC * MIN(DABS(V(p,p)),DABS(V(q,q)))
1.1 bertrand 1367: V(p,q) = - DSIGN( TEMP1, V(q,p) )
1368: 8971 CONTINUE
1369: 8970 CONTINUE
1370: ELSE
1371: CALL DLASET( 'L',NR-1,NR-1,ZERO,ZERO,V(2,1),LDV )
1372: END IF
1373: * Now, compute R2 = L3 * Q3, the LQ factorization.
1374: CALL DGELQF( NR, NR, V, LDV, WORK(2*N+N*NR+1),
1.6 bertrand 1375: $ WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR )
1.1 bertrand 1376: * .. and estimate the condition number
1377: CALL DLACPY( 'L',NR,NR,V,LDV,WORK(2*N+N*NR+NR+1),NR )
1378: DO 4950 p = 1, NR
1379: TEMP1 = DNRM2( p, WORK(2*N+N*NR+NR+p), NR )
1380: CALL DSCAL( p, ONE/TEMP1, WORK(2*N+N*NR+NR+p), NR )
1381: 4950 CONTINUE
1382: CALL DPOCON( 'L',NR,WORK(2*N+N*NR+NR+1),NR,ONE,TEMP1,
1.6 bertrand 1383: $ WORK(2*N+N*NR+NR+NR*NR+1),IWORK(M+2*N+1),IERR )
1.1 bertrand 1384: CONDR2 = ONE / DSQRT(TEMP1)
1385: *
1386: IF ( CONDR2 .GE. COND_OK ) THEN
1387: * .. save the Householder vectors used for Q3
1.20 bertrand 1388: * (this overwrites the copy of R2, as it will not be
1.1 bertrand 1389: * needed in this branch, but it does not overwritte the
1390: * Huseholder vectors of Q2.).
1391: CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N )
1392: * .. and the rest of the information on Q3 is in
1393: * WORK(2*N+N*NR+1:2*N+N*NR+N)
1394: END IF
1395: *
1396: END IF
1397: *
1398: IF ( L2PERT ) THEN
1399: XSC = DSQRT(SMALL)
1400: DO 4968 q = 2, NR
1401: TEMP1 = XSC * V(q,q)
1402: DO 4969 p = 1, q - 1
1403: * V(p,q) = - DSIGN( TEMP1, V(q,p) )
1404: V(p,q) = - DSIGN( TEMP1, V(p,q) )
1405: 4969 CONTINUE
1406: 4968 CONTINUE
1407: ELSE
1408: CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV )
1409: END IF
1410: *
1411: * Second preconditioning finished; continue with Jacobi SVD
1412: * The input matrix is lower trinagular.
1413: *
1414: * Recover the right singular vectors as solution of a well
1415: * conditioned triangular matrix equation.
1416: *
1417: IF ( CONDR1 .LT. COND_OK ) THEN
1418: *
1419: CALL DGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U,
1.6 bertrand 1420: $ LDU,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,INFO )
1.1 bertrand 1421: SCALEM = WORK(2*N+N*NR+NR+1)
1422: NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
1423: DO 3970 p = 1, NR
1424: CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 )
1425: CALL DSCAL( NR, SVA(p), V(1,p), 1 )
1426: 3970 CONTINUE
1427:
1428: * .. pick the right matrix equation and solve it
1429: *
1.6 bertrand 1430: IF ( NR .EQ. N ) THEN
1.1 bertrand 1431: * :)) .. best case, R1 is inverted. The solution of this matrix
1432: * equation is Q2*V2 = the product of the Jacobi rotations
1433: * used in DGESVJ, premultiplied with the orthogonal matrix
1434: * from the second QR factorization.
1435: CALL DTRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,LDV )
1436: ELSE
1437: * .. R1 is well conditioned, but non-square. Transpose(R2)
1438: * is inverted to get the product of the Jacobi rotations
1439: * used in DGESVJ. The Q-factor from the second QR
1440: * factorization is then built in explicitly.
1441: CALL DTRSM('L','U','T','N',NR,NR,ONE,WORK(2*N+1),
1.6 bertrand 1442: $ N,V,LDV)
1.1 bertrand 1443: IF ( NR .LT. N ) THEN
1444: CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
1445: CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
1446: CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
1447: END IF
1448: CALL DORMQR('L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1.6 bertrand 1449: $ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
1.1 bertrand 1450: END IF
1451: *
1452: ELSE IF ( CONDR2 .LT. COND_OK ) THEN
1453: *
1454: * :) .. the input matrix A is very likely a relative of
1455: * the Kahan matrix :)
1456: * The matrix R2 is inverted. The solution of the matrix equation
1457: * is Q3^T*V3 = the product of the Jacobi rotations (appplied to
1458: * the lower triangular L3 from the LQ factorization of
1459: * R2=L3*Q3), pre-multiplied with the transposed Q3.
1460: CALL DGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
1.6 bertrand 1461: $ LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
1.1 bertrand 1462: SCALEM = WORK(2*N+N*NR+NR+1)
1463: NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
1464: DO 3870 p = 1, NR
1465: CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 )
1466: CALL DSCAL( NR, SVA(p), U(1,p), 1 )
1467: 3870 CONTINUE
1468: CALL DTRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,LDU)
1469: * .. apply the permutation from the second QR factorization
1470: DO 873 q = 1, NR
1471: DO 872 p = 1, NR
1472: WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
1473: 872 CONTINUE
1474: DO 874 p = 1, NR
1475: U(p,q) = WORK(2*N+N*NR+NR+p)
1476: 874 CONTINUE
1477: 873 CONTINUE
1478: IF ( NR .LT. N ) THEN
1479: CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
1480: CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
1481: CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
1482: END IF
1483: CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1.6 bertrand 1484: $ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
1.1 bertrand 1485: ELSE
1486: * Last line of defense.
1487: * #:( This is a rather pathological case: no scaled condition
1488: * improvement after two pivoted QR factorizations. Other
1489: * possibility is that the rank revealing QR factorization
1490: * or the condition estimator has failed, or the COND_OK
1491: * is set very close to ONE (which is unnecessary). Normally,
1492: * this branch should never be executed, but in rare cases of
1493: * failure of the RRQR or condition estimator, the last line of
1494: * defense ensures that DGEJSV completes the task.
1495: * Compute the full SVD of L3 using DGESVJ with explicit
1496: * accumulation of Jacobi rotations.
1497: CALL DGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
1.6 bertrand 1498: $ LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
1.1 bertrand 1499: SCALEM = WORK(2*N+N*NR+NR+1)
1500: NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
1501: IF ( NR .LT. N ) THEN
1502: CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
1503: CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
1504: CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
1505: END IF
1506: CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1.6 bertrand 1507: $ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
1.1 bertrand 1508: *
1509: CALL DORMLQ( 'L', 'T', NR, NR, NR, WORK(2*N+1), N,
1.6 bertrand 1510: $ WORK(2*N+N*NR+1), U, LDU, WORK(2*N+N*NR+NR+1),
1511: $ LWORK-2*N-N*NR-NR, IERR )
1.1 bertrand 1512: DO 773 q = 1, NR
1513: DO 772 p = 1, NR
1514: WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
1515: 772 CONTINUE
1516: DO 774 p = 1, NR
1517: U(p,q) = WORK(2*N+N*NR+NR+p)
1518: 774 CONTINUE
1519: 773 CONTINUE
1520: *
1521: END IF
1522: *
1523: * Permute the rows of V using the (column) permutation from the
1524: * first QRF. Also, scale the columns to make them unit in
1525: * Euclidean norm. This applies to all cases.
1526: *
1527: TEMP1 = DSQRT(DBLE(N)) * EPSLN
1528: DO 1972 q = 1, N
1529: DO 972 p = 1, N
1530: WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
1531: 972 CONTINUE
1532: DO 973 p = 1, N
1533: V(p,q) = WORK(2*N+N*NR+NR+p)
1534: 973 CONTINUE
1535: XSC = ONE / DNRM2( N, V(1,q), 1 )
1536: IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1.6 bertrand 1537: $ CALL DSCAL( N, XSC, V(1,q), 1 )
1.1 bertrand 1538: 1972 CONTINUE
1539: * At this moment, V contains the right singular vectors of A.
1540: * Next, assemble the left singular vector matrix U (M x N).
1541: IF ( NR .LT. M ) THEN
1542: CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
1543: IF ( NR .LT. N1 ) THEN
1544: CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
1545: CALL DLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),LDU)
1546: END IF
1547: END IF
1548: *
1549: * The Q matrix from the first QRF is built into the left singular
1550: * matrix U. This applies to all cases.
1551: *
1552: CALL DORMQR( 'Left', 'No_Tr', M, N1, N, A, LDA, WORK, U,
1.6 bertrand 1553: $ LDU, WORK(N+1), LWORK-N, IERR )
1.1 bertrand 1554:
1555: * The columns of U are normalized. The cost is O(M*N) flops.
1556: TEMP1 = DSQRT(DBLE(M)) * EPSLN
1557: DO 1973 p = 1, NR
1558: XSC = ONE / DNRM2( M, U(1,p), 1 )
1559: IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1.6 bertrand 1560: $ CALL DSCAL( M, XSC, U(1,p), 1 )
1.1 bertrand 1561: 1973 CONTINUE
1562: *
1563: * If the initial QRF is computed with row pivoting, the left
1564: * singular vectors must be adjusted.
1565: *
1566: IF ( ROWPIV )
1.6 bertrand 1567: $ CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1.1 bertrand 1568: *
1569: ELSE
1570: *
1571: * .. the initial matrix A has almost orthogonal columns and
1572: * the second QRF is not needed
1573: *
1574: CALL DLACPY( 'Upper', N, N, A, LDA, WORK(N+1), N )
1575: IF ( L2PERT ) THEN
1576: XSC = DSQRT(SMALL)
1577: DO 5970 p = 2, N
1578: TEMP1 = XSC * WORK( N + (p-1)*N + p )
1579: DO 5971 q = 1, p - 1
1580: WORK(N+(q-1)*N+p)=-DSIGN(TEMP1,WORK(N+(p-1)*N+q))
1581: 5971 CONTINUE
1582: 5970 CONTINUE
1583: ELSE
1584: CALL DLASET( 'Lower',N-1,N-1,ZERO,ZERO,WORK(N+2),N )
1585: END IF
1586: *
1587: CALL DGESVJ( 'Upper', 'U', 'N', N, N, WORK(N+1), N, SVA,
1.6 bertrand 1588: $ N, U, LDU, WORK(N+N*N+1), LWORK-N-N*N, INFO )
1.1 bertrand 1589: *
1590: SCALEM = WORK(N+N*N+1)
1591: NUMRANK = IDNINT(WORK(N+N*N+2))
1592: DO 6970 p = 1, N
1593: CALL DCOPY( N, WORK(N+(p-1)*N+1), 1, U(1,p), 1 )
1594: CALL DSCAL( N, SVA(p), WORK(N+(p-1)*N+1), 1 )
1595: 6970 CONTINUE
1596: *
1597: CALL DTRSM( 'Left', 'Upper', 'NoTrans', 'No UD', N, N,
1.6 bertrand 1598: $ ONE, A, LDA, WORK(N+1), N )
1.1 bertrand 1599: DO 6972 p = 1, N
1600: CALL DCOPY( N, WORK(N+p), N, V(IWORK(p),1), LDV )
1601: 6972 CONTINUE
1602: TEMP1 = DSQRT(DBLE(N))*EPSLN
1603: DO 6971 p = 1, N
1604: XSC = ONE / DNRM2( N, V(1,p), 1 )
1605: IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1.6 bertrand 1606: $ CALL DSCAL( N, XSC, V(1,p), 1 )
1.1 bertrand 1607: 6971 CONTINUE
1608: *
1609: * Assemble the left singular vector matrix U (M x N).
1610: *
1611: IF ( N .LT. M ) THEN
1.4 bertrand 1612: CALL DLASET( 'A', M-N, N, ZERO, ZERO, U(N+1,1), LDU )
1.1 bertrand 1613: IF ( N .LT. N1 ) THEN
1614: CALL DLASET( 'A',N, N1-N, ZERO, ZERO, U(1,N+1),LDU )
1.4 bertrand 1615: CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),LDU )
1.1 bertrand 1616: END IF
1617: END IF
1618: CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
1.6 bertrand 1619: $ LDU, WORK(N+1), LWORK-N, IERR )
1.1 bertrand 1620: TEMP1 = DSQRT(DBLE(M))*EPSLN
1621: DO 6973 p = 1, N1
1622: XSC = ONE / DNRM2( M, U(1,p), 1 )
1623: IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1.6 bertrand 1624: $ CALL DSCAL( M, XSC, U(1,p), 1 )
1.1 bertrand 1625: 6973 CONTINUE
1626: *
1627: IF ( ROWPIV )
1.6 bertrand 1628: $ CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1.1 bertrand 1629: *
1630: END IF
1631: *
1632: * end of the >> almost orthogonal case << in the full SVD
1633: *
1634: ELSE
1635: *
1636: * This branch deploys a preconditioned Jacobi SVD with explicitly
1637: * accumulated rotations. It is included as optional, mainly for
1.20 bertrand 1638: * experimental purposes. It does perform well, and can also be used.
1.1 bertrand 1639: * In this implementation, this branch will be automatically activated
1640: * if the condition number sigma_max(A) / sigma_min(A) is predicted
1641: * to be greater than the overflow threshold. This is because the
1642: * a posteriori computation of the singular vectors assumes robust
1643: * implementation of BLAS and some LAPACK procedures, capable of working
1644: * in presence of extreme values. Since that is not always the case, ...
1645: *
1646: DO 7968 p = 1, NR
1647: CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1648: 7968 CONTINUE
1649: *
1650: IF ( L2PERT ) THEN
1651: XSC = DSQRT(SMALL/EPSLN)
1652: DO 5969 q = 1, NR
1653: TEMP1 = XSC*DABS( V(q,q) )
1654: DO 5968 p = 1, N
1655: IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 )
1.6 bertrand 1656: $ .OR. ( p .LT. q ) )
1657: $ V(p,q) = DSIGN( TEMP1, V(p,q) )
1658: IF ( p .LT. q ) V(p,q) = - V(p,q)
1.1 bertrand 1659: 5968 CONTINUE
1660: 5969 CONTINUE
1661: ELSE
1662: CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1663: END IF
1664:
1665: CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1.6 bertrand 1666: $ LWORK-2*N, IERR )
1.1 bertrand 1667: CALL DLACPY( 'L', N, NR, V, LDV, WORK(2*N+1), N )
1668: *
1669: DO 7969 p = 1, NR
1670: CALL DCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 )
1671: 7969 CONTINUE
1672:
1673: IF ( L2PERT ) THEN
1674: XSC = DSQRT(SMALL/EPSLN)
1675: DO 9970 q = 2, NR
1676: DO 9971 p = 1, q - 1
1.13 bertrand 1677: TEMP1 = XSC * MIN(DABS(U(p,p)),DABS(U(q,q)))
1.1 bertrand 1678: U(p,q) = - DSIGN( TEMP1, U(q,p) )
1679: 9971 CONTINUE
1680: 9970 CONTINUE
1681: ELSE
1682: CALL DLASET('U', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
1683: END IF
1684:
1685: CALL DGESVJ( 'G', 'U', 'V', NR, NR, U, LDU, SVA,
1.6 bertrand 1686: $ N, V, LDV, WORK(2*N+N*NR+1), LWORK-2*N-N*NR, INFO )
1.1 bertrand 1687: SCALEM = WORK(2*N+N*NR+1)
1688: NUMRANK = IDNINT(WORK(2*N+N*NR+2))
1689:
1690: IF ( NR .LT. N ) THEN
1691: CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
1692: CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
1693: CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
1694: END IF
1695:
1696: CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1.6 bertrand 1697: $ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
1.1 bertrand 1698: *
1699: * Permute the rows of V using the (column) permutation from the
1700: * first QRF. Also, scale the columns to make them unit in
1701: * Euclidean norm. This applies to all cases.
1702: *
1703: TEMP1 = DSQRT(DBLE(N)) * EPSLN
1704: DO 7972 q = 1, N
1705: DO 8972 p = 1, N
1706: WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
1707: 8972 CONTINUE
1708: DO 8973 p = 1, N
1709: V(p,q) = WORK(2*N+N*NR+NR+p)
1710: 8973 CONTINUE
1711: XSC = ONE / DNRM2( N, V(1,q), 1 )
1712: IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1.6 bertrand 1713: $ CALL DSCAL( N, XSC, V(1,q), 1 )
1.1 bertrand 1714: 7972 CONTINUE
1715: *
1716: * At this moment, V contains the right singular vectors of A.
1717: * Next, assemble the left singular vector matrix U (M x N).
1718: *
1.4 bertrand 1719: IF ( NR .LT. M ) THEN
1720: CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
1721: IF ( NR .LT. N1 ) THEN
1722: CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),LDU )
1723: CALL DLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),LDU )
1.1 bertrand 1724: END IF
1725: END IF
1726: *
1727: CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
1.6 bertrand 1728: $ LDU, WORK(N+1), LWORK-N, IERR )
1.1 bertrand 1729: *
1730: IF ( ROWPIV )
1.6 bertrand 1731: $ CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1.1 bertrand 1732: *
1733: *
1734: END IF
1735: IF ( TRANSP ) THEN
1736: * .. swap U and V because the procedure worked on A^t
1737: DO 6974 p = 1, N
1738: CALL DSWAP( N, U(1,p), 1, V(1,p), 1 )
1739: 6974 CONTINUE
1740: END IF
1741: *
1742: END IF
1743: * end of the full SVD
1744: *
1745: * Undo scaling, if necessary (and possible)
1746: *
1747: IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
1748: CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
1749: USCAL1 = ONE
1750: USCAL2 = ONE
1751: END IF
1752: *
1753: IF ( NR .LT. N ) THEN
1754: DO 3004 p = NR+1, N
1755: SVA(p) = ZERO
1756: 3004 CONTINUE
1757: END IF
1758: *
1759: WORK(1) = USCAL2 * SCALEM
1760: WORK(2) = USCAL1
1761: IF ( ERREST ) WORK(3) = SCONDA
1762: IF ( LSVEC .AND. RSVEC ) THEN
1763: WORK(4) = CONDR1
1764: WORK(5) = CONDR2
1765: END IF
1766: IF ( L2TRAN ) THEN
1767: WORK(6) = ENTRA
1768: WORK(7) = ENTRAT
1769: END IF
1770: *
1771: IWORK(1) = NR
1772: IWORK(2) = NUMRANK
1773: IWORK(3) = WARNING
1774: *
1775: RETURN
1776: * ..
1777: * .. END OF DGEJSV
1778: * ..
1779: END
1780: *
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgejsv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack