/*
   LSReg.c

   Here is a program that solves for two least-squares
   regression coefficients based on four experimentally
   designed tests.  It uses Matrix.library.

   This code is FREELY DISTRIBUTABLE but NOT PUBLIC DOMAIN
   (written by Randy Finch)
*/

#include "Test.h"

/*---------------
      main()
-----------------*/

void main(argc,argv)
   int argc;
   char **argv;
{
   DOUBLE *yield;	/* RHS vector */
   DOUBLE *yieldinter;	/* RHS intermediate vector */
   DOUBLE *b;		/* Coefficient vector */
   DOUBLE *predyield;	/* Predicted yields vector, data x b */
   DOUBLE **data;	/* Data matrix */
   DOUBLE **datatr;	/* Transpose of data */
   DOUBLE **data2;	/* data-transpose x data */
   DOUBLE **data2inv;	/* Inverse of data2 */
   DOUBLE **testinv;	/* Check the inverse */

   puts("\nOpening Matrix.library\n");
   Matrix = OpenLibrary("Matrix.library",1);
   printf("\nMatrix = %x\n",Matrix);

   if(!Matrix) {
      printf("\nCould not open Matrix.library!\n");
      exit(0);
   }

   /* Allocate necessary vectors and matrices */
   if (!(yield = (DOUBLE *)AllocDVector(4))) {
	 printf("\nCould not allocate YIELD vector.\n");
	 exit(0);
   }

   if(!(data = (DOUBLE **)AllocDMatrix(4,3))) {
	 printf("\nCould not allocate DATA matrix.\n");
	 exit(0);
   }

   /* Assign values to yield vector and data matrix */
   yield[0] = 40.5;
   yield[1] = 49.9;
   yield[2] = 61.8;
   yield[3] = 78.2;

   (data[0])[0] = 1.0;
   (data[0])[1] = 50.0;
   (data[0])[2] = 5.0;
   (data[1])[0] = 1.0;
   (data[1])[1] = 70.0;
   (data[1])[2] = 5.0;
   (data[2])[0] = 1.0;
   (data[2])[1] = 50.0;
   (data[2])[2] = 9.0;
   (data[3])[0] = 1.0;
   (data[3])[1] = 70.0;
   (data[3])[2] = 9.0;

   /* In the following statements DATATR, DATA2, YIELDINTER,
      DATA2INV, TESTINV, B, and PREDYIELD are allocated by
      the functions */

   /* Calculate the transpose of DATA */
   datatr = TransposeDMatrix(data,NULL,4,3);

   /* Multiply DATATR and DATA */
   data2 = MultDMatrices(datatr,data,NULL,3,4,3);

   /* Multiply DATATR and YIELD to get YIELDINTER */
   yieldinter = MultDMatrixVector(datatr,yield,NULL,3,4);

   /* Calculate the inverse of (DATA-transpose x DATA) */
   data2inv = InvertDMatrix(data2,NULL,NULL,3);
   if (!data2inv) {
      printf("\nMatrix inversion of DATA2 failed.\n");
      exit(0);
   }

   /* Multiply DATA2INV and DATA2 for check */
   testinv = MultDMatrices(data2inv,data2,NULL,3,3,3);

   /* Multiply DATA2INV and YIELDINTER to get B */
   b = MultDMatrixVector(data2inv, yieldinter, NULL, 3, 3);

   /* Multiply DATA and B to get PREDYIELD, compare with YIELD */
   predyield = MultDMatrixVector(data, b, NULL, 4, 3);

   /* Display all the vectors and matrices */
   printf("\nMatrix DATA:");
   DisplayDMatrix(data,4,3);
   printf("\nVector YIELD:");
   DisplayDVector(yield,4);
   printf("\nMatrix DATA-transpose:");
   DisplayDMatrix(datatr,3,4);
   printf("\nMatrix DATA2 (DATA-transpose x DATA):");
   DisplayDMatrix(data2,3,3);
   printf("\nVector YIELDINTER (DATA-transpose x YIELD):");
   DisplayDVector(yieldinter,3);
   printf("\nMatrix DATA2-inverse:");
   DisplayDMatrix(data2inv,3,3);
   printf("\nMatrix TESTINV (DATA2INV x DATA2):");
   DisplayDMatrix(testinv,3,3);
   printf("\nCoefficient vector B:");
   DisplayDVector(b,3);
   printf("\nVector PREDYIELD (= DATAxB) Compare with YIELD:");
   DisplayDVector(predyield,4);

   /* Free vectors and matrices */
   FreeDVector(yield,4);
   FreeDVector(yieldinter,3);
   FreeDVector(b,3);
   FreeDVector(predyield,4);
   FreeDMatrix(data,4,3);
   FreeDMatrix(datatr,3,4);
   FreeDMatrix(data2,3,3);
   FreeDMatrix(data2inv,3,3);

   /* Finally, close the library */
   CloseLibrary(Matrix);

} /* main */