Krylov.hpp

/*
 * This file is part of bogus, a C++ sparse block matrix library.
 *
 * Copyright 2013 Gilles Daviet <gdaviet@gmail.com>
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
*/

#ifndef BOGUS_KRYLOV_HPP
#define BOGUS_KRYLOV_HPP

#include "BlockSolverBase.hpp"
#include "Preconditioners.hpp"
#include "KrylovMethods.hpp"

#include <vector>

namespace bogus
{



template < typename BlockMatrixType,
           template< typename BlockMatrixT > class PreconditionerType >
class Krylov : public BlockSolverBase< BlockMatrixType >
{
public:
    typedef BlockSolverBase< BlockMatrixType > Base ;
        typedef PreconditionerType< BlockObjectBase< BlockMatrixType > > PreconditionerImplType ;

    typedef typename Base::GlobalProblemTraits GlobalProblemTraits ;
    typedef typename GlobalProblemTraits::Scalar Scalar ;

    Krylov( ) ;
    explicit Krylov( const BlockObjectBase< BlockMatrixType > & matrix ) ;

    Krylov& setMatrix( const BlockObjectBase< BlockMatrixType > & matrix ) ;

    // For each value of the krylov::Method enum, create a MethodNameType typedef
    // and the asMethodName() -> MethodNameType and solve_MethodName -> Scalar methods
#define BOGUS_PROCESS_KRYLOV_METHOD( MethodName )\
    typedef krylov::solvers::MethodName<  \
        BlockMatrixType, PreconditionerType< BlockObjectBase< BlockMatrixType > >, \
        GlobalProblemTraits > MethodName##Type ; \
    \
    MethodName##Type as##MethodName() const { \
        return MethodName##Type(  \
                Base::m_matrix->derived(), \
                Base::m_maxIters, Base::m_tol, \
                &m_preconditioner, \
                &this->m_callback  \
                ) ; } \
    \
    template < typename RhsT, typename ResT > \
    Scalar solve_##MethodName( const RhsT &b, ResT &x ) const  \
    { return as##MethodName().solve( b, x ) ; }

BOGUS_KRYLOV_METHODS
#undef BOGUS_PROCESS_KRYLOV_METHOD

    template < typename RhsT, typename ResT >
    Scalar solve(  const RhsT &b, ResT &x,
                   krylov::Method method = krylov::CG ) const ;

        const PreconditionerImplType& preconditioner() const {return m_preconditioner ; }
        PreconditionerImplType& preconditioner() { return m_preconditioner ; }

protected:

    PreconditionerImplType m_preconditioner ;
} ;

} //namesoace bogus

#endif