/* File: dense.c

   Implementation of the Dense Matrix Operations
*/

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "dense.h"

/* local prototypes */
void DCrashOnNull(void *ptr, char *mesg) ;
dm_ptr NewDenseMatrix(int rows, int columns) ;


/* Generic Error Routine for Pointers */

void DCrashOnNull(void *ptr, char *mesg) {
   if (ptr == NULL) {
      fprintf(stderr, "Crash: %s\n", mesg) ;
      exit(1) ;
   }
}


/* Initialize a New Matrix */

dm_ptr NewDenseMatrix(int rows, int columns) {
   dm_ptr M ;
   double **temp2D, *temp1D ;
   int i, j ;

   M = (dm_ptr) malloc(sizeof(dense_matrix)) ;
   DCrashOnNull(M, "Could not make new matrix") ;

   temp2D = (double **) malloc(rows * sizeof(double *)) ;
   DCrashOnNull(temp2D, "Could not make new matrix: create row pointers") ;

   M->rows = rows ;
   M->columns = columns ;
   M->matrix = temp2D ;

   for (i = 0 ; i < rows ; i++) {
      temp1D = (double *) malloc(columns * sizeof(double)) ;
      DCrashOnNull(temp1D, "Could not make new matrix: create row") ;
      for (j = 0 ; j < columns ; j++ ) {
         temp1D[j] = 0.0 ;
      }
      M->matrix[i] = temp1D ;
   }

   return M ;
}



/* Read in data from a file */
dm_ptr LoadDense(char *filename) {
   FILE *ifile ;
   dm_ptr M ;
   int rows, columns, r, i, j ;
   double value ;

   ifile = fopen(filename, "r") ;
   if (ifile == NULL) return NULL ;

   r = fscanf(ifile, "%d %d", &rows, &columns) ;
   if (r < 2) {
      fprintf (stderr, "Bad file format\n") ;
      exit(1) ;
   }

   M = NewDenseMatrix(rows, columns) ;
   while(1) {
      r = fscanf(ifile, "%d %d %lf", &i, &j, &value) ;
      if (r < 3) break ;
      M->matrix[i][j] = value ;
   }
   fclose(ifile) ;

   return M ;
}


/* Store data in a file */
int StoreDense(char *filename, dm_ptr M) {
   int i, j ;
   double value ;
   FILE *ofile ;

   ofile = fopen(filename, "w") ;
   if (ofile == NULL) return 0 ;

   fprintf(ofile, "%d %d\n", M->rows, M->columns) ;

   for (i = 0 ; i < M->rows ; i++) {
      for (j = 0 ; j < M->columns ; j++) {
         value = M->matrix[i][j]  ;
         if (value != 0.0) {
            fprintf(ofile, "%d %d %lf\n", i, j, value) ;
         }
      }
   }

   fclose(ofile) ;
   return 1 ;
}



/* Duplicate entire data structure */
dm_ptr CopyDense (dm_ptr M) {
   int i, j ;
   dm_ptr New ;

   New = NewDenseMatrix(M->rows,M->columns) ;

   for (i = 0 ; i < M->rows ; i++) {
      for (j = 0 ; j < M->columns ; j++) {
         New->matrix[i][j] = M->matrix[i][j] ;
      }
   }
   
   return New ;
}



/* Free entire data structure */
void FreeDense (dm_ptr M) {
   int i ;

   assert(M != NULL) ;

   for (i = 0 ; i < M->rows ; i++) {
      free( M->matrix[i] ) ;
   }
   
   free(M->matrix) ;
   free(M) ;
}



/* Do the addition */
dm_ptr AddDense (dm_ptr A, dm_ptr B) {
   int i, j, rows, columns ;
   dm_ptr New ;

   if (A->rows != B->rows || A->columns != B->columns)
      return NULL ;

   rows = A->rows ;
   columns = A->columns ;

   New = NewDenseMatrix(rows, columns) ;

   for (i = 0 ; i < rows ; i++) {
      for (j = 0 ; j < columns ; j++) {
          New->matrix[i][j] = A->matrix[i][j] + B->matrix[i][j] ;
      }
   }

   return New ;
}


/* Do the subtraction */
dm_ptr SubtractDense (dm_ptr A, dm_ptr B) {
   int i, j, rows, columns ;
   dm_ptr New ;

   if (A->rows != B->rows || A->columns != B->columns)
      return NULL ;

   rows = A->rows ;
   columns = A->columns ;

   New = NewDenseMatrix(rows, columns) ;

   for (i = 0 ; i < rows ; i++) {
      for (j = 0 ; j < columns ; j++) {
         New->matrix[i][j] = A->matrix[i][j] - B->matrix[i][j] ;
      }
   }

   return New ;
}


/* Compare absolute value differences */
double CompareDense (dm_ptr A, dm_ptr B) {
   int i, j, rows, columns ;
   double diff, sum ;

   if (A->rows != B->rows || A->columns != B->columns)
      return NULL ;

   rows = A->rows ;
   columns = A->columns ;
   sum = 0.0 ;

   for (i = 0 ; i < rows ; i++) {
      for (j = 0 ; j < columns ; j++) {
         diff = A->matrix[i][j] - B->matrix[i][j] ;
         if (diff > 0) {
            sum += diff ;
         } else {
            sum -= diff ;
         }
      }
   }

   return sum ;
}


/* Do the multiplication */
dm_ptr MultiplyDense(dm_ptr A, dm_ptr B) {
   dm_ptr New ;
   int rows, columns, Acolumns ;
   int i, j, k ;
   double sum ;

   if (A->columns != B->rows) return NULL ;

   rows = A->rows ;
   Acolumns = A->columns ;
   columns = B->columns ;

   New = NewDenseMatrix(rows, columns) ;

   for (i = 0 ; i < rows ; i++) {

      for (j = 0 ; j < columns ; j++) {

         sum = 0.0 ;
         for (k = 0 ; k < Acolumns ; k++) {
            sum += A->matrix[i][k] * B->matrix[k][j] ;
         }
         New->matrix[i][j] = sum ;
      }
   }

   return New ;
}


/* Stupid iterative implementation */
/* Your project should do this recursively */

dm_ptr PowerDense(dm_ptr M, int n) {
   dm_ptr Product, temp ;
   int i ;

   if (n <= 0) return NULL ;

   Product = CopyDense(M) ;
   
   for (i = 2 ; i <= n ; i++){
      temp = MultiplyDense(Product, M) ;
      FreeDense(Product) ;
      Product = temp ;
   }

   return Product ;
}