/*
 *  Sample code for the DIS Matrix Stressmark
 *
 * This source code is the completely correct source code based on 
 * the example codes provided by Atlantic Aerospace Division, Titan 
 * Systems Corporation, 2000.
 * 
 * If you just compile and generate the executables from this source 
 * code, this code would be enough. However, if you wish to get a complete 
 * understanding of this stressmark, it is strongly suggested that you
 * read the Benchmark Analysis and Specifications Document Version 1.0
 * before going on since the detailed comments are given in this documents.
 * the comments are not repeated here.
 */

/*
 *  The Sparse Matrix Storage is implemented by Compact Row Storage Scheme
 *  In the code, the data is first generated by randomNonzeroFloat()
 *  the data is first stored in a full-space matrix with size of dim*dim
 *  then the data is transfered to the Compact Row Matrix, 
 *  the data value is kept in *value,
 *  the columns corresponding to the value are stored in *col_ind,
 *  the start element of each row is stored in *row_start.
 */
 
/* 
 * Please note: 
 * the total number of data is numberNonzero +dim
 * among which, NumberNonzero because this is symmetric matrix
 * dim because the diagonal elements
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <assert.h>
#include "DISstressmarkRNG.h"

#define MIN_SEED -2147483647
#define MAX_SEED -1
#define MIN_DIM  1
#define MAX_DIM  32768
#define MAX_ITERATIONS 65536
#define MIN_TOLERANCE 0.000007
#define MAX_TOLERANCE 0.5
#define MIN_NUMBER   -3.4e10/dim
#define MAX_NUMBER 3.4e10/dim
#define EPSI   1.0e-10
#define MIN_DIG_NUMBER 1.0e-10
#define MAX_DIG_NUMBER 3.4e10

/*
 *  External variable, dimension
 */

static int dim;

/*      
 *  matrix * vector     
 */

double  *matrixMulvector(double *value, 
			   int *col_ind, 
			   int *row_start,
			   double  *vector)
{  
  int l, ll;
  double  *out;
  double  sum;
  int tmp_rs, tmp_re;
  
  out = (double *)malloc(dim*sizeof(double));
 
  for (l=0; l<dim; l++){  
     *(out + l) = 0;
     tmp_rs = row_start[l];
    
     if (tmp_rs != -1){
      tmp_re = row_start[l+1];   /*
				  *get the start and ending elements of 
				  *  each row
				 */
      for (ll=tmp_rs; ll<tmp_re; ll++){
	*(out + l) += value[ll]*vector[col_ind[ll]];
      }
    }
  }
  return out;
}


/*
 *    vector1 - vector2
 */

double  *vectorSub(double  *vector1, double  *vector2){

  int l;
  double  *vector;
  vector = (double  *)malloc(dim*sizeof(double ));
  for (l=0; l<dim; l++){
    *(vector + l) = *(vector1 + l) - *(vector2 + l);
  }
  return vector;
}


/*
 * vector1 + vector2
 */

double  *vectorAdd(double  *vector1, double  *vector2){

  int l;
  double  *vector;

  vector = (double  *)malloc(dim*sizeof(double ));

  for (l=0; l<dim; l++){
    *(vector + l) = *(vector1 + l) + *(vector2 + l);
  }
  return vector;
} 

/* 
 * vector1 * vector2
 */

double  vectorMul(double  *vector1, double  *vector2){

  int l;
  double  product;

  product = 0;

  for (l=0; l<dim; l++){
    product += (*(vector1 + l))*(*(vector2 + l));

  }
  return product;
} 

/*
 * /vector/
 */

double  vectorValue(double  *vector){

  double  value;
  int l;

  value = 0;

  for (l=0; l<dim; l++){
    value += (*(vector + l)) * (*(vector + l));
  }

  return (sqrt(value));
}

/*
 * transpose(vector)
 * In fact, we return the original vector here
 */

double  *transpose(double  *vector){

  double  *vect;
  int l;

  vect = (double  *)malloc(dim*sizeof(double ));

  for (l=0; l<dim; l++){
    *(vect+l) = *(vector+l);
  }
  return vect;
}

/*
 * value * <vector>
 */
double  *valueMulvector(double  value, double  *vector){

  int l;
  double  *vect;
    int lll;
    double sum;

  vect = (double  *) malloc(dim * sizeof(double ));

  for (l=0; l<dim; l++){
    *(vect + l) = (*(vector + l)) * value;
  }

  return vect;
}
  
/*
 * generate the data distributed sparsely in matrix
 */

