/**
 * @file   gsrb.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Apr 18 13:08:20 2011
 *
 * @brief  Gauss-Seidel Red Black method
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "base/discretization.hpp"
#include "base/stencil.hpp"
#include "grid/grid.hpp"
#include "smoother/gsrb.hpp"

using namespace VMG;

void GaussSeidelRB::Compute(Grid& sol, Grid& rhs)
{
#ifdef DEBUG
  sol.IsConsistent();
  rhs.IsConsistent();
  sol.IsCompatible(rhs);
#endif

  int i,j,k;
  vmg_float temp;
  const Stencil& A = MG::GetDiscretization()->GetStencil();
  const vmg_float prefactor_inv = 1.0 / MG::GetDiscretization()->OperatorPrefactor(sol);
  const vmg_float diag_inv = 1.0 / A.GetDiag();

  for (i=rhs.Local().Begin().X(); i<rhs.Local().End().X(); i++)
    for (j=rhs.Local().Begin().Y(); j<rhs.Local().End().Y(); j++)
      for (k=rhs.Local().Begin().Z() + ((i+j+1) % 2); k<rhs.Local().End().Z(); k+=2) {

	temp = prefactor_inv * rhs.GetCorrectedVal(i,j,k);

	for (Stencil::iterator iter=A.begin(); iter!=A.end(); iter++)
	  temp -= iter->val * sol.GetVal(i+iter->m, j+iter->n, k+iter->o);

	sol(i, j, k) = temp * diag_inv;

      }

  for (i=rhs.Local().Begin().X(); i<rhs.Local().End().X(); i++)
    for (j=rhs.Local().Begin().Y(); j<rhs.Local().End().Y(); j++)
      for (k=rhs.Local().Begin().Z() + ((i+j) % 2); k<rhs.Local().End().Z(); k+=2) {

	temp = prefactor_inv * rhs.GetCorrectedVal(i,j,k);

	for (Stencil::iterator iter=A.begin(); iter!=A.end(); iter++)
	  temp -= iter->val * sol.GetVal(i+iter->m, j+iter->n, k+iter->o);

	sol(i, j, k) = temp * diag_inv;

      }
}