/*Copyright (C) 1992, 1996 by Thomas Glen Smith.  All Rights Reserved.*/
/* regrest APL2 V1.0.0 *************************************************
* Called from regress when memory is allocated.                        *
***********************************************************************/
#define INCLUDES APLCB+APLMEM
#include "includes.h"
double *regrest(left,rite,lrows,lcols,rrows,rcols,beta,xpx,xpy)
Aplcb left,rite,xpx,xpy;
double *beta;
int lrows,lcols,rrows,rcols;
{
	Invert;
	extern int aplerr;
	int i,j,k;
	double *dp,*ret,*xp,*xpxi,*xr,*yp,*yr;

	xp = rite->aplptr.apldata; /* input x matrix */
	dp = xpx->aplptr.apldata; /* first output location */
	for(i=0; i<rcols; i++) /* obtain XPX = X'*X */
		for(j = 0; j < rcols; j++) {
			*dp = 0.0; /* zero next output location */
			for(k = 0; k < rrows; k++) {
				xr = xp + k * rcols; /* proper row */
				*dp += *(xr + i) * *(xr + j); }
			dp++; /* point to next output location */ }
	xp = rite->aplptr.apldata; /* input x matrix */
	if (left == NULL) yp = NULL; /* use identity matrix for left */
	else yp = left->aplptr.apldata; /* input y matrix */
	dp = xpy->aplptr.apldata; /* first output location */
	for(i = 0; i < rcols; i++) /* calculate XPY = X'*Y */
		for(j = 0; j < lcols; j++) {
			*dp = 0.0; /* zero next output location */
			for(k = 0; k < lrows; k++) {
				xr = xp + k * rcols; /* proper row in x */
				if (left == NULL) *dp += *(xr + i) * (j == k);
				else {
					yr = yp + k * lcols; /* proper row in y */
					*dp += *(xr + i) * *(yr + j);
			}	}
			dp++; /* point to next output location */ }
	xpxi = invert(xpx->aplptr.apldata,rcols); /* obtain inverse of xpx */
	if (aplerr == 0) { /* now obtain BETA = XPXI * XPY */
		xp = xpxi; /* left matrix */ 
		yp = xpy->aplptr.apldata; /* right matrix */
		dp = beta; /* first output location */
		for(i = 0; i < rcols; i++) {
			xr = xp + i * rcols; /* proper row in XPXI */
			for(j = 0; j < lcols; j++) {
				yr = yp + j; /* proper col in XPY */
				*dp = 0.0; /* zero next output location */
				for(k = 0; k < rcols; k++)
					*dp += *(xr + k) * *(yr + k*lcols);
				dp++; /* point to next output location */
	}	}	}
	if (xpxi != NULL) free(xpxi); 
	endoper(xpx); endoper(xpy);
	if (aplerr) { free(beta); beta=NULL; }
	return(beta);
}