GaussSeidel.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_BLOCK_GAUSS_SEIDEL_HPP
#define BOGUS_BLOCK_GAUSS_SEIDEL_HPP

#include "GaussSeidelBase.hpp"
#include "Coloring.hpp"

#include <vector>

namespace bogus
{


template < typename BlockMatrixType >
class GaussSeidel : public GaussSeidelBase< GaussSeidel<BlockMatrixType>, BlockMatrixType >
{
public:
    typedef GaussSeidelBase< GaussSeidel, BlockMatrixType > Base ;

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

    GaussSeidel( ) : Base() { }
    explicit GaussSeidel( const BlockObjectBase< BlockMatrixType > & matrix ) : Base()
    {  setMatrix( matrix ) ; }

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


    template < typename NSLaw, typename RhsT, typename ResT >
    Scalar solve( const NSLaw &law, const RhsT &b, ResT &x, bool tryZeroAsWell = true ) const ;

    template < typename NSLaw, typename RhsT, typename ResT, typename LSDerived, typename HDerived >
    Scalar solveWithLinearConstraints( const NSLaw &law,
                                       const BlockObjectBase< LSDerived >& Cinv,
                                       const BlockObjectBase<  HDerived >& H,
                                       const Scalar alpha,
                                       const RhsT &b, const RhsT &c,
                                       ResT &x,
                                       bool tryZeroAsWell = true, unsigned solveEvery = 1) const ;

    template < typename NSLaw, typename RhsT, typename ResT, typename WDerived >
    Scalar solveWithLinearConstraints( const NSLaw &law,
                                       const BlockObjectBase < WDerived >& W,
                                       const RhsT &b, ResT &x,
                                       bool tryZeroAsWell = true, unsigned solveEvery = 1 ) const ;


    Coloring& coloring( ) { return m_coloring ; }

protected:
    void updateLocalMatrices() ;

    template < typename NSLaw,  typename RhsT, typename ResT >
    void innerLoop (
        bool parallelize, const NSLaw &law, const RhsT& b,
        std::vector< unsigned char > &skip, Scalar &ndxRef,
        ResT &x ) const ;

    typedef typename Base::Index Index ;

    using Base::m_matrix ;
    using Base::m_maxIters ;
    using Base::m_tol ;
    using Base::m_scaling ;
    using Base::m_maxThreads ;
    using Base::m_evalEvery ;
    using Base::m_skipTol ;
    using Base::m_skipIters ;
    using Base::m_localMatrices ;
    using Base::m_regularization ;

    Coloring m_coloring ;
} ;

} //namespace bogus


#endif