sqr.c

该程序是用vc开发的对动态数组进行管理的DLL

			}
		}
		qraux[l] = 0.0;
		if (l!=n-1) {
		
			/*
			 * compute the Householder transformation
			 * for column l
			 */
			nrmxl = snrm2(n-l,&x[l][l],1);
			if (nrmxl!=0.0) {
				if (x[l][l]!=0.0)
					nrmxl = (x[l][l]>0.0) ?
						ABS(nrmxl) :
						-ABS(nrmxl);
				sscal(n-l,1.0f/nrmxl,&x[l][l],1);
				x[l][l] += 1.0f;
				
				/*
				 * apply the transformation to the remaining
				 * columns, updating the norms
				 */
				 for (j=l+1; j<p; j++) {
					 t = -sdot(n-l,&x[l][l],1,&x[j][l],1)/
					 	x[l][l];
					saxpy(n-l,t,&x[l][l],1,&x[j][l],1);
					if (j>=pl && j<=pu && qraux[j]!=0.0) {
						tt = ABS(x[j][l])/qraux[j];
						tt = 1.0f-tt*tt;
						tt = (float)MAX(tt,0.0);
						t = tt;
						ttt = qraux[j]/work[j];
						tt = 1.0f+0.05f*tt*ttt*ttt;
						if (tt!=1.0) {
							qraux[j] *= (float)sqrt(t);
						} else {
							qraux[j] = snrm2(n-l-1,
								&x[j][l+1],1);
							work[j] = qraux[j];
						}
					}
				}
				
				/* save the transformation */
				qraux[l] = x[l][l];
				x[l][l] = -nrmxl;
			}
		}
	}
}	
							
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)
/*****************************************************************************
Use 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.
******************************************************************************
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.
******************************************************************************
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;
	
	/* 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;
}