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

static int log2(int probsize) {
  int retval;

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

  return retval;
}


static void CreateRandomArray(SparseArray X,int probsize) {
  int i;
  int j;
  double density;
  double randval;
  static int num=0;

  if (num == 0) {
    srandom(1);
  } else {
    srandom(2);
  }
  density = ((double)log2(probsize))/probsize;
  for (i=0;i<probsize;i++) {
    for (j=0;j<probsize;j++) {
      randval = ((double)random())/RAND_MAX;
      if (randval < density) {
	Store(X,i,j,1.0+(1.0/(probsize*probsize))*(i*probsize + j));
      }
    }
  }
  if (probsize < 200) {
    char *name;

    name = (char *)malloc(32*sizeof(char));
    sprintf(name,"check.%.3d.in.%d.C.mat",probsize,num+1);

    {
      FILE *outfile;

      outfile = fopen(name,"w");
      WriteDense(X,outfile);
      fclose(outfile);
    }

    free(name);
  }
  num = (num+1)%2;
}


static void CreateNotSoRandomArray(SparseArray X,int probsize) {
  int i;
  int j;
  static int num=0;

  for (i=0;i<probsize;i++) {
    for (j=0;j<probsize;j++) {
      if (num%2 == 0) {
	if (i==j) {
	  Store(X,i,j,1.0);
	}
      } else {
	Store(X,i,j,i*probsize+j+1);
      }
    }
  }
  num++;
}


static void MatMult(SparseArray Res,SparseArray X,SparseArray Y,int n) {
  int rowX;
  int colX;
  int rowY;
  int colY;
  double Xval;
  double Yval;
  double Resval;
  double *Resvalptr;

  for (rowX=0; rowX<n; rowX++) {
    colX = SetCursorRowStart(X,rowX);
    while (colX != -1) {
      Xval = Read(X,rowX,colX);

      rowY = colX;
      colY = SetCursorRowStart(Y,rowY);
      while (colY != -1) {
	Yval = Read(Y,rowY,colY);
	Resval = Xval*Yval;

	Resvalptr = Access(Res,rowX,colY);
	if (Resvalptr == NULL) {
	  Store(Res,rowX,colY,Resval);
	} else {
	  (*Resvalptr) += Resval;
	}

	colY = NextInRow(Y);
      }

      colX = NextInRow(X);
    }
  }
}


static void ColMatMult(SparseArray Res,SparseArray X,SparseArray Y,int n) {
  int rowX;
  int colX;
  int rowY;
  int colY;
  double Xval;
  double Yval;
  double Resval;
  double *Resvalptr;

  for (colX=0; colX<n; colX++) {
    rowX = SetCursorColStart(X,colX);
    while (rowX != -1) {
      Xval = Read(X,rowX,colX);

      colY = rowX;
      rowY = SetCursorColStart(Y,colY);
      while (rowY != -1) {
	Yval = Read(Y,rowY,colY);
	Resval = Xval*Yval;

	Resvalptr = Access(Res,rowX,colY);
	if (Resvalptr == NULL) {
	  Store(Res,rowX,colY,Resval);
	} else {
	  *Resvalptr = Resval;
	}

	rowY = NextInCol(Y);
      }

      rowX = NextInCol(X);
    }
  }
}


static double RunTiming(int probsize) {
  SparseArray A;
  SparseArray B;
  SparseArray C;
  double setuptime;
  double runtime;
  double multtime;

  A = NewSparseArray(probsize,probsize,0.0);
  B = NewSparseArray(probsize,probsize,0.0);
  C = NewSparseArray(probsize,probsize,0.0);
  ResetTimer();

  CreateRandomArray(A,probsize);
  CreateRandomArray(B,probsize);

  setuptime = CheckTimer();

  MatMult(C,A,B,probsize);

  runtime = CheckTimer();
  multtime = runtime-setuptime;

  printf("%9d : %7.3f (%7.3f + %7.3f) seconds\n",probsize,runtime,setuptime,
	 multtime);

  if (probsize < 200) {
    char *name;

    name = (char *)malloc(32*sizeof(char));
    sprintf(name,"check.%.3d.out.C.mat",probsize);

    {
      FILE *outfile;

      outfile = fopen(name,"w");
      WriteDense(C,outfile);
      fclose(outfile);
    }

    free(name);
  }

  DeleteSparseArray(A);
  DeleteSparseArray(B);
  DeleteSparseArray(C);

  return runtime;
}


static void RunColTiming(int probsize) {
  SparseArray A;
  SparseArray B;
  SparseArray C;
  double setuptime;
  double runtime;
  double multtime;

  A = NewSparseArray(probsize,probsize,0.0);
  B = NewSparseArray(probsize,probsize,0.0);
  C = NewSparseArray(probsize,probsize,0.0);

  ResetTimer();

  CreateRandomArray(A,probsize);
  CreateRandomArray(B,probsize);

  setuptime = CheckTimer();

  MatMult(C,A,B,probsize);

  runtime = CheckTimer();
  multtime = runtime-setuptime;

  printf("(column)  : %7.3f (%7.3f + %7.3f) seconds\n",runtime,setuptime,
	 multtime);

  if (probsize < 200) {
    char *name;

    name = (char *)malloc(32*sizeof(char));
    sprintf(name,"check.%.3d.out.col.C.mat",probsize);

    {
      FILE *outfile;

      outfile = fopen(name,"w");
      WriteDense(C,outfile);
      fclose(outfile);
    }

    free(name);
  }

  DeleteSparseArray(A);
  DeleteSparseArray(B);
  DeleteSparseArray(C);
}


void main() {
  int probsize;
  double runtime;

  probsize = 10;

  printf(" probsize   runtime (  setup + multiply )\n");
  printf("------------------------------------------\n");
  do {
    runtime = RunTiming(probsize);
    RunColTiming(probsize);
    probsize *= 2;
  } while (runtime < 120.0);
}