/* test1_sparse.c */

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


static void PrintM(double M[], int N)
{
  int i, j;
  for(i=0; i<N; i++)
    for(j=0; j<=N; j++)
      printf("i=%d, j=%d, M=%g \n", i, j, M[i*N+j]);
} /* end PrintM */

static void PrintV(double V[], int N)
{
  int i;
  for(i=0; i<N; i++)
    printf("i=%d , V=%g \n", i, V[i]);
} /* end PrintV */

int main(int argc, char * argv[])
{
  double Val = 0.0;
  double M[7*8]; /* multiple use */
  double Y[7];   /* multiple use */
  double Yc[7];  /* multiple use, check */
  double X[7];   /* multiple use */
  int Nrow, i, j;
  struct sparse A; /* size established by sparseInit */

  printf("test1_sparse.c \n");
  Nrow = 1;
  sparseInit(&A, Nrow); /* the sparse simultaneous equation, 1 */
  sparseWrite_All(A);
  printf("Init should have 0, 0 \n");
  printf("\n");

  sparsePut(A,0,0,2.0);
  sparsePut(A,0,1,6.0);
  sparseWrite_All(A);
  printf("Put should have 2, 6 \n");
  printf("\n");

  sparseExport(A,M);
  PrintM(M,Nrow);
  printf("Export should have M 2, 6 \n");
  printf("\n");

  sparseGetRHS(A,Y);
  PrintV(Y,Nrow);
  printf("GetRHS, should have V 6 \n");
  printf("\n");

  sparsePut(A,0,0,2.0);
  sparsePut(A,0,1,6.0);
  sparseWrite_All(A);
  printf("Put should have 2, 6 \n");
  printf("\n");

  sparsePut(A,0,1,7.0);
  sparseWrite_All(A);
  printf("Put should have 2, 7 \n");
  printf("\n");

  sparsePut(A,0,0,8.0);
  sparseWrite_All(A);
  printf("Put should have 8, 7 \n");
  printf("\n");

  sparseInit(&A, Nrow);
  M[0] = 5.0;
  M[1] = 7.0;
  sparseImport(A,M);
  sparseWrite_All(A);
  printf("Import over Init should have be 5,7 \n");
  printf("\n");

  M[0] = 9.0;
  M[1] = 11.0;
  sparseImport(A,M);
  sparseWrite_All(A);
  printf("Import over existing should have be 9,11 \n");
  printf("\n");

  Y[0] = 9;
  sparseSetRHS(A,Y);
  sparseWrite_All(A);
  printf("SetRHS, should have 9 \n");
  printf("\n");

  sparseInit(&A, Nrow); /* the sparse simultaneous equation, 1 */
  sparsePut(A,0,0,2.0);
  sparsePut(A,0,1,6.0);
  sparseWrite_All(A);
  printf("Put should have 2, 6 \n");
  printf("\n");

  sparseSimeq(A,X);
  PrintV(X,Nrow);
  printf("Simeq should be +3.0 \n");
  sparseWrite_All(A);
  printf("Simeq should have 2, 3 \n");
  printf("\n");

  Nrow = 2;
  sparseInit(&A, Nrow); /* the sparse simultaneous equation, 2 */
  printf("test1_sparse.c nrow=%d \n", Nrow);
  sparsePut(A,0,0,2.0);
  sparsePut(A,0,1,4.0);
  sparsePut(A,0,2,16.0);
  sparsePut(A,1,0,6.0);
  sparsePut(A,1,1,3.0);
  sparsePut(A,1,2,21.0);
  sparseWrite_All(A);
  sparseExport(A,M);
  PrintM(M,Nrow);
  sparseSimeq(A,X);
  PrintV(X,Nrow);
  printf("should be +2.0 +3.0 \n");
  sparseWrite_All(A);
  printf("\n");

  Nrow = 6;
  sparseInit(&A, Nrow); /* the sparse simultaneous equation, 6 */
  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(A, i, Nrow, 8.0+4.0*i);
    }
  }
  printf("initial set of equations in A \n");
  sparseWrite_All(A);
  sparseSimeq(A,X);
  PrintV(X,Nrow);
  printf("\n");

  sparseWrite_All(A);
  printf("\n");

  printf("last simeq check \n");
  Nrow = 7;
  sparseInit(&A, Nrow);
  for(i=0; i<Nrow; i++)
  {
    for(j=0; j<Nrow; j++)
      sparsePut(A,i,j,1.0/(double)(i+j+1));
    sparsePut(A,i,Nrow,(double)(i+1)); /* RHS */
  }
  sparseGetRHS(A,Yc); /* for checking */
  sparseSimeq(A,X); /* A * X = Y */

  for(i=0; i<Nrow; i++) /* restore A */
    for(j=0; j<Nrow; j++)
      sparsePut(A,i,j,1.0/(double)(i+j+1));
  sparseMul(A,X,Y); /* Y = A * X check against Yc */

  printf("expected result \n");
  PrintV(Yc,Nrow);
  printf("computed result \n");
  PrintV(Y,Nrow);

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