#include
#define LOC(a,b) ((b) * p + (a))
/*
function gs performs a Gram - Schmidt orthogonalization
of an n x p matrix, x. The orthogonal (q) part of the
decomposition is returned in x, and the lower triangular
part in r. (Both x and r must be allocated by the calling
program, and r should be filled with zeroes.) The matrix
x is assumed to be stored by columns.
*/
void gs(x,r,n,p)
double *x,*r;
long n,p;
{
long i,j,k;
double t;
double *xnow;
for(j=0;j