/* test_inverse.c  tough matrix to invert accurately */
#include <stdio.h>
#include "inverse.h"
#define  abs(a) ((a)<0.0?(-(a)):(a))

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

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

  printf("test_inverse.c on tough matrix to invert \n");
  n = 1;
  while(n < 1024)
  {
    n=n*2;
    printf("initializing big matrix, n=%d, n*n=%d \n", n, n*n);
    for(i=0; i<n; i++)
    {
      for(j=0; j<n; j++)
      {
        AA[i*n+j] = 1.0/(double)((i+1)+(j+1));
      }
    }
    if(n<=4)
    {
      printf("data \n");
      for(i=0; i<n; i++)
      {
        for(j=0; j<n; j++)
        {
          printf("AA[%d][%d]=%g \n", i, j, AA[i*n+j]);
        }
      }
    } 
    inverse( n, AA , AI );
    if(n<=4)
    {
      printf("inverse \n");
      for(i=0; i<n; i++)
      {
        for(j=0; j<n; j++)
	{
          printf("AI[%d][%d]=%g \n", i, j, AI[i*n+j]);
        }
      }
    } 
    err = 0.0;
    for(i=0; i<n; i++)
    {
      for(j=0; j<n; j++)
      {
        II[i*n+j] = 0.0;
        for(k=0; k<n; k++)
	{
          II[i*n+j] = II[i*n+j] + AA[i*n+k]*AI[k*n+j];
        }
        if(i == j)
          err = err + abs(II[i*n+j]-1.0);
        else
          err = err + abs(II[i*n+j]);
      }
    }
    printf("AI[%d][%d]=%g \n", n, n, AI[(n-1)*n+(n-1)]);
    printf("sum of error=%g, avg error=%g \n", err, err/(double)(n*n));
  }
  printf("finished test_inverse.c \n");
  return 0;
} /* end test_inverse.c */