/*************
**
**  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)
**
**  Note: compile with -lm
**  Adapted for Maisie v. 3.0 on 2-25-97
**
**************/
#include <stdio.h>
#include <math.h>
#include "maisie.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) */

message answer { int k; float row[N];};

entity driver {int argc, char **argv}
{
  FILE *fp[2];
  int i,j;
  float temp, A[N][N];
  ename p_id[N];
  message answer answer;

  if ((fp[0] = fopen (INPUT_FILE, "r")) == NULL)
        {
    fputs("error opening the input file\n",stderr);
                terminate();
        }

  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 init_row {i*NR+j, A[i*NR+j]};

/** 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];
};

message link { ename prev; ename next; ename main;};
message start { };
message end   { };
message eliminator { int k; float row[N];};
message init_row   { int k; float row[N];};

entity processor {int myid}
{
  int i,j,k,count=0;
  float tf;
  ename previd, nextid;
  struct indx_row myrow[NR], temp;
  message link link;
  message eliminator eli;
  message init_row init_row;

  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 init_row=mtype(init_row);
    for (j=0; j<N; j++)
      myrow[i].row[j] = init_row.row[j];
    myrow[i].k = init_row.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};

  /* 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]*eli.row[j];
          myrow[i].row[count] = 0;
        }
        count++;
      } or mtype (end)
        terminate();
    }
  } while (1);
}