gsPreconditioner.h

Go to the documentation of this file.

#pragma once

#include <gsIO/gsOptionList.h>
#include <gsSolver/gsLinearOperator.h>

namespace gismo
{

template<class T>
class gsPreconditionerOp : public gsLinearOperator<T>
{
public:

    typedef memory::shared_ptr<gsPreconditionerOp> Ptr;

    typedef memory::unique_ptr<gsPreconditionerOp> uPtr;

    typedef gsLinearOperator<T> Base;

    typedef typename gsLinearOperator<T>::Ptr BasePtr;

    gsPreconditionerOp() : m_num_of_sweeps(1) {}

    virtual void step(const gsMatrix<T> & rhs, gsMatrix<T> & x) const = 0;

    virtual void stepT(const gsMatrix<T> & rhs, gsMatrix<T> & x) const
    { step(rhs, x); } // Assume symmetry.

    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    {
        x.setZero(this->rows(),input.cols()); // we assume quadratic matrices
        for ( index_t i = 0; i < m_num_of_sweeps; ++i )
            step(input,x);
    }

    virtual BasePtr underlyingOp() const = 0;

    void setNumOfSweeps( index_t n )
    {
        GISMO_ASSERT ( n > 0, "Number of sweeps needs to be positive." );
        m_num_of_sweeps = n;
    }

    index_t numOfSweeps()
    {
        return m_num_of_sweeps;
    }

    static gsOptionList defaultOptions()
    {
        gsOptionList opt = Base::defaultOptions();
        opt.addInt( "NumOfSweeps", "Number of sweeps to be applied in the member function apply", 1 );
        return opt;
    }

    virtual void setOptions(const gsOptionList & opt)
    {
        Base::setOptions(opt);
        m_num_of_sweeps = opt.askInt( "NumOfSweeps", m_num_of_sweeps );
    }

    T estimateLargestEigenvalueOfPreconditionedSystem(index_t steps = 10) const
    {
        gsMatrix<T> x, tmp;
        x.setRandom(this->rows(),1);
        for (index_t i=0; i<steps; ++i)
        {
            x.array() /= math::sqrt( x.col(0).dot(x.col(0)) );
            underlyingOp()->apply(x,tmp);
            x.setZero(this->rows(),1);
            this->step(tmp,x);
        }
        return math::sqrt( x.col(0).dot(x.col(0)) );
    }

protected:
    index_t m_num_of_sweeps;

}; // gsPreconditionerOp

template<class T>
class gsPreconditionerFromOp GISMO_FINAL : public gsPreconditionerOp<T>
{
public:

    typedef memory::shared_ptr<gsPreconditionerFromOp> Ptr;

    typedef memory::unique_ptr<gsPreconditionerFromOp> uPtr;

    typedef gsLinearOperator<T> Base;

    typedef typename gsLinearOperator<T>::Ptr BasePtr;

    gsPreconditionerFromOp( BasePtr underlying, BasePtr preconditioner, T tau = (T)1)
        : m_underlying(give(underlying)), m_preconditioner(give(preconditioner)), m_tau(tau)
    {
        GISMO_ASSERT( m_underlying->rows() == m_underlying->cols()
            && m_preconditioner->rows() == m_preconditioner->cols()
            && m_underlying->rows() == m_preconditioner->rows(),
            "The dimensions do not agree." );
    }

    static uPtr make( BasePtr underlying, BasePtr preconditioner, T tau = (T)1)
    { return uPtr( new gsPreconditionerFromOp(give(underlying), give(preconditioner), tau) ); }

    void step(const gsMatrix<T> & rhs, gsMatrix<T> & x) const
    {
        GISMO_ASSERT( m_underlying->rows() == x.rows() && x.rows() == rhs.rows() && x.cols() == rhs.cols(),
            "The dimensions do not agree." );

        m_underlying->apply(x, m_res);
        m_res -= rhs;
        m_preconditioner->apply(m_res, m_corr);
        x -= m_tau * m_corr;
    }

    // virtual void stepT(const gsMatrix<T> & rhs, gsMatrix<T> & x) const { step( rhs, x); } // Assume symmetry.

    void apply(const gsMatrix<T> & rhs, gsMatrix<T> & x) const
    {
        GISMO_ASSERT( m_underlying->rows() == rhs.rows(), "The dimensions do not agree." );

        // special treatment for the first step
        m_preconditioner->apply(rhs, x);
        if (m_tau != (T)1)
            x *= m_tau;

        for (index_t i=1; i<m_num_of_sweeps; ++i)
        {
            m_underlying->apply(x, m_res);
            m_res -= rhs;
            m_preconditioner->apply(m_res, m_corr);
            x -= m_tau * m_corr;
        }
    }

    void setDamping(const T tau) { m_tau = tau;  }

    T getDamping() const         { return m_tau; }

    static gsOptionList defaultOptions()
    {
        gsOptionList opt = Base::defaultOptions();
        opt.addReal( "Damping", "Damping parameter of the operator preconditioner", 1 );
        return opt;
    }

    void setOptions(const gsOptionList & opt)
    {
        Base::setOptions(opt);
        m_tau = opt.askReal( "Damping", m_tau );
    }

    BasePtr underlyingOp() const { return m_underlying; }
    index_t rows() const { return m_preconditioner->rows(); }
    index_t cols() const { return m_preconditioner->cols(); }

protected:
    BasePtr m_underlying;
    BasePtr m_preconditioner;
    T m_tau;
    mutable gsMatrix<T> m_res;
    mutable gsMatrix<T> m_corr;
    using gsPreconditionerOp<T>::m_num_of_sweeps;
}; // gsPreconditionerFromOp

} // namespace gismo