Singular Value Decomposition of a General Matrix

Calling Sequence

CALL DSVDC	(DA, LDA, N, P, DSVALS, DE, DLSVEC, LDL, DRSVEC, LDR, DWORK, JOB, INFO)
CALL SSVDC	(SA, LDA, N, P, SSVALS, SE, SLSVEC, LDL, SRSVEC, LDR, SWORK, JOB, INFO)
CALL ZSVDC	(ZA, LDA, N, P, ZSVALS, ZE, ZLSVEC, LDL, ZRSVEC, LDR, ZWORK, JOB, INFO)
CALL CSVDC	(CA, LDA, N, P, CSVALS, CE, CLSVEC, LDL, CRSVEC, LDR, CWORK, JOB, INFO)
void dsvdc	(double *da, int lda, int n, int p, double *dsvals, double *de, double *dlsvec, int ldl, double *drsvec, int ldr, int job, int *info)
void ssvdc	(float *sa, int lda, int n, int p, float *ssvals, float *se, float *slsvec, int ldl, float *srsvec, int ldr, int job, int *info)
void zsvdc	(doublecomplex *za, int lda, int n, int p, doublecomplex *zsvals, doublecomplex *ze, doublecomplex *zlsvec, int ldl, doublecomplex *zrsvec, int ldr, int job, int *info)
void csvdc	(complex *ca, int lda, int n, int p, complex *csvals, complex *ce, complex *clsvec, int ldl, complex *crsvec, int ldr, int job, int *info)

Arguments

xA	Matrix A.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA max(1,N).
N	Number of rows of the matrix A. N 0.
P	Number of columns of the matrix A. P 0.
xSVALS	On exit, the singular values of A arranged in descending order of magnitude.
xE	On exit, normally contains zero, but see INFO below.
xLSVEC	On exit, the matrix of left singular vectors; not referenced if the a digit of JOB = 0.
LDL	Leading dimension of LSVEC as specified in a dimension or type statement. LDL max(1,N).
xRSVEC	On exit, the matrix of right singular vectors; not referenced if the b digit of JOB = 0.
LDR	Leading dimension of RSVEC as specified in a dimension or type statement. LDR max(1,P).
xWORK	Scratch array with a dimension of N.
JOB	Integer in the form ab; determines operation subroutine will perform:
	a = 0	do not compute the left singular vectors
	a = 1	return the n left singular vectors in LSVEC
	a 2	return the first min(N,P) singular vectors in LSVEC
	b = 0	do not compute the right singular vectors
	b = 1	return the right singular vectors in RSVEC
INFO	On exit:
	The singular values and their corresponding singular vectors SVALS(INFO+1), SVALS(INFO+2),...,SVALS(min(N,P)) are correct. If INFO = 0 then all singular values and singular vectors are correct. The matrix B defined as LSVECT × A × RSVEC is the bidiagonal matrix with S on its diagonal and E on its superdiagonal. Therefore the singular values of A and B are the same.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER JOB, LDL, LDR, LDX, N, P

PARAMETER (JOB = 21)

PARAMETER (N = 3)

PARAMETER (P = 3)

PARAMETER (LDL = N)

PARAMETER (LDR = P)

PARAMETER (LDX = N)

C

DOUBLE PRECISION EPSLON, EXCEPT(P), LSVALS(LDL,N)

DOUBLE PRECISION RSVALS(LDR,P), SVALS(N), WORK(N), X(LDX,P)

INTEGER I, ICOL, INFO, IRANK, IROW

C

EXTERNAL DSVDC

INTRINSIC ABS, SQRT

C

C Initialize the array X to store the matrix X shown below.

C

C 1 1 3

C X = 0 1 1

C 1 0 1

C

DATA X / 1.0D0, 0.0D0, 1.0D0, 1.0D0, 1.0D0, 0.0D0,

$ 3.0D0, 1.0D0, 1.0D0 /

C

PRINT 1000

PRINT 1010, ((X(IROW,ICOL), ICOL = 1, N), IROW = 1, N)

CALL DSVDC (X, LDX, N, P, SVALS, EXCEPT, LSVALS, LDL,

$ RSVALS, LDR, WORK, JOB, INFO)

PRINT 1020

PRINT 1030, SVALS

C

C Compute the unit roundoff

C

EPSLON = 1.0D0

10 IF (DBLE (1.0D0 + EPSLON) .NE. 1.0D0) THEN

EPSLON = EPSLON / 2.0D0

GO TO 10

END IF

C

C Make a conservative estimate of the rank of A.

C

IRANK = 0

EPSLON = SQRT (EPSLON + EPSLON)

DO 20, I = 1, N

IF (ABS(SVALS(I)) .GT. EPSLON) THEN

IRANK = IRANK + 1

END IF

20 CONTINUE

PRINT 1040, IRANK

C

1000 FORMAT (1X, 'X:')

1010 FORMAT (3(3X, F4.1))

1020 FORMAT (/1X, 'Singular values of X:')

1030 FORMAT (3X, F4.1)

1040 FORMAT (/1X, 'The rank of X is ', I1)

C

END

Sample Output

X:
1.0 1.0 3.0
0.0 1.0 1.0
1.0 0.0 1.0
Singular values of X:
3.7
1.0
0.3
The rank of X is 3