/* tropical matrix multiplication */

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

void random_ternary_matrix ( int nbrows, int nbcols, double **matrix );
/*
 * DESCRIPTION :
 *   Returns a random tropical matrix with three types of elements
 *   -inf, 0, or 1.  The matrix is stored row wise.
 *
 * ON ENTRY :
 *   nbrows   the number of rows of the matrix;
 *   nbcols   the number of columns of the matrix;
 *   matrix   allocated space for nbrows of pointers to rows of doubles,
 *            with every row allocated for nbcols doubles.
 *
 * ON RETURN :
 *   matrix   matrix with as many rows as the value of nbrows,
 *            and as many columns as the value of nbcols,
 *            with elements -inf, 0.0, or 1.0. */

void write_matrix ( int nbrows, int nbcols, double **matrix );
/*
 * DESCRIPTION :
 *   Writes the matrix with number of rows in nbrows
 *   and number of columns in nbcols. */

double max ( double x, double y);
/*
 * DESCRIPTION :
 *   Returns the maximum of x and y. */

void multiply
 ( int nbArows, int nbAcols, int nbBcols,
   double **Amatrix, double **Bmatrix, double **Cmatrix );
/*
 * DESCRIPTION :
 *   Multiplies two tropical matrices, Amatrix and Bmatrix.
 *   The result is in the Cmatrix.
 *
 * ON ENTRY :
 *   nbArows  number of rows of the Amatrix;
 *   nbAcols  number of columns of the Amatrix,
 *            number of rows of the Bmatrix;
 *   nbBcols  number of columsn of the Bmatrix.
 *   Amatrix  matrix with as many rows as nbArows
 *            and as many columns as nbAcols;
 *   Bmatrix  matrix with as many rows as nbAcols
 *            and as many columns as nbBcols.
 *
 * ON RETURN :
 *   Cmatrix  result of the multiplication of Amatrix with Bmatrix,
 *            with as many rows as nbArows and as many columns as nbBcols. */

int main ( int argc, char *argv[] )
{
   const int nbArows = 3;
   const int nbAcols = 5;
   const int nbBrows = nbAcols;
   const int nbBcols = 4;
   double *Amatrix[nbArows];
   double *Bmatrix[nbBrows];
   double *Cmatrix[nbArows];
   int i;

   for(i=0; i<nbArows; i++)
      Amatrix[i] = (double*) calloc(nbAcols, sizeof(double));
   for(i=0; i<nbBrows; i++)
      Bmatrix[i] = (double*) calloc(nbBcols, sizeof(double));
   for(i=0; i<nbArows; i++)
      Cmatrix[i] = (double*) calloc(nbBcols, sizeof(double));

   srand(time(NULL)); /* initialize seed with the time */
   random_ternary_matrix(nbArows, nbAcols, Amatrix);
   printf("The matrix A :\n");
   write_matrix(nbArows, nbAcols, Amatrix);

   random_ternary_matrix(nbBrows, nbBcols, Bmatrix);
   printf("The matrix B :\n");
   write_matrix(nbBrows, nbBcols, Bmatrix);

   multiply(nbArows, nbAcols, nbBcols, Amatrix, Bmatrix, Cmatrix);
   printf("The matrix A*B :\n");
   write_matrix(nbArows, nbBcols, Cmatrix);
   
   return 0;
}

void random_ternary_matrix ( int nbrows, int nbcols, double **matrix )
{
   int i,j;

   for(i=0; i<nbrows; i++)
      for(j=0; j<nbcols; j++)
         matrix[i][j] = log2(rand() % 3);
}

void write_matrix ( int nbrows, int nbcols, double **matrix )
{
   int i,j;

   for(i=0; i<nbrows; i++, printf("\n"))
      for(j=0; j<nbcols; j++)
         printf(" %.2f", matrix[i][j]);
}

double max ( double x, double y )
{
   if(x > y)
      return x;
   else
      return y;
}

void multiply
 ( int nbArows, int nbAcols, int nbBcols,
   double **Amatrix, double **Bmatrix, double **Cmatrix )
{
   int i,j,k;
   double sab;

   for(i=0; i<nbArows; i++)
      for(j=0; j<nbBcols; j++)
      { 
         Cmatrix[i][j] = -INFINITY;
         for(k=0; k<nbAcols; k++)
         {
            sab = Amatrix[i][k] + Bmatrix[k][j];
            Cmatrix[i][j] = max(Cmatrix[i][j], sab);
         }
      }
}