/***************************************************************************
CANNON'S ALGORITHM FOR PARALLEL MATRIX MULTIPLICATION
 computes square matrix multiplication of size upto 64x64.
 This memory efficient algorithm uses checkerboard partitioning such that
 every process stores only one block of each operand at anytime.
 Using dynamic allocation will surely exploit the memory efficiency of this
 algorithm. However, this program uses static arrays to make it easier to
 implement the algorithm (note: some "strange" problems with "drifting"
 pointers forces the static array implementation)
 The program also measures the computation time and the total running time
 by using ftime routine of libbsd.a library (uses -lbsd option).

COMPILATION:
 mc -sync mpc -lbsd -o executable_name this_file_name

PARAMETERS:
 -procs n => uses n processors to run the program.
          This parameter is added automatically by the compiler.
          The default value is set in MP_PROCS environment variable.
 -input inputfile => reading the data of the matrix from text file inputfile
                  The first integer in the textfile is the number of rows
                  in both operand matrices.
 -output outputfile => send the output to outputfile.
                    If not given, output will be sent to output.dat
 -block b => uses block of size b in computation.
          The block size is limited to 16.
 -pro p   => to tell the program to use p processors .
          This parameter is necessary because the compiler-added -procs isn't 
          accessible to the program. 
          If not given, MP_PROCS will be used.
          p should be the same number as the number given to -procs parameter.
****************************************************************************/
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<sys/types.h>
#include<sys/timeb.h>
#include"maisie.h"

#define MAXSIZE 64                       /* maximum size of matrix is 64x64 */
#define MAXBLOCKSIZE 16                      /* maximum block size is 16x16 */

/****************************************************************************
entity block_processor
 represents one process of the parallel matrix multiplication.
INPUT:
 - myA and myB: the matrices
 - starting_row, starting_col:index of the the block in the big result matrix
 - num_of_blocks: total number of blocks that the matrix is divided into.
 - blocksize: the size of the block.
 - print_processor: printing entity.
OUTPUT:
 a block of the result of the matrix multiplication, sent to print_processor
****************************************************************************/
entity block_processor{myA, myB, starting_row, starting_col, num_of_blocks,
                       blocksize, print_processor}
int myA[MAXBLOCKSIZE][MAXBLOCKSIZE];
int myB[MAXBLOCKSIZE][MAXBLOCKSIZE];
int starting_row;
int starting_col;
int num_of_blocks;
int blocksize;
ename print_processor;
{
int i, j, k, shifting,
    myC[MAXBLOCKSIZE][MAXBLOCKSIZE];
ename left_processor;
ename up_processor;
message row_shift {int new_A[MAXBLOCKSIZE][MAXBLOCKSIZE];} rs;
message col_shift {int new_B[MAXBLOCKSIZE][MAXBLOCKSIZE];} cs;
message your_neighbor {ename left_processor; ename up_processor;} initmsg;


/* WAITING FOR THE NEIGHBORS' ADDRESS */
wait until mtype (your_neighbor)
 {
  left_processor = msg.your_neighbor.left_processor;
  up_processor   = msg.your_neighbor.up_processor;
 }

/* COMPUTATION OF INITIAL DATA */
for(i=0; i<blocksize; i++)
 for(j=0; j<blocksize; j++)
  {
   myC[i][j] = 0;
   for(k=0; k<blocksize; k++)
     myC[i][j] += myA[i][k] * myB[k][j];
  }

/* SHIFTING STATES */
for(shifting=0; shifting < num_of_blocks-1; shifting++)
 {
  invoke left_processor with row_shift{myA};
  invoke up_processor   with col_shift{myB};
  wait until mtype(row_shift) 
    for(i=0; i<blocksize; i++)
     for(j=0; j<blocksize; j++)
       myA[i][j] = msg.row_shift.new_A[i][j];
  wait until mtype(col_shift)
    for(i=0; i<blocksize; i++)
     for(j=0; j<blocksize; j++)
       myB[i][j] = msg.col_shift.new_B[i][j];
  for(i=0; i<blocksize; i++)
   for(j=0; j<blocksize; j++)
    for(k=0; k<blocksize; k++)
      myC[i][j] += myA[i][k] * myB[k][j];
 }

/* SEND THE RESULT */
invoke print_processor with result{starting_row, starting_col, myC};
return 1;

}/* end of entity block_processor */


