---

GenBiCGStabSmootherImplem.H

Go to the documentation of this file.

/* _______              __
  / ___/ /  ___  __ _  / /  ___
 / /__/ _ \/ _ \/  ' \/ _ \/ _ \
 \___/_//_/\___/_/_/_/_.__/\___/
*/
//
// This software is copyright (C) by the Lawrence Berkeley
// National Laboratory.  Permission is granted to reproduce
// this software for non-commercial purposes provided that
// this notice is left intact.
//
// It is acknowledged that the U.S. Government has rights to
// this software under Contract DE-AC03-765F00098 between
// the U.S.  Department of Energy and the University of
// California.
//
// This software is provided as a professional and academic
// contribution for joint exchange. Thus it is experimental,
// is provided ``as is'', with no warranties of any kind
// whatsoever, no support, no promise of updates, or printed
// documentation. By using this software, you acknowledge
// that the Lawrence Berkeley National Laboratory and
// Regents of the University of California shall have no
// liability with respect to the infringement of other
// copyrights by any part of this software.
//

#ifndef _GENBICGSTABSMOOTHERIMPLEM_H_
#define _GENBICGSTABSMOOTHERIMPLEM_H_

#include "LayoutIterator.H"

template <class T>
GenBiCGStabSmoother<T>::GenBiCGStabSmoother()
{
  m_maxIter = 23;
  m_solverTol = 1.0e-7;
  m_small = 1.0e-60;

  m_convergeSmall = 1.0e-2;
  // m_convergeSmall = 1.0e-15;
  // m_convergeSmall = 1.0e-60;

  m_enableRestart = true;

  m_verbose = false;
  // m_verbose = true;
}

template <class T>
GenBiCGStabSmoother<T>::~GenBiCGStabSmoother()
{
}

template <class T>
GenBiCGStabSmoother<T>* GenBiCGStabSmoother<T>::newBottomSmoother() const
{
  GenBiCGStabSmoother* newsmoother = new GenBiCGStabSmoother();
  newsmoother->setVerbose(m_verbose);

  if (newsmoother == NULL)
  {
    MayDay::Error("Out of Memory in GenBiCGStabSmoother::newBottomSmoother");
  }

  return newsmoother;
}

template <class T>
void GenBiCGStabSmoother<T>::setMaxIter(int a_maxIter)
{
  m_maxIter = a_maxIter;
}

template <class T>
void GenBiCGStabSmoother<T>::setSolverTol(Real a_solverTol)
{
  m_solverTol = a_solverTol;
}

template <class T>
void GenBiCGStabSmoother<T>::setEnableRestart(bool a_enableRestart)
{
  m_enableRestart = a_enableRestart;
}

template <class T>
void GenBiCGStabSmoother<T>::setVerbose(bool a_verbose)
{
  m_verbose = a_verbose;
}

