DGESVJ(1LAPACK routine (version 3.2)DGESVJ(1)NAME
DGESVJ - computes the singular value decomposition (SVD) of a real M-
by-N matrix A, where M >= N
SYNOPSIS
SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
+ LDV, WORK, LWORK, INFO )
IMPLICIT NONE
INTEGER INFO, LDA, LDV, LWORK, M, MV, N
CHARACTER*1 JOBA, JOBU, JOBV
DOUBLE PRECISION A( LDA, * ), SVA( N ), V( LDV, * ),
+ WORK( LWORK )
PURPOSE
DGESVJ computes the singular value decomposition (SVD) of a real M-by-N
matrix A, where M >= N. The SVD of A is written as
[++] [xx] [x0] [xx]
A = U * SIGMA * V^t, [++] = [xx] * [ox] * [xx]
[++] [xx]
where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
matrix, and V is an N-by-N orthogonal matrix. The diagonal elements of
SIGMA are the singular values of A. The columns of U and V are the left
and the right singular vectors of A, respectively.
Further Details
The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
rotations. The rotations are implemented as fast scaled rotations of
Anda and Park [1]. In the case of underflow of the Jacobi angle, a mod‐
ified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
column interchanges of de Rijk [2]. The relative accuracy of the com‐
puted singular values and the accuracy of the computed singular vectors
(in angle metric) is as guaranteed by the theory of Demmel and Veselic
[3]. The condition number that determines the accuracy in the full
rank case is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
spectral condition number. The best performance of this Jacobi SVD pro‐
cedure is achieved if used in an accelerated version of Drmac and
Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
Some tunning parameters (marked with [TP]) are available for the imple‐
menter.
The computational range for the nonzero singular values is the machine
number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even denor‐
malized singular values can be computed with the corresponding gradual
loss of accurate digits.
Contributors
~~~~~~~~~~~~
Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
References
~~~~~~~~~~
SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
singular value decomposition on a vector computer.
SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
value computation in floating point arithmetic.
SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
LAPACK Working note 169.
SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
LAPACK Working note 170.
QSVD, (H,K)-SVD computations.
Department of Mathematics, University of Zagreb, 2008. Bugs, Exam‐
ples and Comments
~~~~~~~~~~~~~~~~~~~~~~~~~~~
Please report all bugs and send interesting test examples and comments
to drmac@math.hr. Thank you.
ARGUMENTS
JOBA (input) CHARACTER* 1
Specifies the structure of A. = 'L': The input matrix A is
lower triangular;
= 'U': The input matrix A is upper triangular;
= 'G': The input matrix A is general M-by-N matrix, M >= N.
JOBU (input) CHARACTER*1
Specifies whether to compute the left singular vectors (columns
of U):
= 'U': The left singular vectors corresponding to the nonzero
singular values are computed and returned in the leading col‐
umns of A. See more details in the description of A. The
default numerical orthogonality threshold is set to approxi‐
mately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E'). = 'C':
Analogous to JOBU='U', except that user can control the level
of numerical orthogonality of the computed left singular vec‐
tors. TOL can be set to TOL = CTOL*EPS, where CTOL is given on
input in the array WORK. No CTOL smaller than ONE is allowed.
CTOL greater than 1 / EPS is meaningless. The option 'C' can be
used if M*EPS is satisfactory orthogonality of the computed
left singular vectors, so CTOL=M could save few sweeps of
Jacobi rotations. See the descriptions of A and WORK(1). =
'N': The matrix U is not computed. However, see the description
of A.
JOBV (input) CHARACTER*1
Specifies whether to compute the right singular vectors, that
is, the matrix V:
= 'V' : the matrix V is computed and returned in the array V
= 'A' : the Jacobi rotations are applied to the MV-by-N array
V. In other words, the right singular vector matrix V is not
computed explicitly, instead it is applied to an MV-by-N matrix
initially stored in the first MV rows of V. = 'N' : the matrix
V is not computed and the array V is not referenced
M (input) INTEGER
The number of rows of the input matrix A. M >= 0.
N (input) INTEGER
The number of columns of the input matrix A. M >= N >= 0.
A (input/output) REAL array, dimension (LDA,N)
On entry, the M-by-N matrix A. On exit :
If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' :
If INFO .EQ. 0 : RANKA orthonormal columns of U are returned in
the leading RANKA columns of the array A. Here RANKA <= N is
the number of computed singular values of A that are above the
underflow threshold DLAMCH('S'). The singular vectors corre‐
sponding to underflowed or zero singular values are not com‐
puted. The value of RANKA is returned in the array WORK as
RANKA=NINT(WORK(2)). Also see the descriptions of SVA and WORK.
The computed columns of U are mutually numerically orthogonal
up to approximately TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS
(JOBU.EQ.'C'), see the description of JOBU. If INFO .GT. 0 :
the procedure DGESVJ did not converge in the given number of
iterations (sweeps). In that case, the computed columns of U
may not be orthogonal up to TOL. The output U (stored in A),
SIGMA (given by the computed singular values in SVA(1:N)) and V
is still a decomposition of the input matrix A in the sense
that the residual ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
If JOBU .EQ. 'N' : If INFO .EQ. 0 : Note that the left singular
vectors are 'for free' in the one-sided Jacobi SVD algorithm.
However, if only the singular values are needed, the level of
numerical orthogonality of U is not an issue and iterations are
stopped when the columns of the iterated matrix are numerically
orthogonal up to approximately M*EPS. Thus, on exit, A contains
the columns of U scaled with the corresponding singular values.
If INFO .GT. 0 : the procedure DGESVJ did not converge in the
given number of iterations (sweeps).
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
SVA (workspace/output) REAL array, dimension (N)
On exit :
If INFO .EQ. 0 :
depending on the value SCALE = WORK(1), we have:
If SCALE .EQ. ONE :
SVA(1:N) contains the computed singular values of A. During
the computation SVA contains the Euclidean column norms of the
iterated matrices in the array A. If SCALE .NE. ONE :
The singular values of A are SCALE*SVA(1:N), and this factored
representation is due to the fact that some of the singular
values of A might underflow or overflow. If INFO .GT. 0 : the
procedure DGESVJ did not converge in the given number of itera‐
tions (sweeps) and SCALE*SVA(1:N) may not be accurate.
MV (input) INTEGER
If JOBV .EQ. 'A', then the product of Jacobi rotations in
DGESVJ is applied to the first MV rows of V. See the descrip‐
tion of JOBV.
V (input/output) REAL array, dimension (LDV,N)
If JOBV = 'V', then V contains on exit the N-by-N matrix of the
right singular vectors; If JOBV = 'A', then V contains the
product of the computed right singular vector matrix and the
initial matrix in the array V. If JOBV = 'N', then V is not
referenced.
LDV (input) INTEGER
The leading dimension of the array V, LDV .GE. 1. If JOBV .EQ.
'V', then LDV .GE. max(1,N). If JOBV .EQ. 'A', then LDV .GE.
max(1,MV) .
WORK (input/workspace/output) REAL array, dimension max(4,M+N).
On entry :
If JOBU .EQ. 'C' : WORK(1) = CTOL, where CTOL defines the
threshold for convergence. The process stops if all columns of
A are mutually orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). It
is required that CTOL >= ONE, i.e. it is not allowed to force
the routine to obtain orthogonality below EPSILON. On exit :
WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
are the computed singular values of A. (See description of
SVA().) WORK(2) = NINT(WORK(2)) is the number of the computed
nonzero singular values. WORK(3) = NINT(WORK(3)) is the number
of the computed singular values that are larger than the under‐
flow threshold. WORK(4) = NINT(WORK(4)) is the number of
sweeps of Jacobi rotations needed for numerical convergence.
WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
This is useful information in cases when DGESVJ did not con‐
verge, as it can be used to estimate whether the output is stil
useful and for post festum analysis. WORK(6) = the largest
absolute value over all sines of the Jacobi rotation angles in
the last sweep. It can be useful for a post festum analysis.
LWORK length of WORK, WORK >= MAX(6,M+N)
INFO (output) INTEGER
= 0 : successful exit.
< 0 : if INFO = -i, then the i-th argument had an illegal value
> 0 : DGESVJ did not converge in the maximal allowed number
(30) of sweeps. The output may still be useful. See the
description of WORK.
LAPACK routine (version 3.2) November 2008 DGESVJ(1)
