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

Calling Sequence

CALL DGELSX	(M, N, NRHS, DA, LDA, DB, LDB, JPIVOT, DRCOND, IRANK, DWORK, INFO)
CALL SGELSX	(M, N, NRHS, SA, LDA, SB, LDB, JPIVOT, SRCOND, IRANK, SWORK, INFO)
CALL ZGELSX	(M, N, NRHS, ZA, LDA, ZB, LDB, JPIVOT, DRCOND, IRANK, ZWORK, DWORK2, INFO)
CALL CGELSX	(M, N, NRHS, CA, LDA, CB, LDB, JPIVOT, SRCOND, IRANK, CWORK, SWORK2, INFO)
void dgelsx	(int m, int n, int nrhs, double *da, int lda, double *db, int ldb, int *jpivot, double rcond, int *rank, int *info)
void sgelsx	(int m, int n, int nrhs, float *sa, int lda, float *sb, int ldb, int *jpivot, float rcond, int *rank, int *info)
void zgelsx	(int m, int n, int nrhs, doublecomplex *za, int lda, doublecomplex *zb, int ldb, int *jpivot, double rcond, int *rank, int *info)
void cgelsx	(int m, int n, int nrhs, complex *ca, int lda, complex *cb, int ldb, int *jpivot, 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 matrices B and X. NRHS 0.
xA	On entry, the matrix A. On exit, A has been overwritten by details of its complete orthogonal factorization.
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.
LDB	Leading dimension of the array B as specified in a dimension or type statement. LDB max(1, M, N).
JPIVOT	On entry, if JPIVOT(i) 0 for 1 i N then the ith column of A is an initial column, otherwise it is a free column. Before the QR factorization of A, all initial columns are permuted to the leading positions. Only the remaining free columns are moved as a result of column pivoting during the factorization. On exit, if JPIVOT(i) k, then the ith column of A was the kth column of A.
xRCOND	Determines the effective rank of A, defined as the order of the largest leading triangular submatrix R11 in the QR factorization with pivoting of A, whose estimated condition number < 1 / RCOND.
IRANK	On exit, the effective rank of A, equal to the order of the submatrix R11 to which the description of RCOND refers.
xWORK	Scratch array whose dimension is the larger of:
	for real subroutines: (min(M, N) + 3 × N) and (2 × min(M, N) + NRHS).
	for complex subroutines: (min(M, N) + N) and (min(M, N) + 2 × min(M, N) + NRHS).
xWORK2	Scratch array with a dimension of 2 × N for complex subroutines.
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO < 0	The ith argument, where i = |INFO|, had an illegal value.

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 = 4 * N)

PARAMETER (NRHS = 1)

C

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

DOUBLE PRECISION WORK(LDWORK)

INTEGER ICOL, INFO, IPIVOT(N), IROW, IRANK

C

EXTERNAL DGELSX

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 /

DATA IPIVOT / N*1 /

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 Solve the minimization problem. Print the solution and the

C numerical rank of A.

C

RCOND = SQRT(2.2D-16)

CALL DGELSX (M, N, NRHS, A, LDA, B, LDB, IPIVOT, RCOND,

$ IRANK, WORK, 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 #', I2,

$ ' in DGELSX.')

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