dgesvj - N matrix A, where M >= N
SUBROUTINE DGESVJ(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO) CHARACTER*1 JOBA, JOBU, JOBV INTEGER INFO, LDA, LDV, LWORK, M, MV, N DOUBLE PRECISION A(LDA,*), SVA(N), V(LDV,*), WORK(LWORK) SUBROUTINE DGESVJ_64(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO) CHARACTER*1 JOBA, JOBU, JOBV INTEGER*8 INFO, LDA, LDV, LWORK, M, MV, N DOUBLE PRECISION A(LDA,*), SVA(N), V(LDV,*), WORK(LWORK) F95 INTERFACE SUBROUTINE GESVJ(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO) INTEGER :: M, N, LDA, MV, LDV, LWORK, INFO CHARACTER(LEN=1) :: JOBA, JOBU, JOBV REAL(8), DIMENSION(:,:) :: A, V REAL(8), DIMENSION(:) :: SVA, WORK SUBROUTINE GESVJ_64(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO) INTEGER(8) :: M, N, LDA, MV, LDV, LWORK, INFO CHARACTER(LEN=1) :: JOBA, JOBU, JOBV REAL(8), DIMENSION(:,:) :: A, V REAL(8), DIMENSION(:) :: SVA, WORK C INTERFACE #include <sunperf.h> void dgesvj (char joba, char jobu, char jobv, int m, int n, double *a, int lda, double *sva, int mv, double *v, int ldv, double *work, int lwork, int *info); void dgesvj_64 (char joba, char jobu, char jobv, long m, long n, double *a, long lda, double *sva, long mv, double *v, long ldv, dou- ble *work, long lwork, long *info);
Oracle Solaris Studio Performance Library dgesvj(3P)
NAME
dgesvj - compute 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)
CHARACTER*1 JOBA, JOBU, JOBV
INTEGER INFO, LDA, LDV, LWORK, M, MV, N
DOUBLE PRECISION A(LDA,*), SVA(N), V(LDV,*), WORK(LWORK)
SUBROUTINE DGESVJ_64(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV,
WORK, LWORK, INFO)
CHARACTER*1 JOBA, JOBU, JOBV
INTEGER*8 INFO, LDA, LDV, LWORK, M, MV, N
DOUBLE PRECISION A(LDA,*), SVA(N), V(LDV,*), WORK(LWORK)
F95 INTERFACE
SUBROUTINE GESVJ(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK,
LWORK, INFO)
INTEGER :: M, N, LDA, MV, LDV, LWORK, INFO
CHARACTER(LEN=1) :: JOBA, JOBU, JOBV
REAL(8), DIMENSION(:,:) :: A, V
REAL(8), DIMENSION(:) :: SVA, WORK
SUBROUTINE GESVJ_64(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV,
WORK, LWORK, INFO)
INTEGER(8) :: M, N, LDA, MV, LDV, LWORK, INFO
CHARACTER(LEN=1) :: JOBA, JOBU, JOBV
REAL(8), DIMENSION(:,:) :: A, V
REAL(8), DIMENSION(:) :: SVA, WORK
C INTERFACE
#include <sunperf.h>
void dgesvj (char joba, char jobu, char jobv, int m, int n, double *a,
int lda, double *sva, int mv, double *v, int ldv, double
*work, int lwork, int *info);
void dgesvj_64 (char joba, char jobu, char jobv, long m, long n, double
*a, long lda, double *sva, long mv, double *v, long ldv, dou-
ble *work, long lwork, long *info);
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 orthogo-
nal 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.
ARGUMENTS
JOBA (input)
JOBA is 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)
JOBU is CHARACTER*1
Specifies whether to compute the left singular vectors (col-
umns 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 vectors. 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 orthogo-
nality 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)
JOBV is 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)
M is INTEGER
The number of rows of the input matrix A.
1/DLAMCH('E') > M >= 0.
N (input)
N is INTEGER
The number of columns of the input matrix A.
M >= N >= 0.
A (input/output)
A is DOUBLE PRECISION 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 corresponding to
underflowed or zero singular values are not computed. 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 com-
puted columns of U may not be orthogonal up to TOL. The out-
put
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)
LDA is INTEGER
The leading dimension of the array A.
LDA >= max(1,M).
SVA (output)
SVA is DOUBLE PRECISION 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 fac-
tored representation is due to the fact that some of the sin-
gular values of A might underflow or overflow.
If INFO .GT. 0 :
the procedure DGESVJ did not converge in the given number of
iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
MV (input)
MV is 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)
V is DOUBLE PRECISION 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)
LDV is 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/output)
WORK is DOUBLE PRECISION array, dimension max(4,M+N).
On entry :
If JOBU .EQ. 'C' :
WORK(1) = CTOL, where CTOL defines the threshold for conver-
gence.
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 EPS.
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 singu-
lar values that are larger than the underflow 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 converge, as it can be used to estimate whether the out-
put 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 (input)
LWORK is INTEGER
length of WORK, WORK >= MAX(6,M+N)
INFO (output)
INFO is 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.
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 denormalized singular values can be computed with the correspond-
ing gradual loss of accurate digits.
[1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling.
SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
[2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
singular value decomposition on a vector computer.
SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
[3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
[4] Z. Drmac: Implementation of Jacobi rotations for accurate singular
value computation in floating point arithmetic.
SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
[5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm
I.
SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
LAPACK Working note 169.
[6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm
II.
SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
LAPACK Working note 170.
[7] Z. Drmac: SIGMA - mathematical software library for accurate SVD,
PSV,
QSVD, (H,K)-SVD computations.
Department of Mathematics, University of Zagreb, 2008.
7 Nov 2015 dgesvj(3P)