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