#include <stdio.h>
#include "sparse.h"

#define ARRAY_SIZE 9

static int dumperrors = 0;

static double tempval;
static int phase = 0;

#define correctval(r,c) (tempval=CORRECTVAL(r,c),tempval)
#define negval(r,c) (tempval=NEGVAL(r,c),tempval)

static double CORRECTVAL(int r,int c) {
  return ((r+1) + (0.1*(c+1)));
}

static double NEGVAL(int r,int c) {
  return -((r+1) + (0.1*(c+1)));
}

static void SetupTridiag(SparseArray& Tridiag) {
  int row;
  int col;
  int diff;

  /* Set up a tridiagonal using a doubly-nested loop */
  for (row=0; row<ARRAY_SIZE; row++) {
    for (col=0; col<ARRAY_SIZE; col++) {
      diff = row-col;
      if (diff >= -1 && diff <= 1) {
	Tridiag.Store(row,col,correctval(row,col));
      }
    }
  }
}


static void SetupTridiag2(SparseArray& Tridiag) {
  int rowcol;

  /* Set up a tridiagonal using a singly-nested loop */
  Tridiag.Store(ARRAY_SIZE-1,ARRAY_SIZE-1,correctval(ARRAY_SIZE-1,ARRAY_SIZE-1));
  for (rowcol=0; rowcol<ARRAY_SIZE-1; rowcol++) {
    Tridiag.Store(rowcol+1,rowcol,correctval(rowcol+1,rowcol));
    Tridiag.Store(rowcol,rowcol+1,correctval(rowcol,rowcol+1));
    Tridiag.Store(rowcol,rowcol,correctval(rowcol,rowcol));
  }
}


void Dump(SparseArray &X,char *subname) {
  char name[80];
  
  if (phase==0) {
    sprintf(name,"dump.normal.%s.out",subname);
  } else {
    sprintf(name,"dump.unusual.%s.out",subname);
  }
  ofstream outfile(name,ios::out);

  X.WriteDense(outfile);
}


void TestRead(SparseArray& X) {
  int row;
  int col;
  int diff;
  double val;
  int problem;

  problem = 0;

  /* Make sure that values in the tridiagonal are correct using a
     reverse doubly-nested loop */
  for (col=ARRAY_SIZE-1; col>=0; col--) {
    for (row=ARRAY_SIZE-1; row>=0; row--) {
      val = X.Read(row,col);
      diff = row-col;
      if (diff >= -1 && diff <= 1) {
	problem |= (val != correctval(row,col));
      } else {
	problem |= (val != 0.0);
      }
    }
  }

  if (problem) {
    fprintf(stderr,"ERROR: TestRead() failed\n");
    if (dumperrors) {
      Dump(X,"Read");
    }
  } else {
    fprintf(stderr,"       TestRead() passed\n");
  }
}


void TestStore(SparseArray& X) {
  int rowcol;
  int row;
  int col;
  int diff;
  int problem;
  double val;

  problem = 0;

  /* Negate the main diagonal using a reverse single loop */
  for (rowcol=ARRAY_SIZE-1; rowcol>=0; rowcol--) {
    X.Store(rowcol,rowcol,negval(rowcol,rowcol));
  }

  /* Check that the main diagonal is properly negated using a
     doubly-nested loop */
  for (row=0; row<ARRAY_SIZE; row++) {
    for (col=0; col<ARRAY_SIZE; col++) {
      val = X.Read(row,col);
      diff = row-col;
      if (diff == -1 || diff == 1) {
	problem |= (val != correctval(row,col));
      } else if (diff == 0) {
	problem |= (val != negval(row,col));
      } else {
	problem |= (val != 0.0);
      }
    }
  }

  if (problem) {
    if (dumperrors) {
      Dump(X,"Store-halfway");
    }
  }

  /* Fix the main diagonal using a forwards loop */
  for (rowcol=0; rowcol<ARRAY_SIZE; rowcol++) {
    X.Store(rowcol,rowcol,correctval(rowcol,rowcol));
  }

  /* Check that it's correct using a reverse doubly-nested loop */
  for (col=ARRAY_SIZE-1; col>=0; col--) {
    for (row=ARRAY_SIZE-1; row>=0; row--) {
      val = X.Read(row,col);
      diff = row-col;
      if (diff >= -1 && diff <= 1) {
	problem |= (val != correctval(row,col));
      } else {
	problem |= (val != 0.0);
      }
    }
  }

  if (problem) {
    fprintf(stderr,"ERROR: TestStore() failed\n");
    if (dumperrors) {
      Dump(X,"Store");
    }
  } else {
    fprintf(stderr,"       TestStore() passed\n");
  }
}


