gsIterativeSolver.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsExport.h>
#include <gsCore/gsLinearAlgebra.h>
#include <gsSolver/gsMatrixOp.h>
#include <gsIO/gsOptionList.h>

namespace gismo
{
template<class T=real_t>
class gsIterativeSolver
{
public:
    typedef memory::shared_ptr<gsIterativeSolver> Ptr;
    typedef memory::unique_ptr<gsIterativeSolver> uPtr;
    typedef T                                     ScalarType;
    typedef gsMatrix<T>                           VectorType;
    typedef typename gsLinearOperator<T>::Ptr     LinOpPtr;

    gsIterativeSolver( const LinOpPtr& mat,
                       const LinOpPtr& precond)
    : m_mat(mat),
      m_precond(precond),
      m_max_iters(1000),
      m_tol(1e-10),
      m_num_iter(-1),
      m_rhs_norm(-1),
      m_error(-1)
    {
        GISMO_ASSERT(m_mat->rows()     == m_mat->cols(),     "The matrix is not square."                     );

        if (!m_precond) m_precond = gsIdentityOp<T>::make(m_mat->rows());
        GISMO_ASSERT(m_precond->rows() == m_precond->cols(), "The preconditioner is not square."             );
        GISMO_ASSERT(m_precond->rows() == m_mat->rows(),     "The preconditioner does not match the matrix: "
                                                             <<m_precond->rows()<<"!="<<m_mat->rows()        );
    }

    template<typename Derived>
    gsIterativeSolver( const gsEigen::EigenBase<Derived> & mat,
                       const LinOpPtr& precond)
    : m_mat(makeMatrixOp(mat.derived())),
      m_precond(precond),
      m_max_iters(1000),
      m_tol(1e-10),
      m_num_iter(-1),
      m_rhs_norm(-1),
      m_error(-1)
    {
        GISMO_ASSERT(m_mat->rows()     == m_mat->cols(),     "The matrix is not square."                     );

        if (!m_precond) m_precond = gsIdentityOp<T>::make(m_mat->rows());
        GISMO_ASSERT(m_precond->rows() == m_precond->cols(), "The preconditioner is not square."             );
        GISMO_ASSERT(m_precond->rows() == m_mat->rows(),     "The preconditioner does not match the matrix: "
                                                             <<m_precond->rows()<<"!="<<m_mat->rows()        );
    }

    virtual ~gsIterativeSolver()    {}

    static gsOptionList defaultOptions()
    {
        gsOptionList opt;
        opt.addInt   ("MaxIterations"    , "Maximum number of iterations", 1000       );
        opt.addReal  ("Tolerance"        , "Tolerance for the error criteria on the "
                                           "relative residual error",      1e-10      );
        return opt;
    }

    virtual gsIterativeSolver& setOptions(const gsOptionList & opt)
    {
        m_max_iters        = opt.askInt   ("MaxIterations"    , m_max_iters        );
        m_tol              = opt.askReal  ("Tolerance"        , m_tol              );
        return *this;
    }

    void solve( const VectorType& rhs, VectorType& x )
    {
        if (initIteration(rhs, x)) return;

        while (m_num_iter < m_max_iters)
        {
            m_num_iter++;
            if (step(x)) break;
        }

        finalizeIteration(x);
    }

    void solveDetailed( const VectorType& rhs, VectorType& x, VectorType& error_history )
    {
        if (initIteration(rhs, x))
        {
            error_history.resize(1,1); //VectorType is actually gsMatrix
            error_history(0,0) = m_error;
            //gsDebug<<"Solution reached at iteration 0, err="<<m_error<<"\n";
            return;
        }

        std::vector<T> tmp_error_hist;
        tmp_error_hist.clear();
        tmp_error_hist.reserve(m_max_iters / 3 );
        tmp_error_hist.push_back(m_error); // store initial error (as provided by initIteration)

        while (m_num_iter < m_max_iters)
        {
            m_num_iter++;
            //gsDebug<<"Iteration : "<<std::setw(5)<<std::left<< m_num_iter<<"\n";

            if (step(x))
            {
                tmp_error_hist.push_back(m_error);
                //gsDebug<<"            err = "<<m_error<<" --> Solution reached.\n";
                break;
            }

            tmp_error_hist.push_back(m_error);
            //gsDebug<<"            err = "<<m_error<<"\n";
        }

        //if (m_num_iter == m_max_iters) gsDebug<<"Maximum number of iterations reached.\n";

        finalizeIteration(x);
        error_history = gsAsVector<T>(tmp_error_hist);
    }

