/* test_inverse_small.c   small matrix, little input error */
/*                        big difference in output, not seen in error */

#include <stdio.h>
#include "inverse.h"
#define  abs(a) ((a)<0.0?(-(a)):(a))

static double AA[2][2];
static double AI[2][2];
static double II[2][2];

int main(int argc, char * argv[])
{
  int i, j, k, n;
  double err;

  printf("test_inverse_small.c on  \n");
  n = 2;

  printf("err is the error multiplying the solution AI times A=I \n");

  printf("decimal | .835  .667 \n");
  printf("decimal | .333  .266 \n");
  AA[0][0] = 0.835; /* row 0 */
  AA[0][1] = 0.667;
  AA[1][0] = 0.333; /* row 1 */
  AA[1][1] = 0.266;
  printf("flt-pt | %f  %f | \n", AA[0][0], AA[0][1]);
  printf("flt-pt | %f  %f | \n", AA[1][0], AA[1][1]);
  printf("initialization complete, solving \n");


  inverse2(n, AA , AI );
  printf("inverse \n");
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
    {
      printf("AI[%1d][%1d]=%g \n", i, j, AI[i][j]);
    }
  }
 
  err = 0.0;
  for(i=0; i<n; i++)
  {
    for(j=0; j<n; j++)
    {
      II[i][j] = 0.0;
      for(k=0; k<n; k++)
      {
        II[i][j] = II[i][j] + AA[i][k]*AI[k][j];
      }
      if(i == j)
        err = err + abs(II[i][j]-1.0);
      else
        err = err + abs(II[i][j]);
    }
    printf("sum of error=%g, avg error=%g \n", err, err/(double)(n*n));
  }

  printf("finished test_inverse_small.c \n");
  return 0;
} /* end test_inverse_small.c */