Here is the code:

#!/usr/bin/env rdmd
import std.stdio;
import std.concurrency;
import core.thread;
import std.datetime;
import std.conv;
import core.sync.barrier;



immutable int gridSize = 256;
immutable int MAXSTEPS = 50000; /* Maximum number of iterations */ immutable double TOL_VAL =0.00001; /* Numerical Tolerance */
immutable double omega =  0.376;
immutable double one_minus_omega = 1.0 - 0.376;


immutable int numberOfThreads = 2;


double MAX_FUNC(double a, double b)
{
  return a> b? a: b;
}

double ABS_VAL(double a)
{
  return a> 0? a: -a;
}

__gshared Barrier myLocalBarrier = null;
shared double[gridSize+2][gridSize+2] gridInfo;
shared double maxError = 0.0;

void main(string args[])
{

  for(int i=0; i<gridSize+2; i++)
  {
    for(int j=0; j<gridSize+2; j++)
    {
      if(i==0)
        gridInfo[i][j] = 1.0;
      else
        gridInfo[i][j] = 0.0;
    }
  }

  bool shouldCheck = false;
  bool isConverged = false;
  for(int iter = 1; iter <= MAXSTEPS; iter++)
  {
    shouldCheck = false;
    if(iter % 400 ==0)
    {
      shouldCheck = true;
      maxError = 0.0;
    }

      //This is Phase 1
    {
      myLocalBarrier = new Barrier(numberOfThreads+1);
      for (int cc=0; cc<numberOfThreads; cc++)
      {
        spawn(&SolverSlave, thisTid,cc, 0 ,shouldCheck);
      }


      myLocalBarrier.wait();
    }

     //This is Phase 2
    {
      myLocalBarrier = new Barrier(numberOfThreads+1);
      for (int cc=0; cc<numberOfThreads; cc++)
      {
       spawn(&SolverSlave, thisTid,cc, 1 ,shouldCheck );
      }

      myLocalBarrier.wait();
    }

    if( maxError <  TOL_VAL)
      {
        isConverged = true;
        break;
      }

  }
  if(isConverged)
    writeln("It converged");
  else
    writeln("It did not converge");
}



void SolverSlave(Tid owner, int myNumber, int remainder, bool shouldCheckHere)
{

  double sum =0;

  //Divide task among threads
  int iStart = ((myNumber*gridSize)/numberOfThreads) + 1;
  int iEnd =  (((myNumber+1)*gridSize)/numberOfThreads) ;



  for(int i=iStart; i<= iEnd; i++)
  {
    for(int j=1; j< gridSize+1; j++)
    {
      if( ((i+j)%2 ==remainder)) //Phase 1 or 2
      {
        sum = ( gridInfo[i  ][j+1] + gridInfo[i+1][j  ] +
            gridInfo[i-1][j  ] + gridInfo[i  ][j-1] )*0.25;

//Should not check everytime to reduce synchronization overhead
        if(shouldCheckHere)
        {
maxError = MAX_FUNC(ABS_VAL(omega *(sum-gridInfo[i][j])), maxError);
        }
gridInfo[i][j] = one_minus_omega*gridInfo[i][j] + omega*sum;
      }

    }
  }

  myLocalBarrier.wait();
}



Reply via email to