Singular Value Decomposition of a General Matrix

Calling Sequence

CALL DGESVD	(JOBU, JOBVT, M, N, DA, LDA, DSING, DU, LDU, DVT, LDVT, DWORK, LDWORK, INFO)
CALL SGESVD	(JOBU, JOBVT, M, N, SA, LDA, SSING, SU, LDU, SVT, LDVT, SWORK, LDWORK, INFO)
CALL ZGESVD	(JOBU, JOBVT, M, N, ZA, LDA, DSING, ZU, LDU, ZVT, LDVT, ZWORK, LDWORK, DWORK2, INFO)
CALL CGESVD	(JOBU, JOBVT, M, N, CA, LDA, SSING, CU, LDU, CVT, LDVT, CWORK, LDWORK, SWORK2, INFO)
void dgesvd	(char jobu, char jobvt, int m, int n, double *da, int lda, double *dsing, double *du, int ldu, double *dvt, int ldvt, int *info)
void sgesvd	(char jobu, char jobvt, int m, int n, float *sa, int lda, float *ssing, float *su, int ldu, float *svt, int ldvt, int *info)
void zgesvd	(char jobu, char jobvt, int m, int n, doublecomplex *za, int lda, double *dsing, doublecomplex *u, int ldu, doublecomplex *zvt, int ldvt, int *info)
void cgesvd	(char jobu, char jobvt, int m, int n, complex *ca, int lda, float *ssing, complex *u, int ldu, complex *cvt, int ldvt, int *info)

Arguments

JOBU	Indicates which columns of U, which are the left singular vectors will be computed and returned. Legal values of JOBU are shown below. Any values not shown below are illegal.
	'A' or 'a'	All M left singular vectors are returned in U.
	'N' or 'n'	No left singular vectors are computed.
	'O' or 'o'	The first min(M, N) left singular vectors are overwritten on A.
	'S' or 's'	Only the first min(M, N) left singular vectors are returned in U.
	Note that JOBU and JOBVT cannot both direct that their respective results be written to A.
JOBVT	Indicates which rows of VT, which are the right singular vectors, will be computed and returned. Legal values of JOBVT are shown below. Any values not shown below are illegal.
	'A' or 'a'	All N right singular vectors are returned in VT.
	'N' or 'n'	No right singular vectors are computed.
	'O' or 'o'	The first min(M, N) right singular vectors are overwritten on A.
	'S' or 's'	Only the first min(M, N) right singular vectors are returned in VT.
	Note that JOBU and JOBVT cannot both direct that their respective results be written to A.
M	Number of rows of the input matrix A. M 0.
N	Number of columns of the input matrix A. N 0.
xA	On entry, the matrix A. On exit, if JOBU = 'O' or 'o' then A is overwritten with the first min(M, N) columns of the array U, which contain the left singular vectors stored columnwise. On exit, if JOBVT = 'O' or 'o' then A is overwritten with the first min(M, N) rows of VT, which contain the right singular vectors stored rowwise. On exit, if neither JOBU nor JOBVT = 'O' or 'o' then the contents of A are overwritten.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA max(1, M).
xSING	The min(M,N) singular values of A, in descending order. The singular values are the diagonal elements of the diagonal matrix .
xU	On exit, the contents of U are a function of the entry value of JOBU:
	JOBU	Size of U	Contents of U
	'A' or 'a'	M×M	Orthogonal or unitary matrix U
	'N' or 'n'	Not used	U is not used
	'O' or 'o'	Not used	U is not used
	'S' or 's'	M×min(M, N)	First min(M, N) columns of the orthogonal or unitary matrix U
LDU	Leading dimension of the array U as specified in a dimension or type statement. LDU 1.
	If JOBU = 'A', 'a', 'S', or 's' then LDU max(1, M).
xVT	On exit, the contents of VT are a function of the entry value of JOBVT:
	JOBVT	Size of VT	Contents of VT
	'A' or 'a'	N×N	Orthogonal or unitary matrix VT
	'N' or 'n'	Not used	VT is not used
	'O' or 'o'	Not used	VT is not used
	'S' or 's'	M×min(M, N)	First min(M, N) rows of the orthogonal or unitary matrix VT stored rowwise
	Note that VT is returned and not V.
LDVT	Leading dimension of the array VT as specified in a dimension or type statement.
	If JOBVT = 'A' or 'a' then LDVT max(1, N).
	If JOBVT = 'S' or 's' then LDVT min(M, N) and LDVT 1.
	For other values of JOBVT, LDVT 1.
xWORK	Scratch array with a dimension of LDWORK.
LDWORK	Leading dimension of the array WORK as specified in a dimension or type statement. LDWORK is at least 1 but also greater than or equal to:
	for real subroutines: the larger of (3 × min(M, N) + max(M, N)) and (5 × min(M, N) - 4).
	for complex subroutines: (2 × min(M, N) + max(M, N)).
xWORK2	Scratch array with a dimension of 5 × max(M, N) for complex subroutines.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO < 0	The ith argument, where i = |INFO|, had an illegal value.
	INFO > 0	Convergence failure.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDA, LDWORK, M, N

PARAMETER (M = 4)

PARAMETER (N = 4)

PARAMETER (LDA = M)

PARAMETER (LDWORK = 5 * M)

C

DOUBLE PRECISION A(LDA,N), TEMP, SING(M), WORK(LDWORK)

INTEGER ICOL, INFO, IROW

C

EXTERNAL DGESVD

INTRINSIC ABS

C

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

C

C 1/2 1/2 1/2 1/2

C A = 1/2 1/2 -1/2 -1/2

C 1/2 -1/2 1/2 -1/2

C 1/2 -1/2 -1/2 1/2

C

DATA A / 5.0D-1, 5.0D-1, 5.0D-1, 5.0D-1,

$ 5.0D-1, 5.0D-1, -5.0D-1, -5.0D-1,

$ 5.0D-1, -5.0D-1, 5.0D-1, -5.0D-1,

$ 5.0D-1, -5.0D-1, -5.0D-1, 5.0D-1 /

C

C Print the initial value of A.

C

PRINT 1000

PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, M)

C

C Compute and print the singular values of A.

C

CALL DGESVD ('NO LEFT SINGULAR VECTORS',

$ 'NO RIGHT SINGULAR VECTORS', M, N, A, LDA,

$ SING, TEMP, 1, TEMP, 1, WORK, LDWORK, INFO)

IF (INFO .NE. 0) THEN

IF (INFO .LT. 0) THEN

PRINT 1020, ABS (INFO)

STOP 1

ELSE

PRINT 1030

STOP 2

END IF

END IF

PRINT 1040

PRINT 1050, SING

C

1000 FORMAT (1X, 'A:')

1010 FORMAT (4(3X, F8.4))

1020 FORMAT (1X, 'Illegal value for argument #', I2,

$ ' in DGESVD.')

1030 FORMAT (1X, 'Singular value decomposition failed',

$ 1X, 'to converge.')

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

1050 FORMAT (1X, F8.4)

C

END

Sample Output

A:
0.5000 0.5000 0.5000 0.5000
0.5000 0.5000 -0.5000 -0.5000
0.5000 -0.5000 0.5000 -0.5000
0.5000 -0.5000 -0.5000 0.5000
Singular values of A:
1.0000
1.0000
1.0000
1.0000