entity driver{argc, argv}
int argc;
char **argv;
{

FILE *input, *outputfile;
int i, j, k, l, block, 
    starting_row, starting_col, row, col, 
    num_of_blocks, blocksize, num_of_procs, 
    most_digit, digit, output,
    A[MAXSIZE][MAXSIZE], B[MAXSIZE][MAXSIZE], 
    blockA[MAXBLOCKSIZE][MAXBLOCKSIZE], blockB[MAXBLOCKSIZE][MAXBLOCKSIZE];
char format[8];
char *tzstr = "TZ=PST8PDT"; /* pacific standard & daylight savings */
struct timeb Tp;
double start_time, finish_time, prog_start_time;
ename process[MAXSIZE][MAXSIZE];
message result {int row_num; int col_num;
                int data[MAXBLOCKSIZE][MAXBLOCKSIZE];} res;

/* GETTING PROGRAM STARTING TIME */
putenv(tzstr);
tzset();
ftime(&Tp);
prog_start_time = (double) Tp.time + (double) Tp.millitm / 1000;

output=0; 
num_of_procs = atoi(getenv("MP_PROCS"));    /* default number of processors */
blocksize = 1;                                         /* default blocksize */
/* GETTING THE PARAMETERS */
for(i=0; i<argc; i++)
{
 if (!strcmp(argv[i], "-pro"))          /* -pro number_of_processors_to_use */
   num_of_procs=atoi(argv[++i]);
 else if (!strcmp(argv[i], "-input"))               /* -input inputfilename */
        input=fopen(argv[++i], "rt");
      else if (!strcmp(argv[i], "-output"))       /* -output outputfilename */
             output= i+1;
           else if (!strcmp(argv[i], "-block"))         /* -block blocksize */
                  blocksize = atoi(argv[++i]);
}

/* READING DATA FROM THE INPUT FILE */
print("READING THE DATA");
fscanf(input, "%i", &row);
if (row > MAXSIZE)
 {
  print("ERROR: the size of the array is too large!");
  terminate();
 }
if (row % blocksize)
 {
  print("ERROR: array size must be divisible by blocksize!");
  terminate();
 }

print("%ix%i matrix multiplication with %ix%i blocks on %i processors",
       row, row, blocksize, blocksize, num_of_procs);
col = row;
num_of_blocks = row / blocksize;
for(i=0; i<row; i++)
 for(j=0; j<col; j++)
   fscanf(input, "%i", &A[i][j]);
for(i=0; i<row; i++)
 for(j=0; j<col; j++)
   fscanf(input, "%i", &B[i][j]);
fclose(input);

/* GET THE STARTING CALCULATION TIME */
ftime(&Tp);
start_time = (double) Tp.time + (double) Tp.millitm / 1000;
print("COMPUTING");

/* INITIAL DATA ARRANGEMENT */
for(i=0,starting_row=0,block=0; i<num_of_blocks; i++,starting_row+=blocksize)
 for(j=0,starting_col=0;j<num_of_blocks; j++,block++,starting_col+=blocksize)
  {
   for(k=0; k<blocksize; k++)             /* k and l is index of the block */
    for(l=0; l<blocksize; l++)
     {
      blockA[k][l] = A[starting_row+k][((starting_col+starting_row)%col) + l];
      blockB[k][l] = B[((starting_row+starting_col)%row) + k][starting_col+l];
     } 
   process[i][j] =new block_processor{blockA, blockB,
                                      starting_row, starting_col,
                                      num_of_blocks, blocksize, self}
                                      at (block % num_of_procs);
  }

/* TELLING THE PROCESSES THEIR LEFT & UP NEIGHBORS */
for(i=0; i<num_of_blocks; i++)
 for(j=0; j<num_of_blocks; j++)
   invoke process[i][j] with 
          your_neighbor{process[i][(j+num_of_blocks-1)%num_of_blocks], 
                        process[(i+num_of_blocks-1)%num_of_blocks][j]};

/* GETTING THE RESULTS FROM THE PROCESSES */
most_digit = 0;
for(i=0; i<num_of_blocks; i++)
 for(j=0; j<num_of_blocks; j++)
   wait until mtype(result)
     for(k=0; k<blocksize; k++)
      for(l=0; l<blocksize; l++)
       {
        A[msg.result.row_num+k][msg.result.col_num+l] = msg.result.data[k][l];
        digit = (int) ceil(log10(abs(msg.result.data[k][l])));
        if (digit>most_digit)
          most_digit = digit;
       }

/* RESULTS IS IN => SAVE THE END OF COMPUTATION TIME */
ftime(&Tp);
finish_time = (double) Tp.time + (double) Tp.millitm / 1000;

/* PRINTING THE RESULT TO A TEXT FILE */
print("PRINTING");
most_digit = most_digit + 2;
sprintf(format, "%%%ii", most_digit);
if (output)
  outputfile=fopen(argv[output], "wt");
else
  outputfile=fopen("output.dat", "wt");  /* default output file: output.dat */
for(i=0; i<row; i++)
 {
  for(j=0; j<col; j++)
    fprintf(outputfile, format, A[i][j]);
  fprintf(outputfile, "\n");
 }
fclose(outputfile);

/* TIMING RESULT */
ftime(&Tp);
printf("Computation time = %6.4lf\nTotal running time = %6.4lf", 
      finish_time - start_time,
      (double) Tp.time + (double) Tp.millitm / 1000 - prog_start_time);

terminate();
}/* end of entity driver */