#ifndef JACOBI_THREAD_HPP
#define JACOBI_THREAD_HPP

#include <cassert>
#include <cmath>
#include "barrier.hpp"

template <typename Real, typename Matrix>
class JacobiThread {
   public:
      JacobiThread(Matrix a0, Matrix a1, Real eps,
	       AggregatingBarrier<Real>& barrier) :
	    a({a0, a1}), eps(eps), barrier(barrier) {
	 assert(a0.get_rows() == a1.get_rows() &&
	    a0.get_columns() == a1.get_columns());
      }
      void operator()() {
	 unsigned int rows = a[0].get_rows();
	 unsigned int cols = a[0].get_columns();
	 unsigned int index = 0;
	 Real localmax;
	 do {
	    localmax = 0;
	    Matrix& a0(a[index]);
	    Matrix& a1(a[1-index]);
	    for (unsigned int i = 1; i < rows-1; ++i) {
	       for (unsigned int j = 1; j < cols-1; ++j) {
		  a1(i, j) = (a0(i-1,j) + a0(i,j-1) +
		     a0(i,j+1) + a0(i+1,j)) / 4;
		  Real diff = fabs(a1(i, j) - a0(i, j));
		  if (diff > localmax) {
		     localmax = diff;
		  }
	       }
	    }
	    index = 1 - index;
	 } while (barrier.wait(localmax) > eps);
      }
   private:
      Matrix a[2];
      Real eps;
      AggregatingBarrier<Real>& barrier;
};

#endif