void TestAccess(SparseArray& X) {
  int rowcol;
  int row;
  int col;
  int diff;
  int problem;
  double val;
  double* valptr;

  problem = 0;

  /* Negate the main diagonal using a reverse single loop */
  for (rowcol=ARRAY_SIZE-1; rowcol>=0; rowcol--) {
    valptr = X.Access(rowcol,rowcol);
    if (valptr == NULL) {
      problem = 1;
    } else {
      problem |= (*valptr != correctval(rowcol,rowcol));
      *valptr = negval(rowcol,rowcol);
    }
  }

  if (problem) {
    fprintf(stderr,"ERROR: simple TestAccess() failed\n");
    if (dumperrors) {
      Dump(X,"Access-simple");
    }
  } else {
    fprintf(stderr,"       simple TestAccess() passed\n");
  }

  /* Check that the main diagonal is properly negated using a
     doubly-nested loop */
  for (row=0; row<ARRAY_SIZE; row++) {
    for (col=0; col<ARRAY_SIZE; col++) {
      val = X.Read(row,col);
      valptr = X.Access(row,col);
      diff = row-col;
      if (diff == -1 || diff == 1) {
	problem |= (val != correctval(row,col));
	problem |= (valptr == NULL);
      } else if (diff == 0) {
	problem |= (val != negval(row,col));
	problem |= (valptr == NULL);
      } else {
	problem |= (val != 0.0);
	problem |= (valptr != NULL);
      }
    }
  }

  if (problem) {
    if (dumperrors) {
      Dump(X,"Access-halfway");
    }
  }

  /* Fix the main diagonal using a forwards loop */
  for (rowcol=0; rowcol<ARRAY_SIZE; rowcol++) {
    valptr = X.Access(rowcol,rowcol);
    if (valptr == NULL) {
      problem = 1;
    } else {
      problem |= (*valptr != negval(rowcol,rowcol));
      *valptr = correctval(rowcol,rowcol);
    }
  }

  /* Check that it's correct using a reverse doubly-nested loop */
  for (col=ARRAY_SIZE-1; col>=0; col--) {
    for (row=ARRAY_SIZE-1; row>=0; row--) {
      val = X.Read(row,col);
      valptr = X.Access(row,col);
      diff = row-col;
      if (diff >= -1 && diff <= 1) {
	problem |= (val != correctval(row,col));
	problem |= (valptr == NULL);
      } else {
	problem |= (val != 0.0);
	problem |= (valptr != NULL);
      }
    }
  }

  if (problem) {
    fprintf(stderr,"ERROR: TestAccess() failed\n");
    if (dumperrors) {
      Dump(X,"Access");
    }
  } else {
    fprintf(stderr,"       TestAccess() passed\n");
  }
}


