#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define DIMENSION 7

/*******************************************************************/

void print_matrix(double **A, int d) {
    int i, j;
    for(i=0; i<d ; i++) {
       for(j=0; j<d; j++) {
          printf("% 6.2e  ", A[i][j]);
       }
       printf("\n");
    }
    printf("\n");
}

/*******************************************************************/

void print_vector(double *c, int d) {
   int i;
   for(i = 0; i < d ; i++) {
      printf("% 6.2e  ", c[i]);
   }
   printf("\n");
}

/*******************************************************************/

int main() {
   double **A , *b, *c, **Coef;
   int   i, j, k, l, m, d, nn, rmax, **JCoef;
   double  *aa;    /* storage in coordinate form */
   int     *jr, *jc, *ja, *ia; 
   clock_t t1, t2;
   int   nD, D, *aD, *num, num_tmp;
   double **Diag, *tmp, *A_tmp, b_tmp;
   double *dj;
   int    *jdiag, *idiag, ndiag;

   int z;

   d = DIMENSION;

   A = (double **)malloc( d * sizeof(double*) );
   for(i=0; i<d; i++) 
      A[i] = (double *)malloc( d * sizeof(double) );

   b = (double *)malloc( d * sizeof(double) );
   c = (double *)malloc( d * sizeof(double) );

   srand( (unsigned)time((long *)0) );
   for(i = 0; i < d ; i++) 
      for(j = 0; j < d ; j++) 
         if( rand()/2147483647.0 > 0.6 ) {
            A[i][j] = rand() /2147483647.0;
         }
   for(i = 0; i < d ; i++) {
      b[i] = rand() /2147483647.0;
      c[i] = 0.0;
   }

/*------ full matrix multiplication -------------------------------*/

   t1 = clock();
   for(i=0; i < d ; i++)
      for(j=0; j < d; j++)      
         c[i] += A[i][j]*b[j];
   t2 = clock(); 
   printf("Full matrix multiplication:  %6.2e\n", (double)(t2-t1));

/*------ count the non zero elements ------------------------------*/

   nn = 0;
   for(i = 0; i < d ; i++)
      for(j = 0; j < d ; j++)
         if( A[i][j] != 0.0 ) nn++;

/*------ storage in coordinate form -------------------------------*/
    
   aa = (double *)malloc( nn * sizeof(double) );
   jr = (int *)malloc( nn * sizeof(int) );
   jc = (int *)malloc( nn * sizeof(int) );

   k = 0; 
   for(i = 0; i < d ; i++) {
      for(j = 0; j < d ; j++) {
         if( A[i][j] != 0.0 ) {
            aa[k] = A[i][j];
            jr[k] = i;
            jc[k] = j;
            k++;
         }
      }
   }
   for(i = 0; i < d ; i++)  c[i] = 0.0;
   t1 = clock();
   for(k = 0; k < nn; k++)
      c[jr[k]] +=  aa[k] * b[jc[k]];
   t2 = clock();
   printf("Storage in coordinte form:  %6.2e\n", (double)(t2-t1));

   free(aa);  free(jr);  free (jc); 

/*------- compressed sparse row format ----------------------------*/

   aa = (double *)malloc( (nn+1) * sizeof(double) );
   ja = (int *)malloc(    (nn+1) * sizeof(int) );
   ia = (int *)malloc( (d+1) * sizeof(int) );

   k = 0;
   for(i = 0; i < d ; i++) {
      ia[i] = k;
      for(j = 0; j < d ; j++) {
         if( A[i][j] != 0.0 ) {
            aa[k] = A[i][j];
            ja[k] = j;
            k++;
         }
      }
   }
   ia[d] = nn;

   for(i = 0; i < d ; i++)  c[i] = 0.0;
   t1 = clock();
   for(i = 0; i < d; i++) 
      for(j = ia[i]; j < ia[i+1]; j++)
         c[i] += aa[j] * b[ja[j]];
   t2 = clock();
   printf("Compressed Sparse Row Format:  %6.2e\n", (double)(t2-t1));
   free(aa);  free(ja);  free(ia);

/*------- compressed sparse column format -------------------------*/

   aa = (double *)malloc( (nn+1) * sizeof(double) );
   ia = (int *)malloc(    (nn+1) * sizeof(int) );
   ja = (int *)malloc( (d+1) * sizeof(int) );

   k = 0;
   for(j = 0; j < d ; j++) {
      ja[j] = k;
      for(i = 0; i < d ; i++) {
         if( A[i][j] != 0.0 ) {
            aa[k] = A[i][j];
            ia[k] = i;
            k++;
         }
      }
   }
   ja[d] = nn;

   for(i = 0; i < d ; i++)  c[i] = 0.0;
   t1 = clock();
   for(j = 0; j < d; j++)
      for(i = ja[j]; i < ja[j+1]; i++)
         c[ia[i]] += aa[i] * b[j];
   t2 = clock();
   printf("Compressed Sparse Column Format:  %6.2e\n", (double)(t2-t1));
   free(aa);  free(ja);  free(ia);

/*------- CSR with extraction of the diagonal entries  ------------*/

/* count the non diagonal non zero elements */
   nn = 0;
   for(i = 0; i < d ; i++)
      for(j = 0; j < d ; j++)
         if( (i != j) && (A[i][j] != 0.0 ) ) nn++;

   aa = (double *)malloc( (d+nn+1) * sizeof(double) );
   ja = (int *)malloc(    (d+nn+1) * sizeof(int) );

   k = d+1;
   for(i = 0; i < d ; i++) {
      aa[i] = A[i][i];
      ja[i] = k;
      for(j = 0; j < d ; j++) {
         if( (A[i][j] != 0.0) && (i != j) ) {
            aa[k] = A[i][j];
            ja[k] = j;
            k++;
         }
      }
   }
   aa[d] = 0;
   ja[d] = d + nn +1;

   for(i = 0; i < d ; i++)  c[i] = 0.0;
   t1 = clock();
   for(i = 0; i < d; i++) {
      c[i] = aa[i] * b[i]; 
      for(j = ja[i]; j < ja[i+1]; j++) {
         c[i] += aa[j] * b[ja[j]];
      }
   }
   t2 = clock();
   printf("CSR with extraction of the diagonal etries:  %6.2e\n", (double)(t2-t1));

   free(aa);  free(ja);


/*------- diagonalwise storage ------------------------------------*/

/* find number of diagonals (nD) and diagonal description (D): */
   nD = 0;   
   aD = (int *)malloc( (2*d-1) * sizeof(int) );

   for( D = -(d-1); D <= (d-1); D++) {
      for( j = D; j <= (d-1)+D; j++) {                     /* in our notation: j-i = D */
         i = j - D;
         if( (i >= 0) && (i <= d-1) && (j >= 0) && (j <= d-1) && (A[i][j] != 0) ) {
            aD[nD] = D;
            nD++;
            break;
         }
      }
   }
   Diag = (double **)malloc( d * sizeof(double*) );
   for( i = 0; i < d; i++) 
      Diag[i]= (double *)malloc( nD * sizeof(double) );

   for( i = 0; i < d; i++) 
      for( j = 0; j < nD; j++) 
         Diag[i][j] = 0./0.;   /* Initialization with "Not a Number" */

   for( k = 0; k < nD; k++) {
      for( i = 0; i < d; i++) {
         j = i + aD[k];
         if( (j >= 0) && (j <= d-1) )
            Diag[i][k] =  A[i][j];
      }
   }

   for(i = 0; i < d ; i++)  c[i] = 0.0;
   t1 = clock();
   for(k = 0; k < nD; k++) {
      if(aD[k] < 0) {
         for(i = -aD[k]; i <= (d-1); i++) 
            c[i] += Diag[i][k] * b[i+aD[k]];
      } else {
         for(i = 0; i <= (d-1)-aD[k]; i++) 
            c[i] += Diag[i][k] * b[i+aD[k]];
      } 
   }
   t2 = clock();
   printf("Diagonalwise storage:  %6.2e\n", (double)(t2-t1));

   free( aD ); 
   for( i = 0; i < d; i++)
      free( Diag[i] );
   free( Diag );

/*------- rectangular, rowwise storage scheme ---------------------*/

/* find the longest row of A with "rmax" non-zero elements */
   rmax = 0;
   for( i = 0; i < d; i++) {
      k = 0;  
      for( j = 0; j < d; j++)
         if( A[i][j] != 0) { k++; }
      if(k > rmax) { rmax = k; }
   }
   Coef = (double **)malloc( d * sizeof(double*) );
   JCoef = (int **)malloc( d * sizeof(int*) );
   for( i = 0; i < d; i++) {
      Coef[i] = (double *)malloc( rmax * sizeof(double) );
      JCoef[i]= (int *)malloc( rmax * sizeof(int) );
   }

   for( i = 0; i < d; i++) {
      for( j = 0; j < rmax; j++) {
         Coef[i][j]  = 0.0;            /* Initialization */
         JCoef[i][j] = -1; 
      }
   }
   for( i = 0; i < d; i++) {
      k = 0;
      for( j = 0; j < d; j++) {
         if( A[i][j] != 0.0 ) {
            Coef[i][k]  = A[i][j];
            JCoef[i][k] = j;
            k++;
         }
      }
   }

   for(i = 0; i < d ; i++)  c[i] = 0.0;
   t1 = clock();
   for( i = 0; i < d; i++) {
       for(k = 0; k < rmax; k++) {
          if( JCoef[i][k] != -1 ) {
             c[i] += Coef[i][k] * b[JCoef[i][k]];
          } else {
             break;
          }
      }
   }
   t2 = clock();
   printf("Rectangular, rowwise storage:  %6.2e\n", (double)(t2-t1));

   for( i = 0; i < d; i++) {
      free( Coef[i] );
      free( JCoef[i] );
   }
   free( Coef );  free( JCoef );

/*------- jagged diagonal form ------------------------------------*/
/* not yet implemented */

/*-----------------------------------------------------------------*/


}
/*******************************************************************/