Purpose
To compute a QR factorization of an n-by-m matrix A (A = Q * R), having a p-by-min(p,m) zero triangle in the lower left-hand side corner, as shown below, for n = 8, m = 7, and p = 2: [ x x x x x x x ] [ x x x x x x x ] [ x x x x x x x ] [ x x x x x x x ] A = [ x x x x x x x ], [ x x x x x x x ] [ 0 x x x x x x ] [ 0 0 x x x x x ] and optionally apply the transformations to an n-by-l matrix B (from the left). The problem structure is exploited. This computation is useful, for instance, in combined measurement and time update of one iteration of the time-invariant Kalman filter (square root information filter).Specification
SUBROUTINE MB04ID( N, M, P, L, A, LDA, B, LDB, TAU, DWORK, LDWORK, $ INFO ) C .. Scalar Arguments .. INTEGER INFO, L, LDA, LDB, LDWORK, M, N, P C .. Array Arguments .. DOUBLE PRECISION A(LDA,*), B(LDB,*), DWORK(*), TAU(*)Arguments
Input/Output Parameters
N (input) INTEGER The number of rows of the matrix A. N >= 0. M (input) INTEGER The number of columns of the matrix A. M >= 0. P (input) INTEGER The order of the zero triagle. P >= 0. L (input) INTEGER The number of columns of the matrix B. L >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,M) On entry, the leading N-by-M part of this array must contain the matrix A. The elements corresponding to the zero P-by-MIN(P,M) lower trapezoidal/triangular part (if P > 0) are not referenced. On exit, the elements on and above the diagonal of this array contain the MIN(N,M)-by-M upper trapezoidal matrix R (R is upper triangular, if N >= M) of the QR factorization, and the relevant elements below the diagonal contain the trailing components (the vectors v, see Method) of the elementary reflectors used in the factorization. LDA INTEGER The leading dimension of array A. LDA >= MAX(1,N). B (input/output) DOUBLE PRECISION array, dimension (LDB,L) On entry, the leading N-by-L part of this array must contain the matrix B. On exit, the leading N-by-L part of this array contains the updated matrix B. If L = 0, this array is not referenced. LDB INTEGER The leading dimension of array B. LDB >= MAX(1,N) if L > 0; LDB >= 1 if L = 0. TAU (output) DOUBLE PRECISION array, dimension MIN(N,M) The scalar factors of the elementary reflectors used.Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK) On exit, if INFO = 0, DWORK(1) returns the optimal value of LDWORK. LDWORK The length of the array DWORK. LDWORK >= MAX(1,M-1,M-P,L). For optimum performance LDWORK should be larger. If LDWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the DWORK array, returns this value as the first entry of the DWORK array, and no error message related to LDWORK is issued by XERBLA.Error Indicator
INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value.Method
The routine uses min(N,M) Householder transformations exploiting the zero pattern of the matrix. A Householder matrix has the form ( 1 ), H = I - tau *u *u', u = ( v ) i i i i i ( i) where v is an (N-P+I-2)-vector. The components of v are stored i i in the i-th column of A, beginning from the location i+1, and tau is stored in TAU(i). iNumerical Aspects
The algorithm is backward stable.Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None