Solution to a Least Squares Problem for a Rank-Deficient General Matrix (Simple Driver)

Calling Sequence

CALL DGELSS	(M, N, NRHS, DA, LDA, DB, LDB, DSING, DRCOND, IRANK, DWORK, LDWORK, INFO)
CALL SGELSS	(M, N, NRHS, SA, LDA, SB, LDB, SSING, SRCOND, IRANK, SWORK, LDWORK, INFO)
CALL ZGELSS	(M, N, NRHS, ZA, LDA, ZB, LDB, DSING, DRCOND, IRANK, ZWORK, LDWORK, DWORK2, INFO)
CALL CGELSS	(M, N, NRHS, CA, LDA, CB, LDB, SSING, SRCOND, IRANK, CWORK, LDWORK, SWORK2, INFO)
void dgelss	(int m, int n, int nrhs, double *da, int lda, double *db, int ldb, double *dsing, double rcond, int *rank, int *info)
void sgelss	(int m, int n, int nrhs, float *sa, int lda, float *sb, int ldb, float *ssing, float rcond, int *rank, int *info)
void zgelss	(int m, int n, int nrhs, doublecomplex *za, int lda, doublecomplex *zb, int ldb, double *dsing, double rcond, int *rank, int *info)
void cgelss	(int m, int n, int nrhs, complex *ca, int lda, complex *cb, int ldb, float *ssing, float rcond, int *rank, int *info)

Arguments

M	Number of rows of the matrix A. M 0.
N	Number of columns of the matrix A. N 0.
NRHS	Number of right-hand sides, equal to the number of columns of the matrices B and X. NRHS 0.
xA	On entry, the matrix A. On exit, the first min(M, N) rows of A are overwritten with its right singular vectors, stored rowwise.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA max(1, M).
xB	On entry, the M×NRHS right-hand side matrix B. On exit, the N×NRHS solution matrix X. If M N and IRANK = N then the residual sum-of-squares for the solution in the ith column is given by the sum of squares of elements N+1 through M in that column. The effective rank of A is determined by treating as zero those singular values that are less than RCOND times the largest singular value.
LDB	Leading dimension of the array B as specified in a dimension or type statement. LDB max(1, M, N).
xSING	On exit, the min(M, N) singular values of A in decreasing order. The condition number of A in the 2-norm = SING(1) / SING(min(M, N)). Of course, one must be cautious in computing the condition number because of possible division by zero, overflow, and other potential difficulties.
xRCOND	Determines the effective rank of A. Singular values that satisfy SING(i) RCOND × SING(1) are treated as zero. If RCOND < 0 then machine precision is used instead.
IRANK	On exit, the effective rank of A, equal to the number of singular values which are greater than RCOND × SING(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 1 and must also satisfy:
	for M N	LDWORK 3 × N + max(2 × N, NRHS, M) for real subroutines LDWORK 2 × N + max(NRHS, M) for complex subroutines
	for M < N	LDWORK 3 × M + max(2 × M, NRHS, N) for real subroutines LDWORK 2 × M + max(NRHS, N) for complex subroutines
xWORK2	Scratch array with dimension that is the larger of 1 and (5 × min(M, N) - 4) 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	The algorithm for computing the SVD failed to converge.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER LDA, LDB, LDSING, LDWORK, M, N, NRHS

PARAMETER (M = 4)

PARAMETER (N = 4)

PARAMETER (LDA = M)

PARAMETER (LDB = M)

PARAMETER (LDSING = M)

PARAMETER (LDWORK = 5 * N)

PARAMETER (NRHS = 1)

C

DOUBLE PRECISION A(LDA,N), B(LDB,NRHS), RCOND, SING(LDSING)

DOUBLE PRECISION SMALL, WORK(LDWORK)

INTEGER ICOL, INFO, IROW, IRANK

C

EXTERNAL DGELSS

INTRINSIC ABS, SQRT

C

C Initialize the array A to store the coefficient matrix A

C shown below. Initialize the array B to store the right

C hand side vector b shown below. The element labeled z in

C the last column of A differs from 2 only in the last few

C bits. This makes columns three and four nearly equal.

C

C 1 6 2 z 96

C A = 1 -2 -8 -8 b = 192

C 1 -2 4 4 192

C 1 6 14 14 -96

C

DATA A / 1.0D0, 1.0D0, 1.0D0, 1.0D0,

$ 6.0D0, -2.0D0, -2.0D0, 6.0D0,

$ 2.0D0, -8.0D0, 4.0D0, 1.4D1,

$ 2.0D0, -8.0D0, 4.0D0, 1.4D1 /

DATA B / 9.6D1, 1.92D2, 1.92D2, -9.6D1 /

C

C Print the initial value of the arrays.

C

SMALL = ((((2.0D0 / 3.0D0) + 4.0D0) - 4.0D0) -

$ (2.0D0 / 3.0D0))

A(1,N) = A(1,N-1) + SMALL

PRINT 1000

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

PRINT 1020

PRINT 1030, B

C

C Compute and print the rank of A and the solution to the

C minimization problem.

C

RCOND = SQRT(2.2D-16)

CALL DGELSS (M, N, NRHS, A, LDA, B, LDB, SING, RCOND, IRANK,

$ WORK, LDWORK, INFO)

IF (INFO .NE. 0) THEN

PRINT 1040, ABS(INFO)

STOP 1

END IF

PRINT 1050, IRANK

PRINT 1060

PRINT 1030, B

C

1000 FORMAT (1X, 'A:')

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

1020 FORMAT (/1X, 'b:')

1030 FORMAT (1X, F8.4)

1040 FORMAT (1X, 'Illegal value for argument #', I1,

$ ' in DGELSS.')

1050 FORMAT (/1X, 'Numerical rank of A: ', I1)

1060 FORMAT (/1X, 'min(b - Ax):')

C

END

Sample Output

A:
1.0000 6.0000 2.0000 2.0000
1.0000 -2.0000 -8.0000 -8.0000
1.0000 -2.0000 4.0000 4.0000
1.0000 6.0000 14.0000 14.0000
b:
96.0000
192.0000
192.0000
-96.0000
Numerical rank of A: 3
min(b - Ax):
148.0000
-14.0000
-4.0000
-4.0000