source: CIVL/examples/messagePassing/diffusion_1d.cvl@ b2ddb2e

/* Computing the solver for heat equation in multithreads, then compare the 
 *  result with the solver got in a sequencial way.
 *
 *  Command line example:
 *  civl run -inpuNPROCS=3 diffusion_1d.cvl
 **/ 

#include<civlc.h>
#include "mp_root.cvh"
#define FROM_LEFT 1           /* mpi tag which means the source is send or receive
                                 from left chunk*/
#define FROM_RIGHT 2          /* mpi tag which means the source is send or receive
                                 from right chunk*/
#define Threshold 1     /* the threshold of difference */
#define k  0.3                /* k = alpha^2 * dt/(dx^2) */
#define ELENUM 8              /* the length of the array */
#define RELEASE 0
#define HOLD 1

double doubleAbs(double v){
  if(v < 0)
    return -v;
  else
    return v;
}

/* return boundary values */
double boundaryValues(int x){
  return 4.0;
}

/* report error */
void error(){
  return;
}

/* update temprature function 
   returns the max difference between the previous one and updated one*/
double update(double * u, double * u_new, int start, int end){
    int i;
    int chunk_length = end - start + 1;

    // chunk_length = (the length of u or u_new) - 2
    for(i=1; i<chunk_length+1; i++)
      u_new [i] = u[i] + k * (u[i-1] + u[i+1] - 2*u[i]);

    // compute the difference between u and u_new
    double max_diff = 0;
    double diff;

    for(i=1; i<chunk_length+1; i++){
      diff = doubleAbs(u_new[i] - u[i]);
      if(diff > max_diff)
        max_diff = diff;
      // update u
      u[i] = u_new[i];
    }
    return max_diff;
}

/* Communicate with left chunk and right chunk,
   send the first interior element to left chunk and receive for 
   the last interior element */
void communicate(double * u, int left, int right, int start, 
                 int end, int rank){
  #include "mp_proc_diffusion.cvh"

  int chunk_length = end - start + 1;
  // the most left or right chunks just need to exchange with one side 
    if(left == -1 && right == -1)
      return;
    else if(left == -1){ 
      send(&u[chunk_length], 1, right, FROM_LEFT, rank); 
      recv(&u[chunk_length + 1], 1, right, FROM_RIGHT, rank);
    }
    else if(right == -1){
      send(&u[1], 1, left, FROM_RIGHT, rank); 
      recv(&u[0], 1, left, FROM_LEFT, rank);
    }else{
      send(&u[chunk_length], 1, right, FROM_LEFT, rank);
      send(&u[1], 1, left, FROM_RIGHT, rank);
      recv(&u[chunk_length + 1], 1, right, FROM_RIGHT, rank);
      recv(&u[0], 1, left, FROM_LEFT, rank);
    }
}

/* compute the value of start, end, left chunk and right chunk
   return an array with those values in such order*/
void ChunkManager(int length, int nprocs, int rank, int* results){
  if(length < 3 || nprocs > (length - 2))
    error();
  
    int remainder = (length - 2) % nprocs;
    int chunk_length =(length - 2) / nprocs;

    int start = rank * chunk_length + 1; 

    // the last chunk takes the remainder
    if(rank == nprocs - 1){
      chunk_length = chunk_length + remainder;
    }

    int end = start + chunk_length - 1;

    int left_chunk = rank - 1;
    int right_chunk = rank + 1;
    if(right_chunk >= nprocs)
      right_chunk = -1;

    results[0] = start;
    results[1] = end;
    results[2] = left_chunk;
    results[3] = right_chunk;

  }

/* the numbers of elements in u and u_new are 2 more than chunk_length because of 
   ghost elements */
int initChunk(double * u, double * u_new, int start, int end,
         int left_chunk, int right_chunk){
    int i;
    int chunk_length = end - start + 1; // the number of interior elements

    /* initiate interior array */
    for(i=1; i< chunk_length+1; i++){
      u[i] = 0;
      u_new[i] = 0;
    }

  /* for the most left or right chunks, initiate the first ghost ele and 
     the last ghost ele */
    if(left_chunk == -1){
        u[0] = boundaryValues(0);
        u_new[0] = boundaryValues(0);
    }
    if(right_chunk == -1){
      i = chunk_length + 1;
      u[i] = boundaryValues(i);
      u_new[i] = boundaryValues(i);
    }

    return chunk_length;
}

