/* test_mpf_inverse.c  tough matrix to invert accurately, use multiple */
/*                     precision math package                          */

#include <stdlib.h>
#include "gmp.h"
#include "mpf_inverse.h"

#define digits 200

static mpf_t AA[65535];
static mpf_t AI[65536];
static mpf_t II[65536];

int main(int argc, char * argv[])
{
  int i, j, k, n;
  mpf_t err;
  mpf_t JK, T, T2, zero, one;

  gmp_printf("test_mpf_inverse.c on tough matrix to invert, %d digits \n",
             digits);
  mpf_set_default_prec(digits*3.32); /* digits */
  mpf_init(err);
  mpf_init(JK);
  mpf_init(T);
  mpf_init(T2);
  mpf_init(zero);
  mpf_init(one);

  n = 1;
  mpf_set_si(zero, 0);
  mpf_set_si(one, 1);
  while(n < 128)  /* < 1024 */
  {
    n=n*2;
    gmp_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++)
      {
        mpf_init(II[i*n+j]);
        mpf_init(AI[i*n+j]);
        mpf_init(AA[i*n+j]);
        mpf_set_si(JK,(i+1)+(j+1));
        mpf_div(T, one, JK);
        mpf_set(AA[i*n+j], T);
      }
    }
    if(n<=8)
    {
      gmp_printf("data \n");
      for(i=0; i<n; i++)
      {
        for(j=0; j<n; j++)
        {
          gmp_printf("AA[%d][%d]=%Fg \n", i, j, AA[i*n+j]);
        }
      }
    } 
    mpf_inverse( n, digits, AA , AI);
    if(n<=8)
    {
      gmp_printf("inverse \n");
      for(i=0; i<n; i++)
      {
        for(j=0; j<n; j++)
	{
          gmp_printf("AI[%d][%d]=%Fg \n", i, j, AI[i*n+j]);
        }
      }
    } 
    mpf_set_si(err, 0);
    for(i=0; i<n; i++)
    {
      for(j=0; j<n; j++)
      {
        mpf_set(II[i*n+j], zero);
        for(k=0; k<n; k++)
	{
          mpf_mul(T, AA[i*n+k], AI[k*n+j]);
          mpf_add(II[i*n+j], II[i*n+j], T);
        }
        if(i == j)
	{
          mpf_sub(T, II[i*n+j], one);
          mpf_abs(T, T);
	  mpf_add(err, err, T); 
        }
        else
	{
          mpf_abs(T, II[i*n+j]);
          mpf_add(err, err, T);
	}
      }
    }
    if(n<=8)
    {
      gmp_printf("identity \n");
      for(i=0; i<n; i++)
      {
        for(j=0; j<n; j++)
	{
          gmp_printf("II[%d][%d]=%Fg \n", i, j, II[i*n+j]);
        }
      }
    } 

    gmp_printf("sum of errors=%Fg \n", err);
  }
  gmp_printf("finished test_mpf_inverse.c \n");
  return 0;
} /* end test_mpf_inverse.c */