SQR - Single precision QR decomposition functions adapted from LINPACK FORTRAN:
sqrdc Compute QR decomposition of a matrix.
sqrsl Use QR decomposition to solve for coordinate transformations,
projections, and least squares solutions.
sqrst Solve under-determined or over-determined least squares problems,
with a user-specified tolerance.
Function Prototypes:
void sqrdc (float **x, int n, int p, float *qraux, int *jpvt,
float *work, int job);
void sqrsl (float **x, int n, int k, float *qraux,
float *y, float *qy, float *qty,
float *b, float *rsd, float *xb, int job, int *info);
void sqrst (float **x, int n, int p, float *y, float tol,
float *b, float *rsd, int *k,
int *jpvt, float *qraux, float *work);
sqrdc:
Input:
x matrix[p][n] to decompose (see notes below)
n number of rows in the matrix x
p number of columns in the matrix x
jpvt array[p] controlling the pivot columns (see notes below)
job =0 for no pivoting;
=1 for pivoting
Output:
x matrix[p][n] decomposed (see notes below)
qraux array[p] containing information required to recover the
orthogonal part of the decomposition
jpvt array[p] with jpvt[k] containing the index of the original
matrix that has been interchanged into the k-th column, if
pivoting is requested.
Workspace:
work array[p] of workspace
sqrsl:
Input:
x matrix[p][n] containing output of sqrdc.
n number of rows in the matrix xk; same as in sqrdc.
k number of columns in the matrix xk; k must not be greater
than MIN(n,p), where p is the same as in sqrdc.
qraux array[p] containing auxiliary output from sqrdc.
y array[n] to be manipulated by sqrsl.
job specifies what is to be computed. job has the decimal
expansion ABCDE, with the following meaning:
if A != 0, compute qy.
if B, C, D, or E != 0, compute qty.
if C != 0, compute b.
if D != 0, compute rsd.
if E != 0, compute xb.
Note that a request to compute b, rsd, or xb automatically
triggers the computation of qty, for which an array must
be provided.
Output:
qy array[n] containing Qy, if its computation has been
requested.
qty array[n] containing Q'y, if its computation has
been requested. Here Q' denotes the transpose of Q.
b array[k] containing solution of the least squares problem:
minimize norm2(y - xk*b),
if its computation has been requested. (Note that if
pivoting was requested in sqrdc, the j-th component of
b will be associated with column jpvt[j] of the original
matrix x that was input into sqrdc.)
rsd array[n] containing the least squares residual y - xk*b,
if its computation has been requested. rsd is also the
orthogonal projection of y onto the orthogonal complement
of the column space of xk.
xb array[n] containing the least squares approximation xk*b,
if its computation has been requested. xb is also the
orthogonal projection of y onto the column space of x.
info =0 unless the computation of b has been requested and R
is exactly singular. In this case, info is the index of
the first zero diagonal element of R and b is left
unaltered.
sqrst:
Input:
x matrix[p][n] of coefficients (x is destroyed by sqrst.)
n number of rows in the matrix x (number of observations)
p number of columns in the matrix x (number of parameters)
y array[n] containing right-hand-side (observations)
tol relative tolerance used to determine the subset of
columns of x included in the solution. If tol is zero,
a full complement of the MIN(n,p) columns is used.
If tol is non-zero, the problem should be scaled so that
all the elements of X have roughly the same absolute
accuracy eps. Then a reasonable value for tol is roughly
eps divided by the magnitude of the largest element.
Output:
x matrix[p][n] containing output of sqrdc
b array[p] containing the solution (parameters); components
corresponding to columns not used are set to zero.
rsd array[n] of residuals y - Xb
k number of columns of x used in the solution
jpvt array[p] containing pivot information from sqrdc.
qraux array[p] containing auxiliary information from sqrdc.
Workspace:
work array[p] of workspace.
Notes:
!!! WARNING !!!
These functions have many options, not all of which have been tested!
(Dave Hale, 12/31/89)
This function was adapted from LINPACK FORTRAN. Because two-dimensional
arrays cannot be declared with variable dimensions in C, the matrix x
is actually a pointer to an array of pointers to floats, as declared
above and used below.
Elements of x are stored as follows:
x[0][0] x[1][0] x[2][0] ... x[p-1][0]
x[0][1] x[1][1] x[2][1] ... x[p-1][1]
x[0][2] x[1][2] x[2][2] ... x[p-1][2]
. .
. . .
. . .
. .
x[0][n-1] x[1][n-1] x[2][n-1] ... x[p-1][n-1]
sqrdc:
Uses Householder transformations to compute the QR decomposition of an n by p
matrix x. Column pivoting based on the 2-norms of the reduced columns may be
performed at the user's option.
After decomposition, x contains in its upper triangular matrix R of the QR
decomposition. Below its diagonal x contains information from which the
orthogonal part of the decomposition can be recovered. Note that if
pivoting has been requested, the decomposition is not that of the original
matrix x but that of x with its columns permuted as described by jpvt.
The selection of pivot columns is controlled by jpvt as follows.
The k-th column x[k] of x is placed in one of three classes according to
the value of jpvt[k].
if jpvt[k] > 0, then x[k] is an initial column.
if jpvt[k] == 0, then x[k] is a free column.
if jpvt[k] < 0, then x[k] is a final column.
Before the decomposition is computed, initial columns are moved to the
beginning of the array x and final columns to the end. Both initial and
final columns are frozen in place during the computation and only free
columns are moved. At the k-th stage of the reduction, if x[k] is occupied
by a free column it is interchanged with the free column of largest reduced
norm. jpvt is not referenced if job == 0.
sqrsl:
Uses the output of sqrdc to compute coordinate transformations, projections,
and least squares solutions. For k <= MIN(n,p), let xk be the matrix
xk = (x[jpvt[0]], x[jpvt[1]], ..., x[jpvt[k-1]])
formed from columns jpvt[0], jpvt[1], ..., jpvt[k-1] of the original
n by p matrix x that was input to sqrdc. (If no pivoting was done, xk
consists of the first k columns of x in their original order.) sqrdc
produces a factored orthogonal matrix Q and an upper triangular matrix R
such that
xk = Q * (R)
(0)
This information is contained in coded form in the arrays x and qraux.
The parameters qy, qty, b, rsd, and xb are not referenced if their
computation is not requested and in this case can be replaced by NULL
pointers in the calling program. To save storage, the user may in
some cases use the same array for different parameters in the calling
sequence. A frequently occuring example is when one wishes to compute
any of b, rsd, or xb and does not need y or qty. In this case one may
equivalence y, qty, and one of b, rsd, or xb, while providing separate
arrays for anything else that is to be computed. Thus the calling
sequence
sqrsl(x,n,k,qraux,y,NULL,y,b,y,NULL,110,&info)
will result in the computation of b and rsd, with rsd overwriting y.
More generally, each item in the following list contains groups of
permissible equivalences for a single calling sequence.
1. (y,qty,b) (rsd) (xb) (qy)
2. (y,qty,rsd) (b) (xb) (qy)
3. (y,qty,xb) (b) (rsd) (qy)
4. (y,qy) (qty,b) (rsd) (xb)
5. (y,qy) (qty,rsd) (b) (xb)
6. (y,qy) (qty,xb) (b) (rsd)
In any group the value returned in the array allocated to the group
corresponds to the last member of the group.
sqrst:
Computes least squares solutions to the system
Xb = y
which may be either under-determined or over-determined. The user may
supply a tolerance to limit the columns of X used in computing the solution.
In effect, a set of columns with a condition number approximately bounded
by 1/tol is used, the other components of b being set to zero.
On return, the arrays x, jpvt, and qraux contain the usual output from
sqrdc, so that the QR decomposition of x with pivoting is fully available
to the user. In particular, columns jpvt[0], jpvt[1], ..., jpvt[k-1]
were used in the solution, and the condition number associated with
those columns is estimated by ABS(x[0][0]/x[k-1][k-1]).
Author: Dave Hale, Colorado School of Mines, 12/29/89