source: src/smoother/gsrb.cpp@ ab63b6

/**
 * @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/timer.hpp"
#include "base/stencil.hpp"
#include "comm/comm.hpp"
#include "grid/grid.hpp"
#include "smoother/gsrb.hpp"

using namespace VMG;

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

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

  Timer::Start("SmootherWithoutCommunication");

  Index i;
  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();
  const int off = rhs.Global().BeginLocal().Sum() - rhs.Local().HasHalo1().Sum();

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

        temp = prefactor_inv * rhs.GetCorrectedVal(i);

        for (Stencil::iterator iter=A.begin(); iter!=A.end(); ++iter)
          temp -= iter->Val() * sol.GetVal(i + iter->Disp());

        sol(i) = temp * diag_inv;

      }

  Timer::Stop("SmootherWithoutCommunication");

  MG::GetComm()->CommToGhosts(sol);

  Timer::Start("SmootherWithoutCommunication");

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

        temp = prefactor_inv * rhs.GetCorrectedVal(i);

        for (Stencil::iterator iter=A.begin(); iter!=A.end(); ++iter)
          temp -= iter->Val() * sol.GetVal(i + iter->Disp());

        sol(i) = temp * diag_inv;

      }

  Timer::Stop("SmootherWithoutCommunication");
}