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

#define ARRAY_SIZE 10

static void SetupTridiag(SparseArray<double> & X) {
  int row;
  int col;
  int diff;

  for (row=0; row<ARRAY_SIZE; row++) {
    for (col=0; col<ARRAY_SIZE; col++) {
      diff = row-col;
      if (diff >= -1 && diff <= 1) {
	X.Store(row,col,row+(0.01*col));
      }
    }
  }
}


static void SetupReverseTridiag(SparseArray<double> & X) {
  int row;
  int col;
  int diff;

  for (row=0; row<ARRAY_SIZE; row++) {
    for (col=0; col<ARRAY_SIZE; col++) {
      diff = row - (ARRAY_SIZE-(col+1));
      if (diff >= -1 && diff <= 1) {
	X.Store(row,col,row+(0.01*col));
      }
    }
  }
}


static int log2(int probsize) {
  int retval;

  retval = -1;
  while (probsize) {
    probsize = probsize >> 1;
    retval++;
  }

  return retval;
}


static void SetupRandom(SparseArray<double> & X) {
  int i;
  int j;
  double density;
  double randval;

  srand(1);
  density = ((double)log2(ARRAY_SIZE))/ARRAY_SIZE;
  for (i=0; i<ARRAY_SIZE; i++) {
    for (j=0; j<ARRAY_SIZE; j++) {
      randval = ((double)rand())/RAND_MAX;
      if (randval < density) {
	X.Store(i,j,1.0+(1.0/(ARRAY_SIZE*ARRAY_SIZE))*(i*ARRAY_SIZE + j));
      }
    }
  }
}


static void Add(SparseArray<double> & Sum,SparseArray<double> & X,SparseArray<double> & Y) {
  int xrow,xcol;
  int yrow,ycol;
  int rownorm=1;
  int colnorm=1;

  X.SetupIteration(1,2);
  Y.SetupIteration(1,2);

  X.GetCursorPosition(xrow,xcol);
  Y.GetCursorPosition(yrow,ycol);

  do {
    if (xrow == yrow && xcol == ycol) {
      Sum.Store(xrow,xcol,(X.Read(xrow,xcol)+Y.Read(xrow,xcol)));

      X.Advance();
      Y.Advance();

      X.GetCursorPosition(xrow,xcol);
      Y.GetCursorPosition(yrow,ycol);
    } else if (xrow < yrow || (xrow == yrow && xcol < ycol)) {
      Sum.Store(xrow,xcol,X.Read(xrow,xcol));
      X.Advance();
      X.GetCursorPosition(xrow,xcol);
    } else {
      Sum.Store(yrow,ycol,Y.Read(yrow,ycol));
      Y.Advance();
      Y.GetCursorPosition(yrow,ycol);
    }
    if (xrow == -1) {
      xrow = ARRAY_SIZE;
      xcol = ARRAY_SIZE;
    } 
    if (yrow == -1) {
      yrow = ARRAY_SIZE;
      ycol = ARRAY_SIZE;
    }
  } while (xrow != ARRAY_SIZE || yrow != ARRAY_SIZE);
}


static void Double(SparseArray<double> & Res,SparseArray<double> & X,int indim,int outdim) {
  int xrow,xcol;

  X.SetupIteration(indim,outdim);
  X.GetCursorPosition(xrow,xcol);
  do {
    Res.Store(xrow,xcol,2*X.Read(xrow,xcol));
    X.Advance();
    X.GetCursorPosition(xrow,xcol);
  } while (xrow != -1);
}


void main() {
  SparseArray<double> A(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray<double> B(ARRAY_SIZE,ARRAY_SIZE,0.0);
  SparseArray<double> C(ARRAY_SIZE,ARRAY_SIZE,0.0);

  SetupTridiag(A);
  SetupReverseTridiag(B);
  SetupRandom(C);

  cout.precision(2);
  cout.setf(ios::fixed);

  cout << "A is:" << endl;
  A.WriteDense(cout);
  cout << endl << endl;

  cout << "B is:" << endl;
  B.WriteDense(cout);
  cout << endl << endl;

  cout << "C is:" << endl;
  C.WriteDense(cout);
  cout << endl << endl;

  {
    SparseArray<double> Sum(ARRAY_SIZE,ARRAY_SIZE,0.0);
    cout << "A + A = " << endl;
    Double(Sum,A,-1,2);
    Sum.WriteDense(cout);
    cout << endl << endl;
  }

  {
    SparseArray<double> Sum(ARRAY_SIZE,ARRAY_SIZE,0.0);
    cout << "A + B = " << endl;
    Add(Sum,A,B);
    Sum.WriteDense(cout);
    cout << endl << endl;
  }

  {
    SparseArray<double> Sum(ARRAY_SIZE,ARRAY_SIZE,0.0);
    cout << "A + C = " << endl;
    Add(Sum,A,C);
    Sum.WriteDense(cout);
    cout << endl << endl;
  }

  {
    SparseArray<double> Sum(ARRAY_SIZE,ARRAY_SIZE,0.0);
    cout << "B + B = " << endl;
    Double(Sum,B,-2,1);
    Sum.WriteDense(cout);
    cout << endl << endl;
  }

  {
    SparseArray<double> Sum(ARRAY_SIZE,ARRAY_SIZE,0.0);
    cout << "B + C = " << endl;
    Add(Sum,B,C);
    Sum.WriteDense(cout);
    cout << endl << endl;
  }

  {
    SparseArray<double> Sum(ARRAY_SIZE,ARRAY_SIZE,0.0);
    cout << "C + C = " << endl;
    Double(Sum,C,-2,-1);
    Sum.WriteDense(cout);
    cout << endl << endl;
  }

  {
    SparseArray<double> PartSum(ARRAY_SIZE,ARRAY_SIZE,0.0);
    SparseArray<double> Sum(ARRAY_SIZE,ARRAY_SIZE,0.0);
    cout << "A + B + C = " << endl;
    Add(PartSum,A,B);
    Add(Sum,PartSum,C);
    Sum.WriteDense(cout);
    cout << endl << endl;
  }
}