/***************************************************************/
/*  matmult - parallel matrix multiplication algorithm         */
/*            Written by Mike Patterson for CS133, Winter 96   */
/*            Adapted for Maisie v. 3.0 on 2-27-97             */
/*                                                             */
/* This program uses Cannon's algorithm to order the initial   */
/* A and B matrices so that parallel calculations can occur.   */
/* This algorithm first shifts the elements of row i in A, i   */
/* places to the left (MOD size of the row).  Column i in B    */
/* is also shifted i places up (MOD size of column). After the */
/* initial data shifting, A & B are checkerboard partitioned   */
/* into submatrices and parcelled out to entities running in   */
/* parallel.  After each iteration, data is shifted up one     */
/* row for the A submatrices and left one column for the B     */
/* submatrix.                                                  */
/*                                                             */
/* A and B matrices must be square for this program.  The      */
/* maximum dimension size allowed is 64 (or MAXDIM). For       */
/* parallelization, the initial matrices are divided into      */
/* sqr submatrices of dimension blocksize.  Maximum allowed    */
/* blocksize is 16 (defined as MAXSUBDIM).  Entities are       */
/* created to manage an A and a B submatrix, and return the    */
/* results to the driver routine at the end of their           */
/* computation. For this program, the matrices can only hold   */
/* integer data.                                               */
/*                                                             */
/* For each iteration:                                         */
/*                                                             */
/* Each mmult entity handles multiplication of the elements    */
/* in their submatrices and the message passing to their UP    */
/* and LEFT neighbors (following Cannon's algorithm).  They    */
/* also receive data from their DOWN and RIGHT neighbors.      */
/*                                                             */
/* Message traffic between the mmult entities consists of rows */
/* from the A submatrix and cols from the B submatrix, one per */
/* iteration, as the data is shifted for the calculations.     */
/*                                                             */
/* Command line args:                                          */
/*     -input filename  : the name of the input file           */
/*     -output filename : where to save resulting C matrix     */
/*     -block #         : dimension (blocksize) of submatrices */
/*                                                             */
/* The format expected in the input file is illustrated below: */
/* although comments have been added, do not use them in file. */
/*                                                             */
/* 3                  # dimension of A (and B) (3x3)           */
/* 1 2 3              # first row of matrix A                  */
/* 4 0 2                                                       */
/* 8 5 1                                                       */
/* 7 1 4              # first row of matrix B                  */
/* 1 1 8                                                       */
/* 0 3 3                                                       */
/*                                                             */
/* The resulting matric C is saved to the putput file, along   */
/* with dimensions of the initial matrices and the subsequent  */
/* submatrix dimension used.  Timing information is also       */
/* printed in the output, measured starting after the data has */
/* been read in by the driver, up until just before the driver */
/* stores the answer in the file.  Therefore all parallel      */
/* computation (including entity creation) is timed, as well   */
/* as the initial data shifting of the A and B matrices.  The  */
/* output file also contains the number of processes used.     */
/*                                                             */
/***************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "maisie.h"

#define UP 0
#define DOWN 1
#define LEFT 2
#define RIGHT 3

#define  SQR(x)  ((x)*(x))

#define MAXDIM 64
#define MAXSUBDIM 16

int Mod( int numerator, int denominator)
{
  int val = numerator % denominator;

  if (val < 0)
    val += denominator;
  
  return(val);
}
  
int FillSubMatrix(
	int matrix[MAXDIM][MAXDIM],
	int begin_i,
	int begin_j,
	int end_i, 
	int end_j,
	int submatrix[MAXSUBDIM][MAXSUBDIM])
{
  int i, j;

  for (i=begin_i;i<end_i;i++)
    for (j=begin_j;j<end_j;j++)
      submatrix[i-begin_i][j-begin_j] = matrix[i][j];
  
  return(1);
}

int ShiftMatrix(
	int matrix[MAXDIM][MAXDIM],
	int dim,
	int shift)
{
  int tmp[MAXDIM];
  int i, j;

  switch (shift) {
    case LEFT:
      for (i=1;i<dim;i++) {
        for (j=0;j<dim;j++)
          tmp[Mod(j-i,dim)] = matrix[i][j];
        for (j=0;j<dim;j++)
          matrix[i][j] = tmp[j];
      }
      break;
    case UP:
      for (j=1;j<dim;j++) {
        for (i=0;i<dim;i++)
          tmp[Mod(i-j,dim)] = matrix[i][j];
        for (i=0;i<dim;i++)
          matrix[i][j] = tmp[i];
      }
      break;
    default:  /* could have RIGHT and DOWN cases too */
      break;
  }
  return(1);
}