    virtual bool initIteration( const VectorType& rhs, VectorType& x )
    {
        GISMO_ASSERT( rhs.cols() == 1,
                      "Iterative solvers only work for single column right-hand side." );
        GISMO_ASSERT( rhs.rows() == m_mat->rows(),
                      "The right-hand side does not match the matrix: "
                      << rhs.rows() <<"!="<< m_mat->rows() );

        m_num_iter = 0;

        m_rhs_norm = rhs.norm();

        if (0 == m_rhs_norm) // special case of zero rhs
        {
            x.setZero(rhs.rows(),rhs.cols()); // for sure zero is a solution
            m_error = 0.;
            return true; // iteration is finished
        }

        if ( 0 == x.size() ) // if no initial solution, start with zeros
            x.setZero(rhs.rows(), rhs.cols());
        else
        {
            GISMO_ASSERT( x.cols() == 1,
                      "Iterative solvers only work for single right-hand side and solution." );
            GISMO_ASSERT( x.rows() == m_mat->cols(),
                      "The initial guess does not match the matrix: "
                      << x.rows() <<"!="<< m_mat->cols() );
        }
        return false; // iteration is not finished
    }

    virtual bool step( VectorType& x ) = 0;                   
    virtual void finalizeIteration( VectorType& ) {}          

    index_t size() const                                       { return m_mat->rows(); }

    void setPreconditioner(const LinOpPtr & precond)           { m_precond = precond; }

    LinOpPtr preconditioner() const                            { return m_precond; }

    LinOpPtr underlying() const                                { return m_mat; }

    void setMaxIterations(index_t max_iters)                   { m_max_iters = max_iters; }

    void setTolerance(T tol)                                   { m_tol = tol; }

    index_t iterations() const                                 { return m_num_iter; }

    T error() const                                            { return m_error; }

    T tolerance() const                                        { return m_tol; }

    virtual std::ostream &print(std::ostream &os) const = 0;

    virtual std::string detail() const
    {
        std::ostringstream os;
        print(os);
        os << " Tolerance            : " << tolerance() << "\n";
        os << " Solver error         : " << error() << "\n";
        os << " Number of iterations : " << iterations() << " (max="<<m_max_iters<<")\n";
        return os.str();
    }

protected:
    const LinOpPtr     m_mat;             
    LinOpPtr           m_precond;         
    index_t            m_max_iters;       
    gsOptionList::Real m_tol;             
    index_t            m_num_iter;        
    T                  m_rhs_norm;        
    T                  m_error;           
};

template<class T>
std::ostream &operator<<(std::ostream &os, const gsIterativeSolver<T>& b)
{return b.print(os); }

template <class SolverType>
class gsIterativeSolverOp GISMO_FINAL : public gsLinearOperator<typename SolverType::ScalarType>
{
    typedef typename SolverType::ScalarType T;
    typedef typename gsLinearOperator<T>::Ptr LinOpPtr;
public:
    typedef memory::shared_ptr<gsIterativeSolverOp> Ptr;

    typedef memory::unique_ptr<gsIterativeSolverOp> uPtr;

    template<class MatrixType>
    gsIterativeSolverOp(const MatrixType& matrix, const LinOpPtr& preconditioner = LinOpPtr())
        : m_solver(matrix, preconditioner) {}

    template<class MatrixType>
    static uPtr make(const MatrixType& matrix, const LinOpPtr& preconditioner = LinOpPtr())
    { return uPtr( new gsIterativeSolverOp(matrix,preconditioner) ); }

    static uPtr make(SolverType solver) { return memory::make_unique( new gsIterativeSolverOp(solver) ); }

    void apply(const gsMatrix<T> & input, gsMatrix<T> & res) const
    {
        res.setZero(input.rows(), input.cols());
        m_solver.solve(input, res);
        gsDebug << ( (m_solver.error() < m_solver.tolerance())
                     ? "gsIterativeSolverOp reached desired error bound after "
                     : "msIterativeSolverOp did not reach desired error bound within " )
                << m_solver.iterations() << " iterations.\n";
        gsDebug << input.rows() << "; " << input.norm() << "; " << res.norm() << "\n";
    }

    index_t rows() const { return m_solver.underlying()->rows(); }

    index_t cols() const { return m_solver.underlying()->cols(); }

    SolverType& solver()                { return m_solver; }

    const SolverType& solver() const    { return m_solver; }

private:
    mutable SolverType m_solver;  // The iterative solvers change their state during solving
};


} // namespace gismo