amgstep.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef ISTL_AMG_STEP_HH
#define ISTL_AMG_STEP_HH

/** \file 
    \brief A wrapper class for the ISTL AMG implementation
 */

#include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include <dune/istl/paamg/amg.hh>

/** \brief A wrapper class for the ISTL AMG implementation
 */
template <class MatrixType, class VectorType>
class AMGStep
    : public LinearIterationStep<MatrixType, VectorType>
{
    typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;

    typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType,Dune::Amg::RowSum> > Criterion;

    /** \brief Use a sequential SSOR for smoothing */
    typedef Dune::SeqSSOR<MatrixType,VectorType,VectorType> Smoother;
    typedef typename Dune::Amg::template SmootherTraits<Smoother>::Arguments SmootherArgs;
    
    typedef Dune::Amg::AMG<Operator,VectorType,Smoother> AMG;

public:

    /** \brief Default constructor */
    AMGStep () {}

    /** \brief Constructor which initializes and sets up an algebraic hierarchy
        \param smootherArgs  Arguments for the smoother.  See the dune-istl documentation for details
        \param coarseningCriterion  Arguments for the coarsening.  See the dune-istl documentation for details
     */
    AMGStep (const MatrixType* stiffnessMatrix, 
             VectorType& x, 
             VectorType& rhs,
             const SmootherArgs& smootherArgs,
             const Criterion& coarseningCriterion)
    {
        setProblem(*stiffnessMatrix, x, rhs);

        fop_ = std::auto_ptr<Operator>(new Operator(*stiffnessMatrix));
        amg_ = std::auto_ptr<AMG>(new AMG(*fop_, coarseningCriterion, smootherArgs, 1, 1, 1, false));
        amg_->pre(*this->x_,*this->rhs_);
    }

    /** \brief Initialize the iteration step but don't actually build the matrix hierarchy yet */
    AMGStep (const MatrixType* stiffnessMatrix, 
             VectorType& x, 
             VectorType& rhs)
    {
        setProblem(*stiffnessMatrix, x, rhs);
    }

    /** \brief Initialize the iteration step but don't actually build the matrix hierarchy yet */
    virtual void setProblem(const MatrixType& stiffnessMatrix, 
             VectorType& x, 
             VectorType& rhs)
    {
        LinearIterationStep<MatrixType, VectorType>::setProblem(stiffnessMatrix, x, rhs);

        fop_ = std::auto_ptr<Operator>(new Operator(stiffnessMatrix));
    }

    /** \brief Sets up an algebraic hierarchy
     */
    virtual void preprocess();

    /** \brief Perform one iteration */
    virtual void iterate();

    /** \brief Get the solution */
    virtual VectorType getSol();

    /** \brief Arguments for the smoother.  See the dune-istl documentation for details */
    SmootherArgs smootherArgs_;

    /** \brief Arguments for the coarsening.  See the dune-istl documentation for details */
    Criterion coarseningCriterion_;

private:

    /** \brief A matrix wrapper demanded by istl */
    std::auto_ptr<Operator> fop_;

    /** \brief The dune-istl AMG step */
    std::auto_ptr<AMG> amg_;
};

template <class MatrixType, class VectorType>
void AMGStep<MatrixType,VectorType>::preprocess()
{
    amg_ = std::auto_ptr<AMG>(new AMG(*fop_, coarseningCriterion_, smootherArgs_, 1, 1, 1, false));
    amg_->pre(*this->x_,*this->rhs_);
}

template <class MatrixType, class VectorType>
void AMGStep<MatrixType,VectorType>::iterate()
{
    amg_->apply(*this->x_,*this->rhs_);
}

template <class MatrixType, class VectorType>
VectorType AMGStep<MatrixType,VectorType>::getSol()
{
    return *this->x_;
}

#endif