/**
 * @file   gs.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Apr 18 13:07:44 2011
 *
 * @brief  Gauss-Seidel method
 *
 */

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

#include <cassert>

#include "base/discretization.hpp"
#include "base/stencil.hpp"
#include "grid/multigrid.hpp"
#include "smoother/gs.hpp"
#include "mg.hpp"

using namespace VMG;

void GaussSeidel::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(); k<rhs.Local().End().Z(); k++) {

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

      	for (Stencil::iterator iter=A.begin(); iter!=A.end(); iter++) {

	  assert(i+iter->m >= 0 && i+iter->m < sol.Local().SizeTotal().X());
	  assert(j+iter->n >= 0 && j+iter->n < sol.Local().SizeTotal().Y());
	  assert(k+iter->o >= 0 && k+iter->o < sol.Local().SizeTotal().Z());

	  temp -= iter->val * sol.GetVal(i+iter->m, j+iter->n, k+iter->o);

      	}

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

      }
}