// test_simeq_hilbert.c  very difficult as size increases Hilbert matrices
#include <stdio.h>
#include <stdlib.h>
// #include "simeq2.h"
#include "simeq2p.h"

static int sz = 100;

int main(int argc, char *argv[])
{
  double *A; // [sz*sz];
  double *AA; // [sz*sz]; // saved
  double Y[sz];
  double YY[sz];    // saved
  double X[sz];
  double err;
  int i, j;
  int n;

  printf("test_simeq_hilbert.c \n");
  printf("err is the maximum error multiplying the solution, X, times matrix A*X=Y \n");
  for(n=3; n<sz; n=n+4)
  {
    A = (double *)malloc(n*n*sizeof(double));
    AA = (double *)malloc(n*n*sizeof(double));
    printf("initialization of Hilbert Matrix A[%d][%d]\n", n, n);
    for(i=0; i<n; i++)
    {
      Y[i] = 0.0;
      for(j=0; j<n; j++)
      {
        A[i*n+j] = 1.0/((double)(i+1) + (double)(j+1));
        AA[i*n+j] = A[i*n+j];
        Y[i] = Y[i] + A[i*n+j]*(double)(j+1);
      }
      YY[i] = Y[i];
    }

    simeq2( n, A , Y , X );

    for(i=0; i<n; i++) // restore A and Y
    {
      for(j=0; j<n; j++)
      {
        A[i*n+j] = AA[i*n+j];
        if(n<4) printf("A[%d][%d] = %f \n", i, j, A[i*n+j]);
      }
      Y[i] = YY[i];
      if(n<4) printf("Y[%d] = %f \n", i, Y[i]);
    }

    for(i=0; i<n; i++)
    {
      err = -Y[i]; // -(double)(i+1); // solution 1.0, 2.0, ...
      for(j=0; j<n; j++)
      {
        err = err + A[i*n+j]*X[j];
      } 
      if(i<5)printf("solution X[%1d]=%g, err=%g \n", i, X[i], err);
    }
    printf("\n");
    free(A);
    free(AA);
  }
  printf("\n test_simeq_hilbert.c finished");

  return 0;
} // end test_simeq_hilbert.c