/*************
**
**  Dense Cholesky Factorization using block striped mapping
**
**
**  N is the number of unknowns (ie. the input matrix should be N x N)
**  NR is the number of rows per processor (ie. N / S == # of processors)
**
**************/
#include <stdio.h>
#include <math.h>
#include "mayc.h"
 
#define INPUT_FILE "in.matrix"
#define OUTPUT_FILE "out.matrix"
#define N 6     /* # of unknown */
#define NR 1    /* # of rows per block */
#define NP N/NR   /* # of processors (blocks) */
 
entity driver {argc,argv}
int argc;
char **argv;
{
  FILE *fp[2];
  int i,j;
  float temp, A[N][N];
  ename p_id[N];
  message answer { int k; float row[N];} answer;
 
  if ((fp[0] = fopen (INPUT_FILE, "r")) == NULL)
    fputs("error opening the input file\n",stderr);
 
  for (i=0; i<N; i++)
    for (j=0; j<N; j++)
    {
      fscanf(fp[0],"%f",&temp);
      if (i>j) A[i][j] = 0.0;
      else A[i][j] = temp;
    }
  fclose (fp[0]);
 
/** create the processors **/
  for (i=0; i<NP; i++)
    p_id[i] = new processor { i} at i;
 
/** establish link between the processors **/
  invoke p_id[0] with link{ p_id[0], p_id[1], self};
  for (i=1; i<NP-1; i++)
    invoke p_id[i] with link{ p_id[i-1], p_id[i+1], self};
  invoke p_id[NP-1] with link{ p_id[NP-2], p_id[NP-1], self};
 
/** send rows to the processors **/
  for (i=0; i<NP; i++)
    for (j=0; j<NR; j++)
      invoke p_id[i] with eliminator {i*NR+j, A[i*NR+j], 0, self};
 
/** start computing **/
  invoke p_id[0] with start;
 
/** wait for the answer **/
  for (i=0; i<N; i++)
  {
    wait until answer = mtype (answer);
    for (j=0; j<N; j++)
      A[answer.k][j] = answer.row[j];
  }
 
  if ((fp[1] = fopen (OUTPUT_FILE, "w")) == NULL)
    fputs("error opening the output file\n",stderr);
  for (i=0; i<N; i++)
    for (j=0; j<N; j++)
    {
      fprintf(fp[1],"%0.5f",A[i][j]);
      if (j==N-1) fprintf(fp[1],"\n");
      else fprintf(fp[1]," ");
    }
  fclose(fp[1]);
 
  printf("answer is written to '%s'\n",OUTPUT_FILE);
 
  for (i=0; i<NP; i++)
    invoke p_id[i] with end;
}
 
struct indx_row
{
  int k;
  float row[N];
};
 
entity processor {myid}
int myid;
{
  int i,j,k,count=0;
  float tf;
  ename previd, nextid;
  struct indx_row myrow[NR], temp;
  message link { ename prev; ename next; ename main;} link;
  message start;
  message end;
  message eliminator { int k; float row[N];} eli;
 
  wait until link = mtype(link);
  previd = link.prev;
  nextid = link.next;
 
/** receive the rows from the driver **/
  for (i=0; i<NR; i++)
  {
    wait until eli=mtype(eliminator);
    for (j=0; j<N; j++)
      myrow[i].row[j] = eli.row[j];
    myrow[i].k = eli.k;
  }
 
  do
  {
    wait until
    {
      mtype (start) st (qempty(eliminator))
      {
        for (k=0; k<NR; k++)
        {
  /* division step */
          myrow[k].row[count] = sqrt(myrow[k].row[count]);
          for (j=count+1; j<N; j++)
          {
            myrow[k].row[j] =
              myrow[k].row[j] / myrow[k].row[count];
          }
 
  /* send eliminator to other processors */
          if (myid != NP-1) invoke nextid
            with eliminator{myrow[k].k,
                   &myrow[k].row[count]::(N-count)*sizeof(float)};
 
  /* do the eliminating step */
          for (i=k+1; i<NR; i++)
          {
            for (j=count; j<N; j++)
              myrow[i].row[j] = myrow[i].row[j] -
                myrow[k].row[i+myid*NR]*myrow[k].row[j];
            myrow[i].row[count] = 0;
          }
          count++;
        }
 
        if (myid != NP-1) invoke nextid with start;
 
        for (i=0; i<NR; i++)
          invoke link.main with answer { myrow[i].k, myrow[i].row};
      } or eli = mtype (eliminator) {
        if (myid != NP-1) invoke nextid
          with eliminator = eli;
 
      /* do the elimination step */
        for (i=0; i<NR; i++)
        {
          for (j=count; j<N; j++)
            myrow[i].row[j] = myrow[i].row[j] -
                              eli.row[i+NR*myid-count]*eli.row[j-count];
          myrow[i].row[count] = 0;
        }
        count++;
      } or mtype (end)
        exit(0);
    }
  } while (1);
}