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