void TestURV(SparseArray &X) {
  int row;
  int col;
  int diff;  
  double val;
  int problem;

  problem = 0;

  val = X.GetURV();
  problem |= (val != 0.0);
  X.SetURV(-1.0);
  val = X.GetURV();
  problem |= (val != -1.0);
  
  /* Check that it's correct using a reverse doubly-nested loop */
  for (col=ARRAY_SIZE-1; col>=0; col--) {
    for (row=ARRAY_SIZE-1; row>=0; row--) {
      val = X.Read(row,col);
      diff = row-col;
      if (diff >= -1 && diff <= 1) {
	problem |= (val != correctval(row,col));
      } else {
	problem |= (val != -1.0);
      }
    }
  }

  if (problem) {
    if (dumperrors) {
      Dump(X,"URV-halfway");
    }
  }

  X.SetURV(0.0);
  val = X.GetURV();
  problem |= (val != 0.0);

  /* Check that it's correct using a reverse doubly-nested loop */
  for (col=0; col<ARRAY_SIZE; col++) {
    for (row=0; row<ARRAY_SIZE; row++) {
      val = X.Read(row,col);
      diff = row-col;
      if (diff >= -1 && diff <= 1) {
	problem |= (val != correctval(row,col));
      } else {
	problem |= (val != 0.0);
      }
    }
  }

  if (problem) {
    fprintf(stderr,"ERROR: TestURV() failed\n");
    if (dumperrors) {
      Dump(X,"URV");
    }
  } else {
    fprintf(stderr,"       TestURV() passed\n");
  }
}


void TestDenseIO(SparseArray &X,char *name) {
  int row;
  int col;
  double val1;
  double val2;
  int problem;
  
  problem = 0;
  {
    ofstream outfile(name,ios::out);
    X.WriteDense(outfile);
  }
  {
    ifstream infile(name,ios::in|ios::nocreate);
    SparseArray Y(infile,ARRAY_SIZE,ARRAY_SIZE,0.0);

    for (col=ARRAY_SIZE-1; col>=0; col--) {
      for (row=ARRAY_SIZE-1; row>=0; row--) {
	val1 = X.Read(row,col);
	val2 = Y.Read(row,col);
	problem |= (val1 != val2);
      }
    }
    if (problem) {
      fprintf(stderr,"ERROR: TestDenseIO() failed\n");
      if (dumperrors) {
	Dump(Y,"Dense");
      }
    } else {
      fprintf(stderr,"       TestDenseIO() passed\n");
    }
  }
}


void TestSparseIO(SparseArray &X,char *name) {
  int row;
  int col;
  double val1;
  double val2;
  int problem;

  problem=0;
  {
    ofstream outfile(name,ios::out);
    X.WriteSparse(outfile);
  }
  {
    ifstream infile(name,ios::in|ios::nocreate);
    SparseArray Y(infile);

    for (col=ARRAY_SIZE-1; col>=0; col--) {
      for (row=ARRAY_SIZE-1; row>=0; row--) {
	val1 = X.Read(row,col);
	val2 = Y.Read(row,col);
	problem |= (val1 != val2);
      }
    }
    if (problem) {
      fprintf(stderr,"ERROR: TestSparseIO() failed\n");
      if (dumperrors) {
	Dump(Y,"Sparse");
      }
    } else {
      fprintf(stderr,"       TestSparseIO() passed\n");
    }
  }
}


void TestIterateRight(SparseArray& X) {
  int row;
  int col;
  int problem;
  int i;
  int lo;
  int hi;

  problem = 0;

  for (row=ARRAY_SIZE-1; row>=0; row--) {
    col = X.SetCursorRowStart(row);
    lo = -1;
    hi = 1;
    if (row == 0) {
      lo = 0;
    } else if (row == ARRAY_SIZE-1) {
      hi = 0;
    }
    for (i=lo;i<=hi;i++) {
      problem |= (col != row+i);
      col = X.NextInRow();
    }
    problem |= (col != -1);
  }

  if (problem) {
    fprintf(stderr,"ERROR: TestIterateRight() failed\n");
  } else {
    fprintf(stderr,"       TestIterateRight() passed\n");
  }
}


void TestIterateLeft(SparseArray& X) {
  int row;
  int col;
  int problem;
  int i;
  int lo;
  int hi;

  problem = 0;

  for (row=ARRAY_SIZE-1; row>=0; row--) {
    col = X.SetCursorRowEnd(row);
    lo = -1;
    hi = 1;
    if (row == 0) {
      lo = 0;
    } else if (row == ARRAY_SIZE-1) {
      hi = 0;
    }
    for (i=hi;i>=lo;i--) {
      problem |= (col != row+i);
      col = X.PrevInRow();
    }
    problem |= (col != -1);
  }

  if (problem) {
    fprintf(stderr,"ERROR: TestIterateLeft() failed\n");
  } else {
    fprintf(stderr,"       TestIterateLeft() passed\n");
  }
}


