Solution to a Linear System in a QR-Factored General Matrix

Calling Sequence

CALL DQRSL	(DA, LDA, N, K, DQRAUX, DY, DQY, DQTY, DB, DRESID, DAB, JOB, INFO)
CALL SQRSL	(SA, LDA, N, K, SQRAUX, SY, SQY, SQTY, SB, SRESID, SAB, JOB, INFO)
CALL ZQRSL	(ZA, LDA, N, K, ZQRAUX, CY, CQY, ZQTY, ZB, ZRESID, ZAB, JOB, INFO)
CALL CQRSL	(CA, LDA, N, K, CQRAUX, ZY, ZQY, CQTY, CB, CRESID, CAB, JOB, INFO)
void dqrsl	(double *da, int lda, int n, int k, double *dqraux, double *dy, double *dqy, double *dqty, double *db, double *dresid, double *dab, int job, int *info)
void sqrsl	(float *sa, int lda, int n, int k, float *sqraux, float *sy, float *sqy, float *sqty, float *sb, float *sresid, float *sab, int job, int *info)
void zqrsl	(doublecomplex *za, int lda, int n, int k, doublecomplex *zqraux,doublecomplex *zy, doublecomplex *zqy, doublecomplex *zqty, doublecomplex *zb, doublecomplex *zresid, doublecomplex *zab, int job, int *info)
void cqrsl	(complex *ca, int lda, int n, int k, complex *cqraux, complex *cy, complex *cqy, complex *cqty, complex *cb, complex *cresid, complex *cab, int job, int *info)

Arguments

xA	Part of the QR factorization of matrix A as computed by xSRDC.
LDA	Leading dimension of the array A as specified in a dimension or type statement. LDA max(1,N).
N	Number of rows in the matrix AK where AK is described below. N 0.
K	Number of columns in the matrix AK where AK is described below. K 0.
xQRAUX	Auxiliary output from xQRDC.
xY	Vector to be manipulated by xQRSL.
xQY	On exit, QY contains Q × Y if its computation has been requested in JOB; QY is not referenced if its computation is not requested.
xQTY	On exit, QTY contains QT × Y if its computation has been requested in JOB; QTY is not referenced if its computation is not requested.
xB	On entry, the right-hand side vector b. On exit, the solution vector x. B is not referenced if its computation is not requested.
xRESID	On exit, RESID contains the least squares residual y - AK × b if its computation has been requested. RESID also is the orthogonal projection of y onto the orthogonal complement of the column space of AK. RESID is not referenced if its computation is not requested.
xAB	On exit, AB contains the least squares approximation AK × b if its computation has been requested. AB is the orthogonal projection of y onto the column space of x. AB is not referenced if its computation is not requested.
JOB	Integer in the form abcde; determines which operation or operations the subroutine will perform:
	a 0	compute QY
	b, c, d, or e 0	compute QTY
	c 0	compute B
	d 0	compute RESID
	e 0	compute AB
INFO	On exit:
	INFO = 0	Subroutine completed normally.
	INFO > 0	Returns a value k if the computation of B has been requested and R is singular;
		the value of k is then the index of the first zero element of R.

Sample Program

PROGRAM TEST

IMPLICIT NONE

C

INTEGER IDOB, IDORSD, IDOXB, LDA, N, NCOLA, NOPIV

INTEGER NROWA

PARAMETER (IDOB = 100)

PARAMETER (IDORSD = 10)

PARAMETER (IDOXB = 1)

PARAMETER (N = 3)

PARAMETER (LDA = N)

PARAMETER (NCOLA = 2)

PARAMETER (NOPIV = 0)

PARAMETER (NROWA = N)

C

DOUBLE PRECISION A(LDA,NCOLA), B(NCOLA), NULL(N), QRAUX(N)

DOUBLE PRECISION RESID(N), WORK(N), Y(N)

INTEGER ICOL, INFO, IROW, JOB, JPIVOT

C

EXTERNAL DQRDC, DQRSL

C

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

C Initialize the array Y to store the vector y shown below.

C

C 1 1 1

C A = 1 0 y = 0

C 0 1 -5

C

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

DATA Y / 1.0D0, 0.0D0, -5.0D0 /

C

PRINT 1000

PRINT 1010, ((A(IROW,ICOL), ICOL = 1, NCOLA), IROW = 1,

$ NROWA)

PRINT 1020

PRINT 1030, Y

JOB = NOPIV

CALL DQRDC (A, LDA, NROWA, NCOLA, QRAUX, JPIVOT, WORK, JOB)

JOB = IDOB + IDORSD + IDOXB

CALL DQRSL (A, LDA, NROWA, NCOLA, QRAUX, Y, NULL, NULL, B,

$ RESID, NULL, JOB, INFO)

IF (INFO .EQ. 0) THEN

PRINT 1040

PRINT 1050, B

PRINT 1060

PRINT 1050, RESID

ELSE

PRINT 1070

END IF

C

1000 FORMAT (1X, 'A:')

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

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

1030 FORMAT (3X, F4.1)

1040 FORMAT (/1X, 'Least squares solution:')

1050 FORMAT (3X, F4.1)

1060 FORMAT (/1X, 'Residual:')

1070 FORMAT (1X, 'A is singular.')

C

END

Sample Output

A:
1.0 1.0
1.0 0.0
0.0 1.0
y:
1.0
0.0
-5.0
Least squares solution:
2.0
-3.0
Residual:
2.0
-2.0
-2.0