/**
 * @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 <sstream>

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

using namespace VMG;

static inline void ComputePartial(Grid& sol, Grid& rhs, const Stencil& mat,
				  const Index& begin, const Index& end,
				  const vmg_float& prefactor, const vmg_float& diag_inv,
				  const int& off)
{
  int i,j,k;
  vmg_float temp;
  Stencil::iterator iter;

  for (i=begin.X(); i<end.X(); ++i)
    for (j=begin.Y(); j<end.Y(); ++j)
      for (k=begin.Z()+((i+j+off)%2); k<end.Z(); k+=2) {

	temp = prefactor * rhs.GetVal(i,j,k);

	for (iter=mat.begin(); iter!=mat.end(); ++iter)
	  temp -= iter->Val() * sol.GetVal(i+iter->Disp().X(),
					   j+iter->Disp().Y(),
					   k+iter->Disp().Z());

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

      }
}

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

#ifdef DEBUG
  sol.IsCompatible(rhs);
#endif

  const Stencil& mat = MG::GetDiscretization()->GetStencil();
  const vmg_float prefactor_inv = 1.0 / MG::GetDiscretization()->OperatorPrefactor(sol);
  const vmg_float diag_inv = 1.0 / mat.GetDiag();
  const int off = rhs.Global().BeginLocal().Sum() - rhs.Local().HasHalo1().Sum();
  Comm& comm = *MG::GetComm();

  /*
   * Compute first halfstep
   */

  // Start asynchronous communication
  comm.CommToGhostsAsyncStart(sol);

  // Smooth part not depending on ghost cells
  ComputePartial(sol, rhs, mat,
		 rhs.Local().Begin()+1, rhs.Local().End()-1,
		 prefactor_inv, diag_inv, off+1);

  // Finish asynchronous communication
  comm.CommToGhostsAsyncFinish();

  /*
   * Smooth near boundary cells
   */

  ComputePartial(sol, rhs, mat,
		 rhs.Local().Begin(),
		 Index(rhs.Local().Begin().X()+1, rhs.Local().End().Y(), rhs.Local().End().Z()),
		 prefactor_inv, diag_inv, off+1);

  ComputePartial(sol, rhs, mat,
		 Index(rhs.Local().End().X()-1, rhs.Local().Begin().Y(), rhs.Local().Begin().Z()),
		 rhs.Local().End(),
		 prefactor_inv, diag_inv, off+1);

  ComputePartial(sol, rhs, mat,
		 Index(rhs.Local().Begin().X()+1, rhs.Local().Begin().Y(), rhs.Local().Begin().Z()),
		 Index(rhs.Local().End().X()-1, rhs.Local().Begin().Y()+1, rhs.Local().End().Z()),
		 prefactor_inv, diag_inv, off+1);

  ComputePartial(sol, rhs, mat,
		 Index(rhs.Local().Begin().X()+1, rhs.Local().End().Y()-1, rhs.Local().Begin().Z()),
		 Index(rhs.Local().End().X()-1, rhs.Local().End().Y(), rhs.Local().End().Z()),
		 prefactor_inv, diag_inv, off+1);

  ComputePartial(sol, rhs, mat,
		 Index(rhs.Local().Begin().X()+1, rhs.Local().Begin().Y()+1, rhs.Local().Begin().Z()),
		 Index(rhs.Local().End().X()-1, rhs.Local().End().Y()-1, rhs.Local().Begin().Z()+1),
		 prefactor_inv, diag_inv, off+1);

  ComputePartial(sol, rhs, mat,
		 Index(rhs.Local().Begin().X()+1, rhs.Local().Begin().Y()+1, rhs.Local().End().Z()-1),
		 Index(rhs.Local().End().X()-1, rhs.Local().End().Y()-1, rhs.Local().End().Z()),
		 prefactor_inv, diag_inv, off+1);

  /*
   * Compute second halfstep
   */

  // Start asynchronous communication
  comm.CommToGhostsAsyncStart(sol);

  // Smooth part not depending on ghost cells
  ComputePartial(sol, rhs, mat,
		 rhs.Local().Begin()+1, rhs.Local().End()-1,
		 prefactor_inv, diag_inv, off);

  // Finish asynchronous communication
  comm.CommToGhostsAsyncFinish();

  /*
   * Smooth near boundary cells
   */

  ComputePartial(sol, rhs, mat,
		 rhs.Local().Begin(),
		 Index(rhs.Local().Begin().X()+1, rhs.Local().End().Y(), rhs.Local().End().Z()),
		 prefactor_inv, diag_inv, off);

  ComputePartial(sol, rhs, mat,
		 Index(rhs.Local().End().X()-1, rhs.Local().Begin().Y(), rhs.Local().Begin().Z()),
		 rhs.Local().End(),
		 prefactor_inv, diag_inv, off);

  ComputePartial(sol, rhs, mat,
		 Index(rhs.Local().Begin().X()+1, rhs.Local().Begin().Y(), rhs.Local().Begin().Z()),
		 Index(rhs.Local().End().X()-1, rhs.Local().Begin().Y()+1, rhs.Local().End().Z()),
		 prefactor_inv, diag_inv, off);

  ComputePartial(sol, rhs, mat,
		 Index(rhs.Local().Begin().X()+1, rhs.Local().End().Y()-1, rhs.Local().Begin().Z()),
		 Index(rhs.Local().End().X()-1, rhs.Local().End().Y(), rhs.Local().End().Z()),
		 prefactor_inv, diag_inv, off);

  ComputePartial(sol, rhs, mat,
		 Index(rhs.Local().Begin().X()+1, rhs.Local().Begin().Y()+1, rhs.Local().Begin().Z()),
		 Index(rhs.Local().End().X()-1, rhs.Local().End().Y()-1, rhs.Local().Begin().Z()+1),
		 prefactor_inv, diag_inv, off);

  ComputePartial(sol, rhs, mat,
		 Index(rhs.Local().Begin().X()+1, rhs.Local().Begin().Y()+1, rhs.Local().End().Z()-1),
		 Index(rhs.Local().End().X()-1, rhs.Local().End().Y()-1, rhs.Local().End().Z()),
		 prefactor_inv, diag_inv, off);
}