int maxIntSize(int val)
{
  int size = 1;     /* even '0' takes one space */

  if (val < 0) {
    val = abs(val);
    size++;         /* for the minus sign */
  }

  while (val > 9) {
    val = val / 10;
    size++;  
  }

  return(size);
}

int PrintMatrix(
	char *filename,
	int matrix[MAXDIM][MAXDIM],
	int dim,
	int blocksize,
	float running_time)
{
  FILE *fp;
  int i, j;
  int maxval;
  int padding;
  char formatstring[10];

  if ((fp = fopen(filename, "w")) == NULL) {
    printf("ERROR! Unable to open output file for writing!\n");
    return;
  }

  fprintf(fp, "Multiplication of %dx%d matrices, using ", dim, dim);
  fprintf(fp, "%dx%d submatrix size.\n", blocksize, blocksize);
  fprintf(fp, "Running time was %6.4f seconds ", running_time);
  fprintf(fp, "using %d processes\n\n", SQR(dim/blocksize));

  maxval = 0;
  for (i=0;i<dim;i++)
    for (j=0;j<dim;j++)
      if (maxval < abs(matrix[i][j])) 
        maxval = matrix[i][j];    /* not the abs value!! */
  
  padding = maxIntSize(maxval) + 1;

  /* This next line creates something like "%4d" for the print
     stmt below, replacing the 4 with the padding value */

  sprintf(formatstring, " %%%dd", padding);

  for (i=0;i<dim;i++) {
    for (j=0;j<dim;j++)
      fprintf(fp, formatstring, matrix[i][j]);
    fprintf(fp, "\n");
  }
  fprintf(fp, "\n");
  fclose(fp);
}

int InputData(
	char *filename,
	int *dim,
	int A[MAXDIM][MAXDIM],
	int B[MAXDIM][MAXDIM])
{
  FILE *fp;
  int i, j;

  if (strlen(filename) < 1) {
    printf("WARNING: No input filename given\n");
    return(0);
  }

  if ((fp = fopen(filename, "r")) == NULL) {
    printf("WARNING: InputData unable to open file\n");
    return(0);
  }

  fscanf(fp, "%d", dim);

  for (i=0;i<(*dim);i++)
    for (j=0;j<(*dim);j++)
      fscanf(fp, "%d", &(A[i][j]));

  for (i=0;i<(*dim);i++)
    for (j=0;j<(*dim);j++)
      fscanf(fp, "%d", &(B[i][j]));

  fclose(fp);
  return(1);
}

int ParseArgv(int argc, char **argv, char *inf, char *outf, int *blocksize)
{
  int i;

  if (argc < 3) {
    printf("USAGE: %s -input infile -output outfile ", argv[0]);
    printf("-block blocksize\n");
    return(0);
  }

  *blocksize = MAXSUBDIM;  /* initialize to safe values */
  inf[0] = '\0';
  outf[0] = '\0';

  for (i=0;i<argc;i++) {
    if (strcasecmp(argv[i], "-input")==0)
      strcpy(inf, argv[i+1]);
    if (strcasecmp(argv[i], "-output")==0)
      strcpy(outf, argv[i+1]);
    if (strcasecmp(argv[i], "-block")==0)
      *blocksize=atoi(argv[i+1]);
  }

  return(1);
}

int FillInMatrix(
	int blockrow,
	int blockcol,
	int matrix[MAXDIM][MAXDIM],
	int submatrix[MAXSUBDIM][MAXSUBDIM],
	int blocksize)
{
  int i, j;

  for (i=blockrow*blocksize; i<(blockrow+1)*blocksize; i++)
    for (j=blockcol*blocksize; j<(blockcol+1)*blocksize; j++)
      matrix[i][j] = submatrix[i-blockrow*blocksize][j-blockcol*blocksize];
}

message init { ename above; ename left; };
message result {int val[MAXSUBDIM][MAXSUBDIM]; int rno; int cno;};