void TestIterateDown(SparseArray& X) {
  int row;
  int col;
  int problem;
  int i;
  int lo;
  int hi;

  problem = 0;

  for (col=ARRAY_SIZE-1; col>=0; col--) {
    row = X.SetCursorColStart(col);
    lo = -1;
    hi = 1;
    if (col == 0) {
      lo = 0;
    } else if (col == ARRAY_SIZE-1) {
      hi = 0;
    }
    for (i=lo;i<=hi;i++) {
      problem |= (row != col+i);
      row = X.NextInCol();
    }
    problem |= (row != -1);
  }

  if (problem) {
    fprintf(stderr,"ERROR: TestIterateDown() failed\n");
  } else {
    fprintf(stderr,"       TestIterateDown() passed\n");
  }
}


void TestIterateUp(SparseArray& X) {
  int row;
  int col;
  int problem;
  int i;
  int lo;
  int hi;

  problem = 0;

  for (col=ARRAY_SIZE-1; col>=0; col--) {
    row = X.SetCursorColEnd(col);
    lo = -1;
    hi = 1;
    if (col == 0) {
      lo = 0;
    } else if (col == ARRAY_SIZE-1) {
      hi = 0;
    }
    for (i=hi;i>=lo;i--) {
      problem |= (row != col+i);
      row = X.PrevInCol();
    }
    problem |= (row != -1);
  }

  if (problem) {
    fprintf(stderr,"ERROR: TestIterateUp() failed\n");
  } else {
    fprintf(stderr,"       TestIterateUp() passed\n");
  }
}


void IterateAround(SparseArray& X) {
  int rowcol;
  int problem;
  int row;
  int col;

  problem = 0;

  row = 0;
  col = X.SetCursorRowStart(row);
  for (rowcol=0; rowcol<ARRAY_SIZE-1; rowcol++) {
    problem |= (row != col);
    col = X.NextInRow();
    problem |= (col != row+1);
    row = X.NextInCol();
  }
  problem |= (row != ARRAY_SIZE-1);
  row = X.NextInCol();
  problem |= (row != -1);

  col = ARRAY_SIZE-1;
  row = X.SetCursorColEnd(col);
  for (rowcol=ARRAY_SIZE-1; rowcol>0; rowcol--) {
    problem |= (row != col);
    col = X.PrevInRow();
    problem |= (col != row-1);
    row = X.PrevInCol();
  }
  problem |= (col != 0);
  col = X.PrevInCol();
  problem |= (col != -1);

  if (problem) {
    fprintf(stderr,"ERROR: IterateAround() failed\n");
  } else {
    fprintf(stderr,"       IterateAround() passed\n");
  }
}


void TestDensity(SparseArray& X) {
  double val;
  int problem;

  problem = 0;

  val = X.Density();
  if (problem) {
    fprintf(stderr,"ERROR: simple TestDensity() failed: %f\n",val);
  } else {
    fprintf(stderr,"       simple TestDensity() passed\n");
  }
  problem |= (val != (double)(2*2 + 3*(ARRAY_SIZE-2))/(ARRAY_SIZE*ARRAY_SIZE));
  X.Store(0,0,-1.0);
  X.Store(0,ARRAY_SIZE-1,-2.0);
  X.Store(ARRAY_SIZE-1,0,-3.0);
  X.Store(ARRAY_SIZE-1,ARRAY_SIZE-1,-4.0);
  val = X.Density();
  problem |= (val != (double)(3*ARRAY_SIZE)/(ARRAY_SIZE*ARRAY_SIZE));

  if (problem) {
    fprintf(stderr,"ERROR: TestDensity() failed: %f\n",val);
    if (dumperrors) {
      Dump(X,"Density");
    }
  } else {
    fprintf(stderr,"       TestDensity() passed\n");
  }
}


