/* test_sparse.c */    

#include <stdio.h>
#include "sparse.h"
#define min(a,b) ((a)<(b)?(a):(b))
#define abs(a)   ((a)<0.0?(-(a)):(a))

#define nrow 6

int main(int argc, char * argv[])
{
  struct sparse A;    /* the sparse simultaneous equations */
                      /* normal A matrix with Y as last column */
  struct sparse Ac;   /* copy for checking */
  double val;
  double X[nrow];  /* solution of A*X=Y simultaneous equations */
  double Y[nrow];  /* extracted Y */
  double Yc[nrow]; /* for checking */
  double AA[nrow][nrow]; /* for testing */
  int i, j;

  printf("test_sparse.c nrow=%d \n",nrow);
  sparseInit(&A, nrow);
  sparseInit(&Ac, nrow);

  printf("test put(0,4,3.14) and get in A \n");
  sparsePut(A,0,4,3.14);
  val = sparseGet(A,0,4);
  printf("test put in row, 3.14 should be %g \n",val);
  printf("test put(5,0,6.28) and get in A \n");
  sparsePut(A,5,0,6.28);
  val = sparseGet(A,5,0);
  printf("test put in row, 6.28 should be %g \n",val);
  sparseWrite_All(A);
  sparseInit(&A, nrow);

  for(i=0; i<nrow; i++)
  {
    for(j=i-1; j<=i+1; j++)
    {
      val = 1.0+(double)(i*i+(j+1)*(j+1));
      if(i>=0 && j>=0 && i<nrow && j<nrow)
      {
        sparseAdd(A, i, j, val);
        sparseAdd(Ac,i, j, val);
      }
      sparseAdd(A, i, nrow, 8.0+4.0*i);
      sparseAdd(Ac,i, nrow, 8.0+4.0*i);
    }
  }
  printf("initial set of equations in A \n");
  sparseWrite_All(A);

  sparseGetRHS(A,Y);      /* for testing */
  for(i=0; i<nrow; i++)
  {
    printf("i=%d, Y[i]=%g \n",i,Y[i]);
  }


  printf("calling sparse.simeq \n");
  sparseSimeq(A,X); /* modifies, destroys, A */
  for(i=0; i<nrow; i++)
    printf("X[%d]=%g \n",i,X[i]);
  printf("A after simeq \n");
  sparseWrite_All(A);

  printf("Ac after simeq, should be initial \n");
  sparseWrite_All(Ac);

  printf("\n");
  printf("test multiply and simeq \n");
  sparseMul(Ac,X,Yc);

  for(i=0; i<nrow; i++)
  {
    if(abs(Y[i]-Yc[i])> 0.00001)
    {
      printf("either multiply or simeq failed \n");
      printf("i=%d, Y[i]=%g yet found \n",i,Yc[i]);
    }
  }

  printf("\n");
  printf("simple test of setRHS, getRHS \n");
  for(i=0; i<nrow; i++) Y[i] = (double)(i+1);
  sparseSetRHS(A,Y);
  sparseWrite_All(A);
  sparseGetRHS(A,Yc);
  for(i=0; i<nrow; i++)
  {
    if(Y[i]!=Yc[i])
    {
      printf("either setRHS or getRHS failed \n");
      printf("i=%d, Y[i]=%g yet found \n",i,X[i]);
    }
  }
  printf("test_sparse finished \n");
  return 0;
} /* end test_sparse.c */