class Data{
  //the slice of the domain we are iterating over
  const D : domain(2);
  // original D expanded by 1 to hold neighbouring rows / halo regions
  const bigD : domain(2) = D.expand(1, 0);
  var temps: [bigD] real;
  var cond: [D] real;
  var k : int;
}

proc do_compute() {
  /* your main function */

  var r          : results;
  var t          : Timer;
  var b       = new Barrier(numLocales, reusable=true);
  var diff       : [tempsD] real;
  var indexes = calculate_indexes();
  var niter      : int;
  var localeData : [LocaleSpace] Data;


  //setup data across locales
  for l in LocaleSpace{
    on Locales[l]{
      const D : domain(2) = {indexes[l].. (indexes[l+1])-1, 1..M};
      localeData[l] = new Data(D);
      const bigD = localeData[l].bigD;
      if(l == 0){
        localeData[l].temps[{indexes[l]..indexes[l+1], 1..M}] = tinit[{indexes[l]..indexes[l+1], 1..M}];
        localeData[l].temps[0,..]  = tinit[1,..];
      }
      else if (l == (numLocales-1)){
        var x : int = bigD.last(1);
        localeData[l].temps[{indexes[l]-1 ..(indexes[l+1])-1, 1..M}] = tinit[{indexes[l]-1.. (indexes[l+1])-1, 1..M}];
        localeData[l].temps[x,..]            = tinit[x-1,..];
      }
      else{
        localeData[l].temps[bigD] = tinit[bigD];
      }
      localeData[l].cond = tcond[D];
    }
  }
  t.start();
  coforall l in LocaleSpace{
    on Locales[l]{
      const D : domain(2)    = localeData[l].D;
      const bigD : domain(2) = localeData[l].bigD;
      const firstRow         = {(indexes[l+1]-1)..(indexes[l+1]-1), 1..M};
      const lastRow          = {(indexes[l])..(indexes[l]) , 1..M};
      const cond             = localeData[l].cond;
      var temps_new          : [D] real;
      var k                  : int = 0;
      var maxdiff            : real;
      do{
        const temps = localeData[l].temps;
        forall (i,j) in D{
          //offset for accessing neighbour in cylindrical surface
          const pos : int = if (j == M) then -(M-1) else 1;
          const neg : int = if (j == 1) then  (M-1) else -1;
          temps_new[i,j] = (temps[i,j] * cond[i,j])
            //direct neighbours
            + (((1 - cond[i,j]) * direct) *
               ((temps[i+1,j] +
                 temps[i-1,j] +
                 temps[i,j+pos] +
                 temps[i,j+neg]) / 4) )
            //diagonal neighbours
            + (((1 - cond[i,j]) * diag) *
               ((temps[i+1, j+pos] +
                 temps[i+1, j+neg] +
                 temps[i-1, j+pos] +
                 temps[i-1, j+neg] ) / 4) );

          // update current max/min values
          diff[i,j] =  abs(temps_new[i,j] - temps[i,j]);
        }
        localeData[l].temps[D] = temps_new;
        // exchange neighbour rows
        if (l != 0) then
          localeData[l-1].temps[lastRow] = temps_new[lastRow];
        if(l != numLocales-1) then
          localeData[l+1].temps[firstRow] = temps_new[firstRow];

        on diff do maxdiff = max reduce diff;

        k += 1;
        b.barrier(); //wait for all locales to exchange rows
      }while(k<I && maxdiff>E);
      localeData[l].k = k;
    }
  }
  t.stop();
}