/*
   SimEqs.c

   Here is a test program that solves three simultaneous
   equations using 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 *x;	     /* Solution vector */
   DOUBLE *b;	     /* RHS vector */
   DOUBLE *btest;    /* AX test vector */
   DOUBLE **a;	     /* Coefficient matrix */
   DOUBLE **ainv;    /* Inverse of a */

   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 (!(b = (DOUBLE *)AllocDVector(3))) {
	 printf("\nCould not allocate B vector.\n");
	 exit(0);
   }

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

   /* Assign values to B vector and A matrix */
   b[0] = 1003.88;
   b[1] = -667.05;
   b[2] = 3.0;

   (a[0])[0] = 1.3;
   (a[0])[1] = 342.7;
   (a[0])[2] = -17.7;
   (a[1])[0] = -222.17;
   (a[1])[1] = 3.0;
   (a[1])[2] = 23.6;
   (a[2])[0] = 1.0;
   (a[2])[1] = -17.4;
   (a[2])[2] = -15.6;

   /* In the following three statements AINV, X, and BTEST are
      allocated by the functions */

   /* Calculate the inverse of A */
   ainv = InvertDMatrix(a,NULL,NULL,3);
   if (!ainv) {
      printf("\nMatrix inversion of A failed.\n");
      exit(0);
   }

   /* Multiply A inverse and B to get X */
   x = MultDMatrixVector(ainv, b, NULL, 3, 3);

   /* Multiply A and X to get B-test, compare with B */
   btest = MultDMatrixVector(a, x, NULL, 3, 3);

   /* Display all the vectors and matrices */
   printf("\nMatrix A:");
   DisplayDMatrix(a,3,3);
   printf("\nVector B:");
   DisplayDVector(b,3);
   printf("\nMatrix A-inverse:");
   DisplayDMatrix(ainv,3,3);
   printf("\nSolution vector X:");
   DisplayDVector(x,3);
   printf("\nVector B-test (= AX) Compare with B:");
   DisplayDVector(btest,3);

   /* Free vectors and matrices */
   FreeDVector(x,3);
   FreeDVector(b,3);
   FreeDVector(btest,3);
   FreeDMatrix(a,3,3);
   FreeDMatrix(ainv,3,3);

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

} /* main */