entity driver{int argc, char **argv}
{
  int dim;                         /* dimension of matrices A, B & C */
  int i, j, k;                     /* generic counters and indices   */
  int A[MAXDIM][MAXDIM]; 
  int B[MAXDIM][MAXDIM]; 
  int C[MAXDIM][MAXDIM];
  int subA[MAXSUBDIM][MAXSUBDIM];  /* passable A submatrix */
  int subB[MAXSUBDIM][MAXSUBDIM];  /* passable B submatrix */
  int subnum;                      /* # of block rows (& cols)  */
  char infile[80], outfile[80];
  int blocksize;                   /* submatrix dimension  */
  int starti, endi, startj, endj;
  int subrow, subcol;

  clock_t st_time;                 /* initial clock time  */
  clock_t end_time;                /* end clock time    */
  float running_time;

  ename pid[MAXDIM][MAXDIM];

  message init Init;
  message result Result;
  
  if (!ParseArgv(argc, argv, infile, outfile, &blocksize))
    terminate();

  if (!InputData(infile, &dim, A, B)) {
    printf("program terminating\n");
    terminate();
  }

  if (blocksize > MAXSUBDIM) {
    printf("WARNING! blocksize [%d] is larger than ", blocksize);
    printf("maximum allowable [%d].\n", MAXSUBDIM);
    printf("Resetting blocksize to %d.\n", MAXSUBDIM);
    blocksize = MAXSUBDIM;
  }

  if (Mod(dim,blocksize) != 0) {
    printf("ERROR! blocksize [%d] is not an integer ", blocksize);
    printf("multiple of matrix dimension [%d]\n", dim);
    terminate();
  }

  subnum = dim/blocksize;        /* calc submatrix dimension */

  st_time = clock();             /* start timing now...      */

  ShiftMatrix(A, dim, LEFT);
  ShiftMatrix(B, dim, UP);

  for (i=0;i<subnum;i++) {
    starti = i*blocksize;
    endi = (i+1)*blocksize;
    for (j=0;j<subnum;j++) {
      startj = j*blocksize;
      endj = (j+1)*blocksize;
      FillSubMatrix(A, starti, startj, endi, endj, subA);
      FillSubMatrix(B, starti, startj, endi, endj, subB);

      pid[i][j] = new mmult{i,j,subA,subB,dim,blocksize,self} at (i*subnum+j);
    }
  }

  for (i=0;i<subnum;i++)
    for (j=0;j<subnum;j++)
      invoke pid[i][j] with
        init{pid[Mod(i-1,subnum)][j], pid[i][Mod(j-1,subnum)]};

  for (i=0; i<subnum; i++)
    for (j=0; j<subnum; j++)
      wait until
        mtype (result) st (msg.result.rno == i && msg.result.cno == j)
          FillInMatrix(i, j, C, msg.result.val, blocksize);
  
  end_time = clock();    /* stop timing */

  running_time = ((float)end_time - st_time)/CLOCKS_PER_SEC;
  PrintMatrix(outfile, C, dim, blocksize, running_time);
  terminate();
}

message newrow {int data[MAXDIM];};
message newcol {int data[MAXDIM];};

entity mmult{ 
	int myrow, int mycol,
	int myA[MAXSUBDIM][MAXSUBDIM],
	int myB[MAXSUBDIM][MAXSUBDIM],
	int dim,
	int blocksize,
	ename spokesperson}
{
  int myC[MAXSUBDIM][MAXSUBDIM];
  int tmpvect[MAXSUBDIM];        /* used to swap rows (cols) w/ neighbors */
  int row, col;
  ename toMyUp, toMyLeft;        /* neighbor id's above and left of this one */
  int iter;

  message init Init;
  message result Result;
  message newrow Newrow;
  message newcol Newcol;

  for (row=0;row<blocksize;row++)
    for (col=0;col<blocksize;col++)
      myC[row][col] = 0;

  wait until
    mtype (init) {
      toMyLeft = msg.init.left;
      toMyUp = msg.init.above;
    }

  for (iter=0;iter<dim;iter++) {

	/* First, perform matrix calculations */

    for (row=0;row<blocksize;row++)
      for (col=0;col<blocksize;col++)
        myC[row][col] += (myA[row][col]*myB[row][col]);

    /* Now shift cols of my A submatrix */

    for (col=0; col<blocksize; col++) {
      for (row=0; row<blocksize; row++)
        tmpvect[row] = myA[row][col];
  
      if (col==0)
        invoke toMyLeft with newcol{tmpvect};
      else
        for (row=0; row<blocksize; row++)
          myA[row][col-1] = tmpvect[row];
    }

    wait until
      mtype (newcol)
        for (row=0; row<blocksize; row++)
          myA[row][blocksize-1] = msg.newcol.data[row];
  
    /* and now shift rows for B submatrix */

    for (row=0; row<blocksize; row++) {
      for (col=0; col<blocksize; col++)
        tmpvect[col] = myB[row][col];
  
      if (row==0)
        invoke toMyUp with newrow{tmpvect};
      else
        for (col=0; col<blocksize; col++)
          myB[row-1][col] = tmpvect[col];
    }
  
    wait until
      mtype (newrow)
        for (col=0; col<blocksize; col++)
          myB[blocksize-1][col] = msg.newrow.data[col];
  }

  /* Done with matrix calculations, pass myC back to spokesperson */

  invoke spokesperson with result{myC, myrow, mycol};
}