ProductGaussSeidel.hpp

/*
 * This file is part of bogus, a C++ sparse block matrix library.
 *
 * Copyright 2016 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_PRODUCT_GAUSS_SEIDEL_HPP
#define BOGUS_PRODUCT_GAUSS_SEIDEL_HPP

#include "GaussSeidelBase.hpp"
#include "ProductGaussSeidelUtils.hpp"

#include <vector>

namespace bogus
{


template < typename BlockMatrixType, typename DiagonalType = typename BlockMatrixType::Scalar,
           bool PrecomputeDMt = !( IsSame<DiagonalType, typename BlockMatrixType::Scalar>::Value ) >
class ProductGaussSeidel
        : public GaussSeidelBase< ProductGaussSeidel<BlockMatrixType, DiagonalType, PrecomputeDMt>, BlockMatrixType >
{
public:
    typedef GaussSeidelBase< ProductGaussSeidel, BlockMatrixType > Base ;

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

    enum {has_trivial_diagonal = IsSame<DiagonalType, typename BlockMatrixType::Scalar>::Value } ;
    typedef block_solvers_impl::DiagonalMatrixWrapper<DiagonalType, !!has_trivial_diagonal > DiagWrapper ;
    typedef block_solvers_impl::DMtStorage<BlockMatrixType, DiagWrapper, PrecomputeDMt> DMtStorage ;

    ProductGaussSeidel()
        : Base()
    { }
    explicit ProductGaussSeidel( const BlockObjectBase< BlockMatrixType > & matrix )
        : Base()
    {  setMatrix( matrix ) ; }
    ProductGaussSeidel( const BlockObjectBase< BlockMatrixType > & matrix, const DiagonalType& diagonal )
        : Base(), m_diagonal( DiagWrapper( diagonal ) )
    {  setMatrix( matrix ) ; }

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

    ProductGaussSeidel& setDiagonal( const DiagonalType &diagonal ) ;

    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 ;

protected:

    DiagWrapper m_diagonal ;
    DMtStorage  m_DMt ;

    void updateLocalMatrices();

    template < typename NSLaw,  typename VecT, typename ResT >
    void innerLoop (
        bool parallelize, const NSLaw &law, const VecT& b,
        std::vector< unsigned char > &skip, Scalar &ndxRef,
        VecT& Mx, 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 ;

} ;

} //namespace bogus


#endif