void TestNumRepresented(SparseArray& X) {
  int val;
  int problem;

  problem = 0;

  val = X.NumRepresented();
  if (problem) {
    fprintf(stderr,"ERROR: simple TestNumRepresented() failed: %d\n",val);
  } else {
    fprintf(stderr,"       simple TestNumRepresented() passed\n");
  }
  problem |= (val != 2*2 + 3*(ARRAY_SIZE-2));
  X.Store(0,0,-1.0);
  X.Store(0,ARRAY_SIZE-1,-2.0);
  X.Store(ARRAY_SIZE-1,0,-3.0);
  X.Store(ARRAY_SIZE-1,ARRAY_SIZE-1,-4.0);
  val = X.NumRepresented();
  problem |= (val != 3*ARRAY_SIZE);

  if (problem) {
    fprintf(stderr,"ERROR: TestNumRepresented() failed: %d\n",val);
    if (dumperrors) {
      Dump(X,"NumRepresented");
    }
  } else {
    fprintf(stderr,"       TestNumRepresented() passed\n");
  }
}


void main(int argc,char *argv[]) {
  SparseArray T1A(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T1B(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T2A(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T2B(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T3A(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T3B(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T4A(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T4B(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T5A(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T5B(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T6A(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T6B(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T7A(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T7B(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T8A(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T8B(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T9A(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray T9B(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray TaA(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray TaB(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray TbA(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray TbB(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray TcA(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray TcB(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray TdA(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray TdB(ARRAY_SIZE,ARRAY_SIZE,0.0);

  if (argc > 1) {
    dumperrors = 1;
  }

  fprintf(stderr,"Setting up all normal arrays...");
  fflush(stderr);
  SetupTridiag(T1A);
  SetupTridiag(T2A);
  SetupTridiag(T3A);
  SetupTridiag(T4A);
  SetupTridiag(T5A);
  SetupTridiag(T6A);
  SetupTridiag(T7A);
  SetupTridiag(T8A);
  SetupTridiag(T9A);
  SetupTridiag(TaA);
  SetupTridiag(TbA);
  SetupTridiag(TcA);
  SetupTridiag(TdA);
  fprintf(stderr,"done.\n");

  fprintf(stderr,"Testing all operations on normal arrays...\n");
  TestRead(T1A);
  TestStore(T2A);
  TestAccess(T3A);
  TestURV(T4A);
  TestIterateRight(T7A);
  TestIterateLeft(T8A);
  TestIterateDown(T9A);
  TestIterateUp(TaA);
  IterateAround(TbA);
  TestDensity(TcA);
  TestNumRepresented(TdA);
  TestDenseIO(T5A,"dense1.out");
  TestSparseIO(T6A,"sparse1.out");
  
  phase = 1;

  fprintf(stderr,"Setting up all unusual arrays...");
  fflush(stderr);
  SetupTridiag(T1B);
  SetupTridiag(T2B);
  SetupTridiag(T3B);
  SetupTridiag(T4B);
  SetupTridiag(T5B);
  SetupTridiag(T6B);
  SetupTridiag(T7B);
  SetupTridiag(T8B);
  SetupTridiag(T9B);
  SetupTridiag(TaB);
  SetupTridiag(TbB);
  SetupTridiag(TcB);
  SetupTridiag(TdB);
  fprintf(stderr,"done.\n");

  fprintf(stderr,"Testing all operations on unusual arrays...\n");
  TestRead(T1B);
  TestStore(T2B);
  TestAccess(T3B);
  TestURV(T4B);
  TestIterateRight(T7B);
  TestIterateLeft(T8B);
  TestIterateDown(T9B);
  TestIterateUp(TaB);
  IterateAround(TbB);
  TestDensity(TcB);
  TestNumRepresented(TdB);
  TestDenseIO(T5B,"dense2.out");
  TestSparseIO(T6B,"sparse2.out");

  fprintf(stderr,"Deleting all arrays and exiting...\n");
}