Purpose
To determine the minimum-norm solution to a real linear least squares problem: minimize || A * X - B ||, using the rank-revealing QR factorization of a real general M-by-N matrix A, computed by SLICOT Library routine MB03OD.Specification
SUBROUTINE MB02QY( M, N, NRHS, RANK, A, LDA, JPVT, B, LDB, TAU, $ DWORK, LDWORK, INFO ) C .. Scalar Arguments .. INTEGER INFO, LDA, LDB, LDWORK, M, N, NRHS, RANK C .. Array Arguments .. INTEGER JPVT( * ) DOUBLE PRECISION A( LDA, * ), B( LDB, * ), DWORK( * ), TAU( * )Arguments
Input/Output Parameters
M (input) INTEGER The number of rows of the matrices A and B. M >= 0. N (input) INTEGER The number of columns of the matrix A. N >= 0. NRHS (input) INTEGER The number of columns of the matrix B. NRHS >= 0. RANK (input) INTEGER The effective rank of A, as returned by SLICOT Library routine MB03OD. min(M,N) >= RANK >= 0. A (input/output) DOUBLE PRECISION array, dimension ( LDA, N ) On entry, the leading min(M,N)-by-N upper trapezoidal part of this array contains the triangular factor R, as returned by SLICOT Library routine MB03OD. The strict lower trapezoidal part of A is not referenced. On exit, if RANK < N, the leading RANK-by-RANK upper triangular part of this array contains the upper triangular matrix R of the complete orthogonal factorization of A, and the submatrix (1:RANK,RANK+1:N) of this array, with the array TAU, represent the orthogonal matrix Z (of the complete orthogonal factorization of A), as a product of RANK elementary reflectors. On exit, if RANK = N, this array is unchanged. LDA INTEGER The leading dimension of the array A. LDA >= max(1,M). JPVT (input) INTEGER array, dimension ( N ) The recorded permutations performed by SLICOT Library routine MB03OD; if JPVT(i) = k, then the i-th column of A*P was the k-th column of the original matrix A. B (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS ) On entry, if NRHS > 0, the leading M-by-NRHS part of this array must contain the matrix B (corresponding to the transformed matrix A, returned by SLICOT Library routine MB03OD). On exit, if NRHS > 0, the leading N-by-NRHS part of this array contains the solution matrix X. If M >= N and RANK = N, the residual sum-of-squares for the solution in the i-th column is given by the sum of squares of elements N+1:M in that column. If NRHS = 0, the array B is not referenced. LDB INTEGER The leading dimension of the array B. LDB >= max(1,M,N), if NRHS > 0. LDB >= 1, if NRHS = 0. TAU (output) DOUBLE PRECISION array, dimension ( min(M,N) ) The scalar factors of the elementary reflectors. If RANK = N, the array TAU is not referenced.Workspace
DWORK DOUBLE PRECISION array, dimension ( LDWORK ) On exit, if INFO = 0, DWORK(1) returns the optimal value of LDWORK. LDWORK INTEGER The length of the array DWORK. LDWORK >= max( 1, N, NRHS ). For good performance, LDWORK should sometimes be larger.Error Indicator
INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value.Method
The routine uses a QR factorization with column pivoting: A * P = Q * R = Q * [ R11 R12 ], [ 0 R22 ] where R11 is an upper triangular submatrix of estimated rank RANK, the effective rank of A. The submatrix R22 can be considered as negligible. If RANK < N, then R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization: A * P = Q * [ T11 0 ] * Z. [ 0 0 ] The minimum-norm solution is then X = P * Z' [ inv(T11)*Q1'*B ], [ 0 ] where Q1 consists of the first RANK columns of Q. The input data for MB02QY are the transformed matrices Q' * A (returned by SLICOT Library routine MB03OD) and Q' * B. Matrix Q is not needed.Numerical Aspects
The implemented method is numerically stable.Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None