/* C++ standard headers */
#include <cassert>
#include <cmath>
#include <cstdlib>
/* HPC library headers */
#include <hpc/aux/hsvcolor.hpp>
#include <hpc/matvec/gematrix.hpp>
#include <hpc/matvec/matrix2pixbuf.hpp>

template<typename T>
const T PI = std::acos(T(-1.0));
template<typename T>
const T E = std::exp(T(1.0));
template<typename T>
const T E_POWER_MINUS_PI = std::pow(E<T>, -PI<T>);

using namespace hpc;

template<typename T, template<typename> typename Matrix,
   Require<Ge<Matrix<T>>> = true>
T jacobi_iteration(const Matrix<T>& A, Matrix<T>& B) {
   T maxdiff{};
   #pragma omp parallel for reduction(max:maxdiff)
   for (std::size_t i = 1; i < B.numRows() - 1; ++i) {
      for (std::size_t j = 1; j < B.numCols() - 1; ++j) {
	 B(i, j) = 0.25 *
	    (A(i - 1, j) + A(i + 1, j) + A(i, j - 1) + A(i, j + 1));
	 T diff = std::fabs(A(i, j) - B(i, j));
	 if (diff > maxdiff) maxdiff = diff;
      }
   }
   return maxdiff;
}

int main(int argc, char** argv) {
   using namespace hpc::aux;
   using namespace hpc::matvec;
   using T = double;
   using Matrix = GeMatrix<T>;

   Matrix A(100, 100, Order::RowMajor);
   Matrix B(A.numRows(), A.numCols(), Order::RowMajor);
   #pragma omp parallel for
   for (std::size_t i = 0; i < A.numRows(); ++i) {
      for (std::size_t j = 0; j < A.numCols(); ++j) {
	 if (j == 0) {
	    A(i, j) = B(i, j) = std::sin(PI<T> * (T(i)/(A.numRows()-1)));
	 } else if (j == A.numCols() - 1) {
	    A(i, j) = B(i, j) = std::sin(PI<T> * (T(i)/(A.numRows()-1))) *
	       E_POWER_MINUS_PI<T>;
	 } else {
	    A(i, j) = B(i, j) = 0;
	 }
      }
   }

   T eps = 1e-6;
   T maxdiff{};
   do {
      maxdiff = jacobi_iteration(A, B);
      maxdiff = jacobi_iteration(B, A);
   } while (maxdiff >= eps);

   auto pixbuf = create_pixbuf(A, [](T val) -> HSVColor<float> {
      return HSVColor<float>((1-val) * 240, 1, 1);
   }, 8);
   gdk_pixbuf_save(pixbuf, "jacobi.jpg", "jpeg", nullptr,
      "quality", "100", nullptr);
}