/* mat-test.c test routines for the mat package.
 * (C) Copyright 1999 by John Halleck
 * All rights reserved.
 */
/* Version of September 6th, 1999 */

#include <stdio.h>
  /* Standard IO */

#include "errors.h"
  /* Common error codes */
#include "mat.h"
  /* Matrix specific stuff */

#include "matdebug.h"
  /* Standard matrix debugging routines. */

#include "reportframe.h"
  /* Standard test reports */

double test1[1][1];
double test5[5][5];

double R[4][2], A[4][3], B[3][2], BT[2][3];
double TR[4][2], TA[4][3], TB[3][2], TBT[2][3];
double TT[4][2];

/* ========================== MAIN ========================================== */

int main (argc, argv)
  int argc; char argv[];
{
  error problem   = 0;

  int i,j; /* Loop indices */

  double memtestarray[4][5], memtestaref[4][5];

  double identtest[6][6], identref[6][6];

  double size5a[5], size5b[5];

  double transtesta[2][3], transtestf[3][2], transref[3][2];

  size5a[0] = 1.0; size5a[1] = 2.0; size5a[2] = 3.0; size5a[3] = 4.0;
  size5a[4] = 5.0;
  size5b[0] = 1.0; size5b[1] = 10.0; size5b[2] = 10.0; 
  size5b[3] = 10.0; size5b[4] = 5.0;

  inittests("mat");

/* ======================== Routines to load constants ====================== */

  newsection ("Constant loads");

/* ------------------------------ test matcon ------------------------------- */

/* matzero (set a matrix to zero) */
  newtest("matzero");

  for (i=0;i<4;i++) for (j=0;j<5;j++) memtestaref[i][j]=0.0;
  if (   ERRsize != (problem = matzero (0, 5, memtestarray))
      || ERRsize != (problem = matzero (4, 0, memtestarray))
      )
     goterrorstat ("bad size accepted?", problem);
  if (ERRnil != (problem = matzero (4, 5, 0))) 
     goterrorstat ("nil pointer accepted?", problem);

  if ((problem = matzero (4, 5, memtestarray)))
     goterrorstat ("failed valid call?", problem);
  if ((problem = matiseq(4, 5, memtestarray, memtestaref))) {
     if (problem == ERRfalse) problem = NoError;
     goterrorstat ("incorrect fillin?", problem);
     matprint ("Result matrix", 4, 5,  memtestarray);
  }

  endtest();
/* ----------------------------- test matidn -------------------------------- */

/* matidn (set a matrix to the identity) */

  newtest ("matidn");

  for (i=0;i<6;i++) for (j=0;j<6;j++) {
     identtest[i][j] = 12.0; identref[i][j] = 0.0; 
  }
  for (i=0;i<6;i++) identref[i][i] = 1.0;

  if (ERRsize != (problem = matidn (0, identtest)))
     goterrorstat ("bad size accepted?", problem);
  if (ERRnil != (problem = matidn (6, 0)))
     goterrorstat ("nil pointer accepted?", problem);

  if ((problem = matidn (6, identtest)))
     goterrorstat ("valid call failed?", problem);

  if ((problem = matiseq(6, 6, identtest, identref))) {
     if (problem == ERRfalse) problem = NoError;
     goterrorstat ("incorrect fillin?", problem);
     matprint ("Result matrix", 6, 6, identtest);
     matprint ("Should be",     6, 6, identref);
  }

  endtest ();
/* ===================== End of load constants tests ======================== */

endsection();

newsection ("Copy routines");

/* ===================== Copy routines ====================================== */

/* --- matcpy --- */
  newtest("matcpy");

  for (i=0;i<6;i++) for (j=0;j<6;j++) identref[i][j] = (float) (6*i+j);
  for (i=0;i<6;i++) for (j=0;j<6;j++) identtest[i][j] = 9999.0; 

  if (ERRsize != (problem = matcpy (0, 6, identref, identtest))
   || ERRsize != (problem = matcpy (6, 0, identref, identtest)) )
     goterrorstat ("bad size accepted?", problem);
  if (ERRnil != (problem = matcpy (6, 6,        0, identtest))
   || ERRnil != (problem = matcpy (6, 6, identref,         0)) )
     goterrorstat ("nil pointer accepted?", problem);
  if ((problem = matcpy (6, 6, identref, identtest)))
     goterrorstat ("valid call accepted?", problem);
  if ((problem = matiseq (6, 6, identtest, identref))) { 
     if (problem == ERRfalse) problem = NoError;
     goterrorstat ("wrong result?", problem);
     matprint ("was",       6, 6, identtest);
     matprint ("should be", 6, 6, identref);
  }

  endtest();

/* --- matcpyt ---- */
/* Matrix transpose copy test. */
  newtest ("matcpyt");

  transtesta[0][0] = transref[0][0] = 1;
  transtesta[0][1] = transref[1][0] = 2;
  transtesta[0][2] = transref[2][0] = 3;
  transtesta[1][0] = transref[0][1] = 4;
  transtesta[1][1] = transref[1][1] = 5;
  transtesta[1][2] = transref[2][1] = 6;

  if (ERRsize != (problem = matcpyt (0, 2, transtestf, transtesta))
   || ERRsize != (problem = matcpyt (3, 0, transtestf, transtesta)))
     goterrorstat ("bad size accepted?", problem);
  if (ERRnil != (problem = matcpyt (3, 2,          0, transtesta))
   || ERRnil != (problem = matcpyt (3, 2, transtestf,          0)) )
     goterrorstat ("nil pointer accepted?", problem);
  if ((problem = matcpyt (3, 2, transtestf, transtesta)))
     goterrorstat ("valid call rejected?", problem);

  if ((problem = matiseq (3, 2, transtestf, transref))) {
     if (problem == ERRfalse) problem = NoError;
     goterrorstat ("wrong result?", problem);
     matprint ("Transpose of", 2, 3, transtesta);
     matprint ("was",       3, 2, transtestf);
     matprint ("should be", 3, 2, transref);
  }

  endtest();

/* ===================== End of copy tests ================================== */

endsection();

/* ===================== Scalar Multiply tests ============================== */

  newsection ("Scalar Multiply");

   newtest("matscl");
   if ((problem = matcpy (4, 2, R, TR)))
      goterrorstat ("*** matcpy failed in matmuls test?", problem);
   if ( ERRnil != (problem = matscl (4, 2, 0, TR,  2.0))
     || ERRnil != (problem = matscl (4, 2, TR, 0,  2.0)) )
      goterrorstat ("nil array accepted?", problem);
   if ( ERRsize != (problem = matscl (0, 2, TR, TR, 2.0))
     || ERRsize != (problem = matscl (4, 0, TR, TR, 2.0)) )
      goterrorstat ("bad size accepted?", problem);
   if ( ERRmeaning != (problem = matscl (4, 2, TR, TR, 0.0)) )
      goterrorstat ("meaningless scaling accepted", problem);

   if ((problem = matscl (4, 2, TR, TR, 2.0)))
      goterrorstat ("valid call rejected?", problem);
   else {
     for (i=0;i<4;i++) {
         for (j=0;j<2;j++) {
            if (TR[i][j] != R[i][j]*2.0) {
               if (!testerrors) {
                  goterror ("*** TR != R*2.0");
                  matprint ("R=", 4, 2, R);
                  matprint ("TR=", 4, 2, R);
               }
            }
         }
     }
   }

   endtest();

endsection();

/* ===================== Multiply tests. ==================================== */

   newsection ("Multiply");

   newtest ("matmul");
   if ( ERRsize != (problem = matmul (0, 2, R, 3, A, B))
     || ERRsize != (problem = matmul (4, 0, R, 3, A, B))
     || ERRsize != (problem = matmul (4, 2, R, 0, A, B)))
        goterrorstat ("bad size accepted?", problem);
   if ( ERRnil != (problem = matmul (4, 2, 0, 3, A, B))
     || ERRnil != (problem = matmul (4, 2, R, 3, 0, B))
     || ERRnil != (problem = matmul (4, 2, R, 3, A, 0))) 
       goterrorstat ("nil pointer accepted?", problem);
   if ( ERRsame != (problem = matmul (4, 2, R, 3, R, B))
     || ERRsame != (problem = matmul (4, 2, R, 3, A, R)))
        goterrorstat ("non-distinct arrays accepted?", problem);
   for (i=0;i<4;i++) for (j=0;j<2;j++) R[i][j] = 99999.0;
   A[0][0] =  1.0;   A[0][1] =  2.0;   A[0][2] =  3.0;
   A[1][0] =  4.0;   A[1][1] =  5.0;   A[1][2] =  6.0;
   A[2][0] =  7.0;   A[2][1] =  8.0;   A[2][2] =  9.0;
   A[3][0] = 10.0;   A[3][1] = 11.0;   A[3][2] = 12.0;
   B[0][0] =  1.0;   B[0][1] =  2.0;
   B[1][0] =  3.0;   B[1][1] =  4.0;
   B[2][0] =  5.0;   B[2][1] =  6.0;
   if ((problem = matmul (4, 2, R, 3, A, B))) 
      goterrorstat ("valid call rejected?", problem);
   else {
      if (R[0][0] !=  22.0 || R[0][1] !=  28.0 
       || R[1][0] !=  49.0 || R[1][1] !=  64.0
       || R[2][0] !=  76.0 || R[2][1] != 100.0
       || R[3][0] != 103.0 || R[3][1] != 136.0
      ) {
       goterror ("Wrong Answer");
       matprint ("A*B = ", 4, 2, R);
       printf ("*** Should have been:\n");
       printf ("22 28\n49 64\n76 100\n103 136\n");
      }
   }
   endtest();

   BT[0][0] = 1.0;     BT[0][1] = 3.0;   BT[0][2] = 5.0;
   BT[1][0] = 2.0;     BT[1][1] = 4.0;   BT[1][2] = 6.0;

   newtest ("matmult");

   if ( ERRsize != (problem = matmult (0, 2, R, 3, A, BT)) 
     || ERRsize != (problem = matmult (4, 0, R, 3, A, BT))
     || ERRsize != (problem = matmult (4, 2, R, 0, A, BT))
     ) goterrorstat ("bad size accepted", problem);
   if ( ERRnil != (problem = matmult (4, 2, 0, 3, A, BT))
     || ERRnil != (problem = matmult (4, 2, R, 3, 0, BT))
     || ERRnil != (problem = matmult (4, 2, R, 3, A,  0))
      ) goterrorstat ("nil pointer accepted?", problem);
   if ( ERRsame != (problem = matmult (4, 2, R, 3, R, BT))
     || ERRsame != (problem = matmult (4, 2, R, 3, A,  R))
     ) goterror ("non-distinct arrays accepted");
   if ((problem = matmult (4, 2, R, 3, A, BT)))
     goterrorstat ("valid call rejected?", problem);
   else {
      if (R[0][0] !=  22.0 || R[0][1] !=  28.0
       || R[1][0] !=  49.0 || R[1][1] !=  64.0
       || R[2][0] !=  76.0 || R[2][1] != 100.0
       || R[3][0] != 103.0 || R[3][1] != 136.0
      ) {
       goterror ("multiplication by transpose has problems???");
       matprint ("A * B' = ", 4, 2, R);
       printf ("*** Should have been:\n");
       printf ("22 28\n49 64\n76 100\n103 136\n");
      }
   }

   endtest();

   /* ----------- */

endsection();

newsection ("addition");

   newtest("matadd");

   if ((problem = matcpy (4, 2, TR, R)))
      goterrorstat ("*** matcpy failed valid call", problem);
   if ( ERRsize != (problem = matadd (0, 2, TT, TR, R))
     || ERRsize != (problem = matadd (4, 0, TT, TR, R)) )
      goterrorstat ("bad size accepted?", problem);
   if ( ERRnil != (problem = matadd(4,2, 0,TR,R))
     || ERRnil != (problem = matadd(4,2,TT, 0,R))
     || ERRnil != (problem = matadd(4,2,TT,TR,0)) )
      goterrorstat ("nil array accepted?", problem);

   if ((problem = matadd (4, 2, TT, TR, R)))
      goterrorstat ("valid call rejected?", problem);
   else {
     for (i=0;i<4;i++) {
         for (j=0;j<2;j++) {
            if (TT[i][j] != R[i][j]*2.0) {
               if (!testerrors) {
                  goterror ("TT != R*3.0");
                  matprint ("R=", 4, 2, R);
                  matprint ("TR=", 4, 2, TR);
                  matprint ("TT=", 4, 2, TT);
               }
            }
         }
     }
   }

  endtest();

  newtest ("matsub");
   for (i=0;i<4;i++) for (j=0;j<2;j++) {
      R[i][j]  = (i+1)*100 + j+1;
      TR[i][j] = R[i][j] * 2;
      TT[i][j] = 0;
   }
   if ( ERRsize != (problem = matsub (0, 2, TT, TR, R))
     || ERRsize != (problem = matsub (4, 0, TT, TR, R)) )
      goterrorstat ("bad size accepted?", problem);
   if ( ERRnil != (problem = matsub(4,2, 0,TR,R))
     || ERRnil != (problem = matsub(4,2,TT, 0,R))
     || ERRnil != (problem = matsub(4,2,TT,TR,0)) )
      goterrorstat ("nil array accepted?", problem);

   if ((problem = matsub (4, 2, TT, TR, R)))
      goterrorstat ("valid call rejected?", problem);
   else {
     for (i=0;i<4;i++) {
         for (j=0;j<2;j++) {
            if (TT[i][j] != R[i][j]) {
               if (!testerrors) {
                  goterror ("TT != R*3.0");
                  matprint ("R=", 4, 2, R);
                  matprint ("TR=", 4, 2, TR);
                  matprint ("TT=", 4, 2, TT);
               }
            }
         }
     }
   }
  endtest();

endsection();


finalizetests();
return progerrors;
}