/***********************/
// This does a GSRB Pre/Conditioned BiCGStab on a level
// for the bottom solver.
/***********************/
template <class T>
void GenBiCGStabSmoother<T>::doBottomSmooth(T&              a_phi,
                                            const T&        a_rhs,
                                            GenSolverOp<T>* a_operatorPtr)
{
  bool alreadyRestarted = false;
  // if we won't allow resarts, then act like we already have...
  if (!m_enableRestart) alreadyRestarted =true;
  Real small2 = 1.0e-8;

  T corremf;
  a_operatorPtr->defineOpData(corremf);

  T residmf;
  a_operatorPtr->defineData(residmf);

  T rtwidmf;
  a_operatorPtr->defineData(rtwidmf);

  T shatmf;
  a_operatorPtr->defineOpData(shatmf);

  T phatmf;
  a_operatorPtr->defineOpData(phatmf);

  T pmf;
  a_operatorPtr->defineData(pmf);

  T vmf;
  a_operatorPtr->defineData(vmf);

  T tmf;
  a_operatorPtr->defineData(tmf);

  a_operatorPtr->setToZero(corremf);
  a_operatorPtr->setToZero(residmf);
  a_operatorPtr->setToZero(rtwidmf);
  a_operatorPtr->setToZero(shatmf);
  a_operatorPtr->setToZero(phatmf);
  a_operatorPtr->setToZero(pmf);
  a_operatorPtr->setToZero(vmf);
  a_operatorPtr->setToZero(tmf);

  // compute initial residual
  a_operatorPtr->applyOpH(a_phi,residmf);
  a_operatorPtr->diff(residmf,a_rhs);
  a_operatorPtr->negate(residmf);

  Real initnorm;
  Real oldNorm;
  Real rnorm;
  Real norms;

  // compute residual and rhs norms
  // finish is norm(residual)/nrom(rhs) < m_solverTol
  initnorm = a_operatorPtr->norm(residmf, 0);
  rnorm = initnorm;
  oldNorm = initnorm;

  if (m_verbose)
  {
    pout() << "Bottom Smoother initNorm = " << initnorm << endl;
  }

  a_operatorPtr->copy(rtwidmf,residmf);

  Vector<Real> rho;
  Vector<Real> alpha;
  Vector<Real> omega;

  omega.push_back(0.0);
  alpha.push_back(0.0);

  // not sure if this is the right thing to do.
  // option 1: if any rhoDots are bad, add in current correction and exit
  // option 2: keep going with other components?
  bool finished = true;

  if (initnorm > m_small)
  {
    finished = false;
  }

  int iter;
  for (iter = 1; iter <= m_maxIter && !finished; iter++)
  {
    Real rhodot;

    bool badDotProd = false;
    rhodot = a_operatorPtr->dotProduct(residmf, rtwidmf);

    if (Abs(rhodot) < m_small)
    {
      badDotProd = true;
    }

    rho.push_back(rhodot);

    if (badDotProd)
    {
#if FIXIT
      for (dit.begin(); dit.ok(); ++dit)
      {
        a_phi[dit()].plus(corremf[dit()],0,0,ncomp);
      }
#endif

      if (m_verbose)
      {
        pout() << "BiCGStabSmoother: BadDotProd -- returning!  iter = "
               << iter << endl;
      }

      return;
    }

    // if we're restarting the computation, act as if it was
    // the first iteration
    if (iter == 1 || alreadyRestarted)
    {
      a_operatorPtr->copy(pmf,residmf);
    }
    else
    {
      Real betaim1 = (rho[iter-1]/rho[iter-2]) * (alpha[iter-1]/omega[iter-1]);
      Real omegaim1 = omega[iter-1];

      T temp;
      a_operatorPtr->defineData(temp);

      a_operatorPtr->copy(temp,vmf);

      a_operatorPtr->scalarMul(temp,omegaim1);
      a_operatorPtr->negate(temp);

      a_operatorPtr->sum(temp,pmf);

      a_operatorPtr->scalarMul(temp,betaim1);
      a_operatorPtr->sum(temp,residmf);

      a_operatorPtr->copy(pmf,temp);
    } // end if i ==  1 's else

    // solve M(phat) = p(i)
#if FIXIT
    a_operatorPtr->levelPreconditioner(phatmf, pmf);
#endif

    // v(i) = A(phat)
    a_operatorPtr->applyOpH(phatmf,vmf);

    bool badCorrection = false;
    Real denom = a_operatorPtr->dotProduct(vmf, rtwidmf);

    // alpha(i) = rho(i-1)/(rtwid^T v(i))
    Real alphai;

    if (Abs(denom) > small2*Abs(rho[iter-1]))
    {
      alphai = rho[iter-1]/denom;
    }
    else
    {
      alphai = 0.0;
      badCorrection = true;
    }

    alpha.push_back(alphai);

    if (!badCorrection)
    {
#if FIXIT
      // increment residual and solution
      for (dit.begin(); dit.ok(); ++dit)
      {
        T& rfab = residmf[dit()];
        T& vfab = vmf[dit()];
        T& phatfab = phatmf[dit()];
        T& corrfab = corremf[dit()];

        for (int comp = 0; comp < ncomp; comp++)
        {
          // increment residual: resid = resid - alpha*v
          rfab.plus(vfab,-1.0*alpha[comp][iter],comp,comp,1);

          // increment correction: corr = corr + alpha*phat
          corrfab.plus(phatfab,alpha[comp][iter],comp,comp,1);
        }

      } // end loop over grids
#endif
    } // end if the first correction step is OK
    else
    {
      // if correction was bad, don't increment anything,
      // but set residual to zero (forces recomputation of resid)
      a_operatorPtr->setToZero(residmf);
    }

    // compute residual norm
    norms = a_operatorPtr->norm(residmf, 0);

    // check norm of s.
    // if small enough, step to end (forces recomputation
    // of residual, etc

    bool allDone = true;
    if (norms > m_solverTol*initnorm)
    {
      allDone = false;
    }

    // if not "done" here, do second correction step
    if (!allDone)
    {

      // solve Mshat = s
#if FIXIT
      a_operatorPtr->levelPreconditioner(shatmf, residmf);
#endif

      // t = A(shat)
      a_operatorPtr->applyOpH(shatmf,tmf);

      // w(i) = tTresid/tTt
      Real numerw = a_operatorPtr->dotProduct(tmf, residmf);
      Real denomw = a_operatorPtr->dotProduct(tmf, tmf);

      Real omegai = 0.0;

      if (Abs(denomw) > m_small)
      {
        omegai = numerw/denomw;
      }

      omega.push_back(omegai);

      // corr(i) = corr(i-1) + omega(i)*shat
      // resid(i) = resid(i) - omega(i)*t

#if FIXIT
      for (dit.begin(); dit.ok(); ++dit)
      {
        Box fabbox = grids.get(dit());
        T& corrfab = corremf[dit()];
        T& shatfab = shatmf[dit()];
        T& tfab = tmf[dit()];

        T ttemp(fabbox, 1);

        for (int comp = 0; comp < ncomp; comp++)
        {
          // increment correction: corr = corr+omega*shat
          corrfab.plus(shatfab,omega[comp][iter],comp,comp,1);

          // increment residual: resid = resid - omega*t
          residmf[dit()].plus(tfab,-1.0*omega[comp][iter],comp,comp,1);
        }

      } // end loop over grids
#endif

      // reset this just in case (if we made it this
      // far, then we can assume that we're not caught in a
      // loop
      // (dfm, 5/29/02) actually, don't reset this -- only
      // allow one restart.
      // alreadyRestarted = false;
    } // end second update loop
    else
    {
      // if we skipped second update step, still need to store something
      // in omega.  try using 0
      omega.push_back(0.0);
    }

    // now check residual.
    allDone = true;

    rnorm = a_operatorPtr->norm(residmf, 0);

    if (m_verbose && (rnorm > 0))
    {
      pout() << "BottomSmoother: after " << iter
             << " iterations, residnorm = " << rnorm << endl;
    }

    if ((rnorm > m_solverTol*initnorm) && (Abs(omega[iter]) > m_small))
    {
      allDone = false;
    }

    // second test is to eliminate solver hanging
    allDone = allDone &&
              (((rnorm < m_solverTol*initnorm) && (Abs(omega[iter]) < m_small))
            || (rnorm > oldNorm*(1.0-m_convergeSmall)));

    oldNorm = rnorm;

    // if residnorm is small enough (or if correction is
    // small enough), then recompute residual from scratch
    // and check again (to avoid roundoff-drift issues).  if
    // we _still_ think we're done, then exit
    if (allDone)
    {
      // recompute residual from scratch -- first increment
      // a_phi with correction -- also reset correction to zero
      a_operatorPtr->setToZero(corremf);

#if FIXIT
      for (dit.begin(); dit.ok(); ++dit)
      {
        a_phi[dit()].plus(corremf[dit()],0,0,ncomp);
      }
#endif

      a_operatorPtr->applyOpH(a_phi,residmf);
      a_operatorPtr->diff(residmf,a_rhs);
      a_operatorPtr->negate(residmf);

      // now recompute residual norms and check for convergence
      allDone = true;

      rnorm = a_operatorPtr->norm(residmf, 0);

      if (rnorm > m_solverTol*initnorm)
      {
        allDone = false;
      }

      // if we've already restarted and haven't gotten
      // anywhere, then we've stalled and should probably
      // exit.  also, if newly-computed residual still indicates
      // convergence, then we're done.
      if (alreadyRestarted || allDone)
      {
        finished = true;
      }
      else
      {
        // otherwise, go back and try again with new residual
        if (m_verbose)
        {
          pout() << "Bottom Smoother -- restarting iterations"
                 << endl;
        }
#if FIXIT
        // also reset \tilde{r}
        for (dit.begin(); dit.ok(); ++dit)
        {
          rtwidmf[dit()].copy(residmf[dit()],0,0,ncomp);
        }
#endif

        // try doing this recursively?
        // BiCGStabSmoother newSmoother;
        // Real newTol = m_solverTol*initnorm[0]/rnorm[0];
        // newSmoother.setSolverTol(newTol);
        // newSmoother.doBottomSmooth(a_phi, a_rhs, a_operatorPtr);

        // finished = true;
        finished = false;

        alreadyRestarted = true;
      }

      if (finished && m_verbose)
      {
        pout() << "Bottom Smoother final Norm = " << rnorm << endl;
      }
    } // end if we're recomputing the residual
  } // end main loop over BiCGStab iterations

  bool converged = true;
  if (rnorm > m_solverTol*initnorm) converged = false;

  if (!converged && m_verbose)
  {
    pout() << "BiCGStabSmoother not Converged!  iter = " << iter << endl;
    pout() << " resid/initRes = " << rnorm/initnorm << endl;
  }
}

#endif

---