/*************
**
**  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 Parsec 1.0 on 9-5-97 by Monnica Terwilliger
**
**************/
#include <stdio.h>
#include <math.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];};
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);

struct indx_row
{
  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];

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

  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 **/
  send Link{ p_id[0], p_id[1], self} to p_id[0] ;
  for (i=1; i<NP-1; i++)
    send Link{ p_id[i-1], p_id[i+1], self} to p_id[i];
  send Link{ p_id[NP-2], p_id[NP-1], self} to p_id[NP-1];

/** send rows to the processors **/
  for (i=0; i<NP; i++)
    for (j=0; j<NR; j++)
      send Init_row {i*NR+j, A[i*NR+j]} to p_id[i];

/** start computing **/
  send Start {} to p_id[0];

/** wait for the answer **/
  for (i=0; i<N; i++)
  {
    receive (Answer 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++)
    send End { } to p_id[i] ;
}

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

  receive (Link tmp)
  {	link=tmp;
  	previd = link.prev;
  	nextid = link.next;
  }

/** receive the rows from the driver **/
  for (i=0; i<NR; i++)
  {
    receive (Init_row init_row)
    {	for (j=0; j<N; j++)
      		myrow[i].row[j] = init_row.row[j];
    	myrow[i].k = init_row.k;
    }
  }

  do
  {
    receive (Start start) when (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) send 
          	Eliminator{myrow[k].k, myrow[k].row} to nextid;

  /* 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) send Start { } to nextid;

        for (i=0; i<NR; i++)
          send Answer { myrow[i].k, myrow[i].row} to link.main;
      } or receive (Eliminator temp) {
        eli=temp;
        if (myid != NP-1) send eli to nextid;

      /* 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 receive (End signal)
        	pc_exit(1);
 
  } while (1);
}