double total_diff = 0.0;        /* shared by all processes */
int reduce_count = 0;           /* used for sync */
short __release = HOLD;         /* used for sync */
double globle_u[ELENUM];

/* combine all separate chucks in all processes  */
double * combineU(int start, int end, double * u){
  globle_u[0] = boundaryValues(0);
  globle_u[ELENUM - 1] = boundaryValues(ELENUM - 1);
  int i = start;
  int j = 1;
 
  for(; i < end+1; i++){
    globle_u[i] = u[j];
    j++;
  }
}

/* synchronization */
void procsHold(int nprocs){
  $atomic{
    if(reduce_count < nprocs -1){
      reduce_count++;
    }else{
      reduce_count++;
      __release = RELEASE;
    }
  }
  $when (__release == RELEASE);
}
/* synchronization */
void procsRelease(){
  $atomic{
    if(reduce_count > 1){
      reduce_count--;
    }else{
      reduce_count--;
      __release = HOLD;
    }
  }
 $when (__release == HOLD);
}

/* do the same thing in a sequential way*/
void seqDiffusion1d(double * seq_u){
  double diff,max_diff = 0.0;
  double u[ELENUM], u_new[ELENUM];
  int i;

  //Initiate                                                              
  u[0] = u_new[0] = 4;//boundaryValues(0);
  u[ELENUM-1] = u_new[ELENUM-1] = 4 ;//boundaryValues(ELENUM - 1);
  for(i=1; i<ELENUM-1; i++){
    u[i] = 0;
    u_new[i] = 0;
  }

  //Jacobi Iteration             
  while(1){
    //update                                         
    for(i=1; i<ELENUM-1; i++){
      u_new[i] = u[i] + k * (u[i-1] + u[i+1] - 2*u[i]);
    }

    for(i=1; i<ELENUM-1; i++){
      diff = doubleAbs(u_new[i] - u[i]);
      if(diff > max_diff)
        max_diff = diff;
      // update u
      u[i] = u_new[i];
    }

    //termination                                    
    if(max_diff*2 <= Threshold)
      break;
    else
      max_diff = 0;
  }
  for(i =0 ;i<ELENUM; i++){
    seq_u[i] = u[i];
  }
}
double seq_u[ELENUM];
void MPI_Process (int rank){
  
  $when (__start);
  double diff, myTotalDiff;;  /* the max difference between previous 
                                      function and updated function*/
  int nprocs;          /* number of processes*/
  int start, end;      /* the index of the chunk*/
  int left_chunk, right_chunk; /* the index of the left chunk 
                                  and right chunk*/
  int chunk_length;   /* number of elements in this chunk*/
  int temp[4];        /*temp buffer for start, end, top_chunk and bottom_chunk**/


  nprocs = NPROCS;
  $atomic{
    ChunkManager(ELENUM, nprocs, rank, temp);
  }
  start = temp[0];
  end = temp[1];
  left_chunk = temp[2];
  right_chunk = temp[3];
  
  chunk_length = end - start + 1;
  
  double u[chunk_length + 2];
  double u_new[chunk_length + 2];
  $atomic{
    initChunk(u, u_new, start, end, left_chunk, right_chunk);
  }

  /* Jacobi Iterations*/
  while(1){
    communicate(u, left_chunk, right_chunk, start, end, rank);
    diff = update(u, u_new, start, end);

    //accumulate diff
    procsHold(nprocs);
    total_diff += diff;
    procsRelease();
    procsHold(nprocs);
    myTotalDiff = total_diff;
    procsRelease();

    if(myTotalDiff <= Threshold){
      break;
    }
    else{
      myTotalDiff = 0;
      total_diff = 0;
    }
  }
  procsHold(nprocs);
  $atomic{
    combineU(start,end,u);
  }
  procsRelease();
  //compare with seq results
  if(rank == 0){
    $atomic{ 
    int i;
    seqDiffusion1d(seq_u);
    $assert (seq_u[0] > 0);
    for(i=0; i<ELENUM; i++){
      double test_dif = doubleAbs((seq_u[i] - globle_u[i]));
      $assert (test_dif < 0.1);
      //      $assert (doubleAbs(seq_u[i] - globle_u[i]) < 1 );
      /* printf("difference between seq and parallel is %f\n",
         (seq_u[i] - globle_u[i])); */
    }
    }
  }
}