gsGMRes.h

Go to the documentation of this file.

#pragma once

#include <gsSolver/gsIterativeSolver.h>

namespace gismo
{

template<class T = real_t>
class gsGMRes : public gsIterativeSolver<T>
{
public:
    typedef gsIterativeSolver<T> Base;

    typedef gsMatrix<T>  VectorType;

    typedef typename Base::LinOpPtr LinOpPtr;

    typedef memory::shared_ptr<gsGMRes> Ptr;
    typedef memory::unique_ptr<gsGMRes> uPtr;

    template< typename OperatorType >
    explicit gsGMRes( const OperatorType& mat, const LinOpPtr& precond = LinOpPtr() )
    : Base(mat, precond) {}

    template< typename OperatorType >
    static uPtr make( const OperatorType& mat, const LinOpPtr& precond = LinOpPtr() )
    { return uPtr( new gsGMRes(mat, precond) ); }

    bool initIteration( const VectorType& rhs, VectorType& x );
    bool step( VectorType& x );
    void finalizeIteration( VectorType& x );

private:

    void solveUpperTriangular(const VectorType& R, const VectorType& gg)
    {
       y = R.template triangularView<gsEigen::Upper>().solve(gg);
    }

    std::ostream &print(std::ostream &os) const
    {
        os << "gsGMRes\n";
        return os;
    }

private:
    using Base::m_mat;
    using Base::m_precond;
    using Base::m_max_iters;
    using Base::m_tol;
    using Base::m_num_iter;
    using Base::m_rhs_norm;
    using Base::m_error;


    gsMatrix<T> tmp, g, g_tmp, h_tmp, y, w;
    gsMatrix<T> residual;
    gsMatrix<T> H_prev, H, Omega, Omega_prev, Omega_tmp, Omega_prev_tmp;
    std::vector< gsMatrix<T> > v;
    T beta;
};

} // namespace gismo

#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsGMRes.hpp)
#endif