GaussSeidelBase.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_BASE_HPP
#define BOGUS_BLOCK_GAUSS_SEIDEL_BASE_HPP

#include "ConstrainedSolverBase.hpp"
#include "Coloring.hpp"

#include <vector>

namespace bogus
{

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

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


    void setMaxThreads( unsigned maxThreads = 0 ) {
        m_maxThreads = maxThreads ;
    }



    void setAutoRegularization( Scalar maxRegul ) { m_autoRegularization = maxRegul ; }

    // Debug


    void setEvalEvery( unsigned evalEvery ) { m_evalEvery = evalEvery ; }
    void setSkipTol  ( Scalar skipTol   ) { m_skipTol   = skipTol   ; }
    void setSkipIters( unsigned skipIters ) { m_skipIters = skipIters ; }


protected:

    GaussSeidelBase() :
        m_maxThreads (  0 ),
        m_evalEvery ( 25  ),
        m_skipIters ( 10 ),
        m_autoRegularization ( 0. )
    {
        m_tol = 1.e-6 ;
        m_maxIters = 250 ;
        m_skipTol  = 1.e-6 ;
    }

    const BlockMatrixBase< BlockMatrixType >& explicitMatrix() const
    {
        // Will cause a compile error if BlockMatrixType does not derive from BlockMatrixBase
        // This is voluntary, GaussSeidel implementations does not handle arbitrary expressions yet
        return m_matrix->derived() ;
    }
    const IterableBlockObject< BlockMatrixType >& iterableMatrix() const
    {
        // Will cause a compile error if BlockMatrixType does not derive from IterableObjectBase
        // This is voluntary, GaussSeidel implementations does not handle arbitrary expressions yet
        return m_matrix->derived() ;
    }

    // Add regularization, compute scalings from diagonal block matrices
    void processLocalMatrices() ;

    template < typename NSLaw, typename ResT >
    Scalar evalAndKeepBest(
            const NSLaw &law, const ResT &x,
            const typename GlobalProblemTraits::DynVector& y,
            typename GlobalProblemTraits::DynVector& x_best, Scalar &err_best ) const ;

    template < typename NSLaw, typename RhsT, typename ResT >
    bool tryZero(
            const NSLaw &law, const RhsT &b, ResT &x,
            typename GlobalProblemTraits::DynVector& x_best, Scalar &err_best ) const ;


    using Base::m_matrix ;
    using Base::m_maxIters ;
    using Base::m_tol ;
    using Base::m_scaling ;

    typedef typename Base::Index Index ;
    typedef typename Base::BlockProblemTraits::Matrix DiagonalBlockType ;
    typename ResizableSequenceContainer< DiagonalBlockType >::Type m_localMatrices ;
    typename GlobalProblemTraits::DynVector m_regularization ;

    unsigned m_maxThreads ;

    unsigned m_evalEvery ;
    Scalar m_skipTol ;
    unsigned m_skipIters ;

    Scalar m_autoRegularization ;

} ;

} //namespace bogus


#endif