void initMatrix(double  *matrix, int dim, int numberNonzero){
  
  int k, l, ll;
  int i, j;

  int lll;
  double sum;

  for (k=0; k< dim*dim; k++){
    *(matrix + k) = 0;
  }

  for (l=0; l<numberNonzero/2; l++){

    i = randomUInt(1, dim-1);
    j = randomUInt(0, i-1);

    while (*(matrix + i*dim + j) != 0){
      
     i++;
       if (i == dim){
       j++;
       if (j == dim-1){
	 j = 0;
	 i = 1;
       }
       else{
	 i = j+1;
       }
     }
    }
  
    if (*(matrix + i*dim + j) == 0){
      *(matrix + i*dim + j) = (double )randomNonZeroFloat(MIN_NUMBER, 
							  MAX_NUMBER, 
							  EPSI);
      *(matrix + j*dim + i) = *(matrix + i*dim + j);
    }
  }
 
  for (ll=0; ll<dim; ll++){


    
    *(matrix + ll*dim + ll) = (double )randomNonZeroFloat(-MAX_DIG_NUMBER,
							  MAX_DIG_NUMBER, 
							  MIN_DIG_NUMBER);
    
    sum = 0;

    for (lll=0; lll<dim; lll++){
      if (lll != ll){
	sum += *(matrix + lll*dim + ll);
      }
    }
    
    if (*(matrix + ll*dim + ll) < sum ){
      *(matrix + ll*dim + ll) += sum;
   }
  }

  return;
}

/*
 * generate the data value in the vectors
 */
 
void initVector(double *vector, int dim){

  int l;
  
  for (l=0; l<dim; l++){
    *(vector + l) = (double )randomFloat (MIN_NUMBER, MAX_NUMBER);
  }

  return;
}

/*
 * make a vector contains value of zero
 */

void zeroVector(double *vector, int dim){
  int l;
  
  for (l=0; l<dim; l++){
    *(vector + l) = 0;
  }
  return;
}

/*
 * return a vector which is the copy of the vect
 */

double *equalVector(double *vect){

  int l;
  double *vect1;
 
  vect1 = (double *)malloc(dim*sizeof(double ));

  for (l=0; l<dim; l++){
    *(vect1+l) = *(vect+l);
  }
  return vect1;
}



void biConjugateGradient(double *value,
			 int *col_ind,
			 int *row_start,
			 double *vectorB, 
			 double *vectorX,
			 double errorTolerance,
			 int maxIterations,
			 double *actualError,
			 int *actualIteration,
			 int dim)
     /* 
      * in the code, we use a lot of temparary vectors and variables
      * this is just for simple and clear
      * you can optimize these temporary variables and vectors 
      * based on your need
      *
      */
{
  double *vectorR;
  double *vectorP, *matrixAvectorP, *nextVectorR;
  double  error;
  int iteration;
  double  alpha, beta;

  double  *temp0, *temp1,*temp2, *temp3;
  double  *temp4, temp5, temp6, *temp7 ;
  double  *temp8, *temp9, temp10, *temp11;
  double  temp12, *temp13, *temp14;
  double  *temp15, *temp16, *temp17;

  int l;
  int ll;

  alpha = 0;
  beta = 0;

  /*
   * vectorR = vectorB - matrixA*vectorX
   */
  temp0 = matrixMulvector(value,col_ind, row_start, vectorX);
  vectorR = vectorSub(vectorB, temp0);

  /*
   * vectorP = vectorR
   */
  vectorP = equalVector(vectorR);

  /*
   * error = |matrixA * vectorX - vectorB| / |vectorB|
   */
  temp2 = vectorSub(temp0, vectorB); 
  error = vectorValue(temp2)/vectorValue(vectorB);

  free(temp0);
  free (temp2);
  
  iteration = 0;
 
 
  while ((iteration < maxIterations) && (error > errorTolerance)){
   
    /* 
     *   alpha = (transpose(vectorR) * vectorR) /
     *           (transpose(vectorP) * (matrixA * vectorP)
     */
    temp1 = matrixMulvector(value, col_ind, row_start, vectorP);    
    temp3 = transpose(vectorR);       
    temp4 = transpose(vectorP);
    temp5 = vectorMul(temp4, temp1);
    temp6 = vectorMul(temp3, vectorR);
    alpha = temp6/temp5;
 
    /* 
     * nextVectorR = vectorR - alpha*(matrixA * vectorP)
     */

    temp7 = valueMulvector(alpha, temp1); 
    temp8 = vectorSub(vectorR, temp7);
    nextVectorR = equalVector(temp8);

    /* 
     * beta = (transpose(nextVectorR) * nextVectorR) /
     *           (transpose(vectorR) * vectorR)
     */
    temp9 = transpose(nextVectorR);
    temp10 = vectorMul(temp9, nextVectorR);
    temp11 = transpose(vectorR);
    temp12 = vectorMul(temp11, vectorR);
    beta = temp10/temp12;

    /*
     * vectorX = vectorX + alpha * vectorP
     */
    temp13 = valueMulvector(alpha, vectorP);       
    vectorX = vectorAdd(vectorX,temp13);

    /* 
     *vectorP = nextVectorR + beta*vectorP
     */       
    temp14 = valueMulvector(beta, vectorP);   
    temp17 = vectorAdd(nextVectorR, temp14);

    for (ll=0; ll<dim; ll++){
      *(vectorP + ll) = *(temp17 + ll);
    }

    /*
     * vectorR = nextVectorR
     */
    
    for (l=0; l<dim; l++){
    *(vectorR+l) = *(nextVectorR+l);
    }

    /* 
     * error = |matrixA * vectorX - vectorB| / |vectorB|
     */
    temp15 = matrixMulvector(value, col_ind,row_start, vectorX);
    temp16 = vectorSub(temp15,vectorB);
    error = vectorValue(temp16)/vectorValue(vectorB);

    free(temp1);
    free(temp3);
    free(temp4);
    free(temp7);
    free(temp8);
    free(temp9);
    free(temp11);
    free(nextVectorR);
    free(temp13);
    free(temp14);
    free(temp15);
    free(temp16);
    free(temp17);

    iteration++;
  }

  *actualError = error;
  *actualIteration = iteration;

  free(vectorR);
  free(vectorP);

  return;
}
  
