1SGESVJ(1)LAPACK routine (version 3.2) SGESVJ(1)
2
6 SGESVJ - computes the singular value decomposition (SVD) of a real M-
7 by-N matrix A, where M >= N
8
10 SUBROUTINE SGESVJ( 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 REAL A( LDA, * ), SVA( N ), V( LDV, * ),
21 + WORK( LWORK )
23
25 SGESVJ 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=SQRT(M), EPS=SLAMCH('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, If JOBU .EQ. 'U' .OR.
118 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 SLAMCH('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=SQRT(M)*EPS (default); or TOL=CTOL*EPS
128 (JOBU.EQ.'C'), see the description of JOBU. If INFO .GT. 0,
129 the procedure SGESVJ 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':
136 If INFO .EQ. 0 : Note that the left singular vectors are 'for
137 free' in the one-sided Jacobi SVD algorithm. However, if only
138 the singular values are needed, the level of numerical orthogo‐
139 nality of U is not an issue and iterations are stopped when the
140 columns of the iterated matrix are numerically orthogonal up to
141 approximately M*EPS. Thus, on exit, A contains the columns of U
142 scaled with the corresponding singular values. If INFO .GT. 0
143 : the procedure SGESVJ did not converge in the given number of
144 iterations (sweeps).
145 LDA (input) INTEGER
147 The leading dimension of the array A. LDA >= max(1,M).
148 SVA (workspace/output) REAL array, dimension (N)
150 On exit, 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 SGESVJ 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 SGESVJ 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, If JOBU .EQ. 'C' : WORK(1) = CTOL, where CTOL defines
181 the threshold for convergence. The process stops if all col‐
182 umns of A are mutually orthogonal up to CTOL*EPS,
183 EPS=SLAMCH('E'). It is required that CTOL >= ONE, i.e. it is
184 not allowed to force the routine to obtain orthogonality below
185 EPSILON. On exit, WORK(1) = SCALE is the scaling factor such
186 that SCALE*SVA(1:N) are the computed singular vcalues of A.
187 (See description of SVA().) WORK(2) = NINT(WORK(2)) is the
188 number of the computed nonzero singular values. WORK(3) =
189 NINT(WORK(3)) is the number of the computed singular values
190 that are larger than the underflow threshold. WORK(4) =
191 NINT(WORK(4)) is the number of sweeps of Jacobi rotations
192 needed for numerical convergence. WORK(5) = max_{i.NE.j}
193 |COS(A(:,i),A(:,j))| in the last sweep. This is useful infor‐
194 mation in cases when SGESVJ did not converge, as it can be used
195 to estimate whether the output is stil useful and for post fes‐
196 tum analysis. WORK(6) = the largest absolute value over all
197 sines of the Jacobi rotation angles in the last sweep. It can
198 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 : SGESVJ 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 SGESVJ(1)