/* qinv-test.c - Routine to test the qinv routines.
 * (C) Copyright 1999, by John Halleck
 * All rights reserved.
 */
/* Version of August 14th, 1999 */

#include <stdio.h>
 /* We have print problem messages  */

#include "errors.h"
 /* We have to print final problem messages */

#include <math.h>
 /* Trig functions can be used to build good tests. */

#include "matdebug.h"
 /* The matrix debugging routines would be handy */

#include "reportframe.h"
 /* Standard form of our test output */

#include "qinv.h"
 /* And we need, of course, the routines we were testing. */

#define FUDGE (1.0 / 100000.0)
/* How close to arrays have to be to "pass" the inverse**2 = identity
 * tests on this machine.
 */

static double result[3][3], given[3][3], working[3][3];
double check[3][3]; /* This must be as large as the largest we test with */
static int problem;

static double S, C; 
static double  trig1[2][2], trig2[2][2], tfinal[2][2], twork[2][2];
/* Trig functions make the nicest tests... */

static int testone (int size, matrix result, matrix given, matrix working) {
  if ( ERRsize != (problem = invpd (   0, result, given, working)) )
       goterrorstat ("bad size passed", problem);
  if ( ERRnil != (problem = invpd (size,      0, given, working))
    || ERRnil != (problem = invpd (size, result,     0, working))
    || ERRnil != (problem = invpd (size, result, given,       0))
      )
     goterrorstat ("nil array passed?", problem);
  if ( ERRsame != (problem = invpd (size, result,  given, result)) )
     goterrorstat ("array alias problems passed?", problem);

  if ((problem = invpd (size, result, given, working)))
     goterrorstat ("valid call failed", problem);

  if ((problem = invpd (size, check, result, working)))
     goterrorstat ("call failed on valid arguments.", problem);
 
  if ((problem = matisaeq (size, size, given, check, FUDGE))) {
     if (ERRfalse == problem) problem = NoError;
     goterrorstat ("inverse of inverse is not original", problem);
     matprint ("Original", size, size, given);
     matprint ("Inverse", size, size, result);
     matprint ("Inverse of Inverse", size, size, check);
  }

  return problem;
}

int main() {

  inittests("qinv");

  newtest("invpd");

  if ((problem = matidn (3, result)) || (problem = matidn (3, given)))
     goterrorstat ("matidn returned bad problem?", problem);
  testone (3, result, given, working); 
  if ((problem = matiseq(3, 3, given, result))) {
    /* Note that this is an EXACT equality test */
     if (ERRfalse == problem) problem = NoError;
     goterrorstat ("can't get exact inverse of identity?", problem);
     printf ("    (Either it is wrong or arithmetic sucks on this machine.)\n");
     matprint ("given",    3, 3, given);
     matprint ("working:", 3, 3, working);
  }

  S = sin (3.1415926/180.0 * 5);
  C = sqrt (1.0 - S*S);
  trig1[0][0] =  C;   trig1[0][1] = S;
  trig1[1][0] = -S;   trig1[1][1] = C;
  testone (2, tfinal, trig1, twork);
  trig2[0][0] =  C;   trig2[0][1] = -S;
  trig2[1][0] =  S;   trig2[1][1] =  C;
  if ((problem = matisaeq(2, 2, trig2, tfinal, FUDGE))) {
     if (ERRfalse == problem) problem = NoError;
     goterrorstat ("rotation inverse is not inverse", problem);
     matprint ("invpd's version:", 2, 2, tfinal);
     matprint ("real version:", 2, 2, trig2);
  }

  endtest();

/* -------------------------- invns ----------------------------------------- */
/* Test the pivoting version of the above. */

  newtest("invns");
  if ( ERRsize != (problem = invns (   0, result, given, working)) )
     goterrorstat ("bad size accepted", problem);
  if ( ERRnil != (problem = invns (3,      0, given, working))
    || ERRnil != (problem = invns (3, result,     0, working))
    || ERRnil != (problem = invns (3, result, given,       0))
      )
     goterrorstat ("didn't notice nil array?", problem);
  if ( ERRsame != (problem = invpd (3, result,  given, result)) )
     goterrorstat ("alias problems ignored?", problem);

  given[0][0] = 0.0;  given[0][1] = 0.0;  given [0][2] = 1.0;
  given[1][0] = 0.0;  given[1][1] = 1.0;  given [1][2] = 0.0;
  given[2][0] = 1.0;  given[2][1] = 0.0;  given [2][2] = 0.0;

  if ((problem = matcpy (3, 3, check, given)))
     goterrorstat ("matcpy failed during invns test?", problem);

  if ((problem = invns (3, result, given, working)))
     goterrorstat ("rejected valid call", problem);

  if ((problem = matiseq (3, 3, result, check))) {
     if (problem == ERRfalse) problem = NoError;
     goterrorstat ("wrong answer?", problem);
     matprint ("Given:", 3, 3, given);
     matprint ("Answer:", 3, 3, result);
  }

  given[1][1] = 0.0;
  if (ERRnumeric != (problem = invns (3, result, given, working))) {
     goterrorstat ("accepted obviously singular matrix call", problem);
     matprint ("Given:",   3, 3, given);
     matprint ("Working:", 3, 3, working);
     matprint ("Answer:",  3, 3, result);
  }

  given[1][1] = 1.0; given [1][2] = 1.0; given [0][1] = 1.0;
  if (ERRnumeric != (problem = invns (3, result, given, working))) {
     goterrorstat ("accepted subtle singular matrix call", problem);
     matprint ("Given:",   3, 3, given);
     matprint ("Working:", 3, 3, working);
     matprint ("Answer:",  3, 3, result);
  }

/* ========================================================================== */

  endtest();
  finalizetests();
  return progerrors;
}