/*
 * This is the function to transfer the data from the matrix of dense storage 
 * to Compact Row Storage
 */
void create_CRS(double *matrixA,
		double *value, 
		int *col_ind, 
		int *row_start,
		int dim,
		int numberNonzero)
{

  int i, j, k;
  int cnt;
  double tmp;

  /* 
   *initialize the row_start
   */
     
  for(k=0; k<dim; k++){
    row_start[k] = -1;
  }
  
  /* 
   * make the end of the last row to be numberNonzero + dim.
   */

  row_start[dim] = numberNonzero+dim;
  
  /*
   * initialize the col_ind
   */

  for (k=0; k<numberNonzero+dim; k++){
    col_ind[k] = -1;
  }


  cnt = 0;

  for (i=0;  (cnt<numberNonzero+dim)&&(i<dim); i++){
    for (j=0; (cnt<numberNonzero+dim)&&(j<dim); j++){
      
      tmp = *(matrixA + i*dim + j);

         if (tmp!=0){

	   value[cnt] = tmp;
	   col_ind[cnt] = j;
       
	    if (row_start[i] == -1)
	      row_start[i] = cnt;
	    
	    cnt += 1;
	 }	
    }
  }
  row_start[i] = cnt;

  return;
}



int main()
{
  int seed;
  int numberNonzero;
  int maxIterations;
  float errorTolerance;
  double  actualError;
  int actualIteration;
  
  time_t beginTime;
  time_t endTime;

  double  *matrixA;
  double  *vectorB;
  double  *vectorX;

  double *value;
  int    *col_ind;
  int    *row_start;
  int sum;
  int k;

  fscanf(stdin, "%d %d %d %d %f", 
	 &seed, &dim, &numberNonzero,&maxIterations,&errorTolerance);
  assert((seed > MIN_SEED) && (seed < MAX_SEED));
  assert((dim > MIN_DIM) && (dim < MAX_DIM));
  assert((numberNonzero > dim) && (numberNonzero < dim*dim));
  assert((maxIterations > 0) && (maxIterations < MAX_ITERATIONS));
  assert((errorTolerance > MIN_TOLERANCE) && (errorTolerance < MAX_TOLERANCE));
  
  matrixA = (double  *)malloc(dim*dim*sizeof(double ));
  vectorB = (double *)malloc(dim*sizeof(double));
  vectorX = (double *)malloc(dim*sizeof(double));

  value = (double *)malloc((numberNonzero+dim)*sizeof(double));
  col_ind = (int *)malloc((numberNonzero+dim)*sizeof(int));
  row_start = (int *)malloc((dim+1)*sizeof(int));

  randInit(seed);

  initMatrix(matrixA, dim, numberNonzero);
  
  create_CRS(matrixA, value, col_ind, row_start, dim, numberNonzero);

  initVector(vectorB, dim);
  zeroVector(vectorX, dim);
  printf(" after init\n");

  beginTime = time(NULL);
  
  actualError = 0;
  actualIteration = 0;
   
  biConjugateGradient(value, col_ind, row_start, vectorB, vectorX, errorTolerance,
		      maxIterations,
		      &actualError, &actualIteration, dim);
  


  endTime = time(NULL) - beginTime;
  


  sum = 0;
  for (k=1; k<dim; k++){
    sum += sum + *(vectorX + k);
  }
  
  fprintf(stdout, "sum = %d, actualError = %e, actualIteration = %d\n", sum, actualError, actualIteration);
  fprintf(stdout, "total time = %u sec. \n", (unsigned int)endTime);

  return(0);
    }