1DGESVJ(1)LAPACK routine (version 3.2) DGESVJ(1)
2
6 DGESVJ - computes the singular value decomposition (SVD) of a real M-
7 by-N matrix A, where M >= N
8
10 SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
11 + LDV, WORK, LWORK, INFO )
13 IMPLICIT NONE
15 INTEGER INFO, LDA, LDV, LWORK, M, MV, N
17 CHARACTER*1 JOBA, JOBU, JOBV
19 DOUBLE PRECISION A( LDA, * ), SVA( N ), V( LDV, * ),
21 + WORK( LWORK )
23
25 DGESVJ computes the singular value decomposition (SVD) of a real M-by-N
26 matrix A, where M >= N. The SVD of A is written as
27 [++] [xx] [x0] [xx]
28 A = U * SIGMA * V^t, [++] = [xx] * [ox] * [xx]
29 [++] [xx]
30 where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
31 matrix, and V is an N-by-N orthogonal matrix. The diagonal elements of
32 SIGMA are the singular values of A. The columns of U and V are the left
33 and the right singular vectors of A, respectively.
34 Further Details
35 The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
36 rotations. The rotations are implemented as fast scaled rotations of
37 Anda and Park [1]. In the case of underflow of the Jacobi angle, a mod‐
38 ified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
39 column interchanges of de Rijk [2]. The relative accuracy of the com‐
40 puted singular values and the accuracy of the computed singular vectors
41 (in angle metric) is as guaranteed by the theory of Demmel and Veselic
42 [3]. The condition number that determines the accuracy in the full
43 rank case is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
44 spectral condition number. The best performance of this Jacobi SVD pro‐
45 cedure is achieved if used in an accelerated version of Drmac and
46 Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
47 Some tunning parameters (marked with [TP]) are available for the imple‐
48 menter.
49 The computational range for the nonzero singular values is the machine
50 number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even denor‐
51 malized singular values can be computed with the corresponding gradual
52 loss of accurate digits.
53 Contributors
54 ~~~~~~~~~~~~
55 Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
56 References
57 ~~~~~~~~~~
58 SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
59 singular value decomposition on a vector computer.
60 SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
61 value computation in floating point arithmetic.
62 SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
63 SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
64 LAPACK Working note 169.
65 SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
66 LAPACK Working note 170.
67 QSVD, (H,K)-SVD computations.
68 Department of Mathematics, University of Zagreb, 2008. Bugs, Exam‐
69 ples and Comments
70 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
71 Please report all bugs and send interesting test examples and comments
72 to drmac@math.hr. Thank you.
73
75 JOBA (input) CHARACTER* 1
76 Specifies the structure of A. = 'L': The input matrix A is
77 lower triangular;
78 = 'U': The input matrix A is upper triangular;
79 = 'G': The input matrix A is general M-by-N matrix, M >= N.
80 JOBU (input) CHARACTER*1
82 Specifies whether to compute the left singular vectors (columns
83 of U):
84 = 'U': The left singular vectors corresponding to the nonzero
85 singular values are computed and returned in the leading col‐
86 umns of A. See more details in the description of A. The
87 default numerical orthogonality threshold is set to approxi‐
88 mately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E'). = 'C':
89 Analogous to JOBU='U', except that user can control the level
90 of numerical orthogonality of the computed left singular vec‐
91 tors. TOL can be set to TOL = CTOL*EPS, where CTOL is given on
92 input in the array WORK. No CTOL smaller than ONE is allowed.
93 CTOL greater than 1 / EPS is meaningless. The option 'C' can be
94 used if M*EPS is satisfactory orthogonality of the computed
95 left singular vectors, so CTOL=M could save few sweeps of
96 Jacobi rotations. See the descriptions of A and WORK(1). =
97 'N': The matrix U is not computed. However, see the description
98 of A.
99 JOBV (input) CHARACTER*1
101 Specifies whether to compute the right singular vectors, that
102 is, the matrix V:
103 = 'V' : the matrix V is computed and returned in the array V
104 = 'A' : the Jacobi rotations are applied to the MV-by-N array
105 V. In other words, the right singular vector matrix V is not
106 computed explicitly, instead it is applied to an MV-by-N matrix
107 initially stored in the first MV rows of V. = 'N' : the matrix
108 V is not computed and the array V is not referenced
109 M (input) INTEGER
111 The number of rows of the input matrix A. M >= 0.
112 N (input) INTEGER
114 The number of columns of the input matrix A. M >= N >= 0.
115 A (input/output) REAL array, dimension (LDA,N)
117 On entry, the M-by-N matrix A. On exit :
118 If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' :
119 If INFO .EQ. 0 : RANKA orthonormal columns of U are returned in
120 the leading RANKA columns of the array A. Here RANKA <= N is
121 the number of computed singular values of A that are above the
122 underflow threshold DLAMCH('S'). The singular vectors corre‐
123 sponding to underflowed or zero singular values are not com‐
124 puted. The value of RANKA is returned in the array WORK as
125 RANKA=NINT(WORK(2)). Also see the descriptions of SVA and WORK.
126 The computed columns of U are mutually numerically orthogonal
127 up to approximately TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS
128 (JOBU.EQ.'C'), see the description of JOBU. If INFO .GT. 0 :
129 the procedure DGESVJ did not converge in the given number of
130 iterations (sweeps). In that case, the computed columns of U
131 may not be orthogonal up to TOL. The output U (stored in A),
132 SIGMA (given by the computed singular values in SVA(1:N)) and V
133 is still a decomposition of the input matrix A in the sense
134 that the residual ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
135 If JOBU .EQ. 'N' : If INFO .EQ. 0 : Note that the left singular
136 vectors are 'for free' in the one-sided Jacobi SVD algorithm.
137 However, if only the singular values are needed, the level of
138 numerical orthogonality of U is not an issue and iterations are
139 stopped when the columns of the iterated matrix are numerically
140 orthogonal up to approximately M*EPS. Thus, on exit, A contains
141 the columns of U scaled with the corresponding singular values.
142 If INFO .GT. 0 : the procedure DGESVJ did not converge in the
143 given number of iterations (sweeps).
144 LDA (input) INTEGER
146 The leading dimension of the array A. LDA >= max(1,M).
147 SVA (workspace/output) REAL array, dimension (N)
149 On exit :
150 If INFO .EQ. 0 :
151 depending on the value SCALE = WORK(1), we have:
152 If SCALE .EQ. ONE :
153 SVA(1:N) contains the computed singular values of A. During
154 the computation SVA contains the Euclidean column norms of the
155 iterated matrices in the array A. If SCALE .NE. ONE :
156 The singular values of A are SCALE*SVA(1:N), and this factored
157 representation is due to the fact that some of the singular
158 values of A might underflow or overflow. If INFO .GT. 0 : the
159 procedure DGESVJ did not converge in the given number of itera‐
160 tions (sweeps) and SCALE*SVA(1:N) may not be accurate.
161 MV (input) INTEGER
163 If JOBV .EQ. 'A', then the product of Jacobi rotations in
164 DGESVJ is applied to the first MV rows of V. See the descrip‐
165 tion of JOBV.
166 V (input/output) REAL array, dimension (LDV,N)
168 If JOBV = 'V', then V contains on exit the N-by-N matrix of the
169 right singular vectors; If JOBV = 'A', then V contains the
170 product of the computed right singular vector matrix and the
171 initial matrix in the array V. If JOBV = 'N', then V is not
172 referenced.
173 LDV (input) INTEGER
175 The leading dimension of the array V, LDV .GE. 1. If JOBV .EQ.
176 'V', then LDV .GE. max(1,N). If JOBV .EQ. 'A', then LDV .GE.
177 max(1,MV) .
178 WORK (input/workspace/output) REAL array, dimension max(4,M+N).
180 On entry :
181 If JOBU .EQ. 'C' : WORK(1) = CTOL, where CTOL defines the
182 threshold for convergence. The process stops if all columns of
183 A are mutually orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). It
184 is required that CTOL >= ONE, i.e. it is not allowed to force
185 the routine to obtain orthogonality below EPSILON. On exit :
186 WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
187 are the computed singular values of A. (See description of
188 SVA().) WORK(2) = NINT(WORK(2)) is the number of the computed
189 nonzero singular values. WORK(3) = NINT(WORK(3)) is the number
190 of the computed singular values that are larger than the under‐
191 flow threshold. WORK(4) = NINT(WORK(4)) is the number of
192 sweeps of Jacobi rotations needed for numerical convergence.
193 WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
194 This is useful information in cases when DGESVJ did not con‐
195 verge, as it can be used to estimate whether the output is stil
196 useful and for post festum analysis. WORK(6) = the largest
197 absolute value over all sines of the Jacobi rotation angles in
198 the last sweep. It can be useful for a post festum analysis.
199 LWORK length of WORK, WORK >= MAX(6,M+N)
201 INFO (output) INTEGER
203 = 0 : successful exit.
204 < 0 : if INFO = -i, then the i-th argument had an illegal value
205 > 0 : DGESVJ did not converge in the maximal allowed number
206 (30) of sweeps. The output may still be useful. See the
207 description of WORK.
208 LAPACK routine (version 3.2) November 2008 DGESVJ(1)