sqr.c

seismic software,very useful

		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.
******************************************************************************
Notes:
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[k-1][0]
x[0][1]    x[1][1]    x[2][1]   ... x[k-1][1]
x[0][2]    x[1][2]    x[2][2]   ... x[k-1][2]
.                                       .
.             .                         .
.                        .              .
.                                       .
x[0][n-1]  x[1][n-1]  x[2][n-1] ... x[k-1][n-1]

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.
******************************************************************************
Author:		Dave Hale, Colorado School of Mines, 12/29/89
*****************************************************************************/
{
	int i,j,ju,cb,cqy,cqty,cr,cxb;
	float t,temp;
	
	/* set info flag */
	*info = 0;
	
	/* determine what is to be computed */
	cqy = job/10000!=0;
	cqty = job%10000!=0;
	cb = (job%1000)/100!=0;
	cr = (job%100)/10!=0;
	cxb = job%10!=0;
	ju = MIN(k,n-1);
	
	/* special action when n=1 */
	if (ju==0) {
		if (cqy) qy[0] = y[0];
		if (cqty) qty[0] = y[0];
		if (cxb) xb[0] = y[0];
		if (cb) {
			if (x[0][0]==0.0)
				*info = 1;
			else
				b[0] = y[0]/x[0][0];
		}
		if (cr) rsd[0] = 0.0;
		return;
	}
	
	/* set up to compute Qy or Q'y */
	if (cqy) scopy(n,y,1,qy,1);
	if (cqty) scopy(n,y,1,qty,1);
	if (cqy) {
	
		/* compute Qy */
		for (j=ju-1; j>=0; j--) {
			if (qraux[j]!=0.0) {
				temp = x[j][j];
				x[j][j] = qraux[j];
				t = -sdot(n-j,&x[j][j],1,&qy[j],1)/x[j][j];
				saxpy(n-j,t,&x[j][j],1,&qy[j],1);
				x[j][j] = temp;
			}
		}
	}
	if (cqty) {
		
		/* compute Q'y */
		for (j=0; j<ju; j++) {
			if (qraux[j]!=0.0) {
				temp = x[j][j];
				x[j][j] = qraux[j];
				t = -sdot(n-j,&x[j][j],1,&qty[j],1)/x[j][j];
				saxpy(n-j,t,&x[j][j],1,&qty[j],1);
				x[j][j] = temp;
			}
		}
	}
	
	/* set up to compute b, rsd, or xb */
	if (cb) scopy(k,qty,1,b,1);
	if (cxb) scopy(k,qty,1,xb,1);
	if (cr && k<n) scopy(n-k,&qty[k],1,&rsd[k],1);
	if (cxb && k<n)
		for (i=k; i<n; i++)
			xb[i] = 0.0;
	if (cr)
		for (i=0; i<k; i++)
			rsd[i] = 0.0;
	if (cb) {
	
		/* compute b */
		for (j=k-1; j>=0; j--) {
			if (x[j][j]==0.0) {
				*info = j;
				break;
			}
			b[j] /= x[j][j];
			if (j!=0) {
				t = -b[j];
				saxpy(j,t,x[j],1,b,1);
			}
		}
	}
	if (cr || cxb) {
	
		/* compute rsd or xb as requested */
		for (j=ju-1; j>=0; j--) {
			if (qraux[j]!=0.0) {
				temp = x[j][j];
				x[j][j] = qraux[j];
				if (cr) {
					t = -sdot(n-j,&x[j][j],1,&rsd[j],1)/
						x[j][j];
					saxpy(n-j,t,&x[j][j],1,&rsd[j],1);
				}
				if (cxb) {
					t = -sdot(n-j,&x[j][j],1,&xb[j],1)/
						x[j][j];
					saxpy(n-j,t,&x[j][j],1,&xb[j],1);
				}
				x[j][j] = temp;
			}
		}
	}
}
							
void
sqrst (float **x, int n, int p, float *y, float tol,
	float *b, float *rsd, int *k,
	int *jpvt, float *qraux, float *work)
/*****************************************************************************
Compute 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.
******************************************************************************
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:
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[k-1][0]
x[0][1]    x[1][1]    x[2][1]   ... x[k-1][1]
x[0][2]    x[1][2]    x[2][2]   ... x[k-1][2]
.                                       .
.             .                         .
.                        .              .
.                                       .
x[0][n-1]  x[1][n-1]  x[2][n-1] ... x[k-1][n-1]

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/31/89
*****************************************************************************/
{
	int info,j,kk,m;
	float t;
	
	/* initialize jpvt so that all columns are free */
	for (j=0; j<p; j++)
		jpvt[j] = 0;
	
	/* decompose x */
	sqrdc(x,n,p,qraux,jpvt,work,1);
	
	/* determine which columns to use */
	m = MIN(n,p);
	for (j=0,kk=0; j<m && ABS(x[j][j])>tol*ABS(x[0][0]); j++)
		kk++;
	
	/* solve the truncated least squares problem */
	if (kk>0) sqrsl(x,n,kk,qraux,y,rsd,rsd,work,rsd,rsd,110,&info);
	
	/* unscramble the components of b corresponding to columns used */
	for (j=0; j<kk; j++)
		b[jpvt[j]] = work[j];
	
	/* zero components of b corresponding to unused columns */
	for (j=kk; j<p; j++)
		b[jpvt[j]] = 0.0;

	/* return number of columns used */
	*k = kk;
}