diluprecond.h

Go to the documentation of this file.

#ifndef TBCI_SOLVER_DILUPRECOND_H
#define TBCI_SOLVER_DILUPRECOND_H  

#include "tbci/solver/precond.h"

NAMESPACE_TBCI

// ********************************************************************************
// ... For general BdMatrices
// ********************************************************************************


template <typename T> class BdMatrix;

template <typename T>
class DILU_BdMatrixPreconditioner : public Preconditioner_Sig<T, BdMatrix<T> >
{
    public:
        DILU_BdMatrixPreconditioner ()  {}
        DILU_BdMatrixPreconditioner (const BdMatrix<T> &A)
                { update(A); }
        ~DILU_BdMatrixPreconditioner ()  {}

        // update function
        void update (const BdMatrix<T> &A);

        // solve functions
        TVector<T> solve (TVector<T> x) const;

        inline TVector<T> solve (const Vector<T> &v)  const 
        {
                TVector<T> tv(v);
                return (solve(tv));
        } 

        /* template <typename T> */
        inline TVector<T> transSolve (const Vector<T> &v) const  
                { return solve(v); } 
        /* template <typename T> */
        inline TVector<T> transSolve ( TVector<T> tv) const  
                { return solve(tv); }

    private:
        BVector<T> pivots;
        unsigned int dimension;
};

// See Templates Book, Fig. 3.3
template <typename T>
void DILU_BdMatrixPreconditioner<T>::update (const BdMatrix<T> &A) 
{
        // copy matrix variables to private vars
        dimension  = A.dim;

        // calculate Pivots
        if (pivots.size() != dimension) 
                pivots.resize(dimension);
        BCHK(pivots.size() == 0, VecErr, "Null ptr in DILUBdMPrecCond::update", 0, );
        const unsigned ds = A.diagconf.size();
        unsigned i; 
        for (i = 0; i < dimension; ++i)
                pivots.set(i) = A.diag[i];
        for (i = 0; i < dimension; ++i) {
                REGISTER T tmp = pivots.get(i);
                if (UNLIKELY(tmp == T(0.0)))
                        tmp = T(1.0);
                else
                        tmp = T(1.0)/tmp;
                pivots.set(i) = tmp;
                for (unsigned j = 0; j < ds; j++) {
                        const int df = A.diagconf.get(j);
                        if (LIKELY(df < (signed)dimension - (signed)i))
                                pivots.set(i+df) -= A.adiag.get(df)[i]*tmp*A.bdiag.get(df)[i];
                }
        }
}

template <typename T>
/*inline*/ TVector<T> DILU_BdMatrixPreconditioner<T>::solve (TVector<T> y) const 
{

        // check dimensions
        BCHKNR(dimension != y.size(), VecErr, "DILU_BdMPreCond::solve: Dimension mismatch", y.size());
        
        // check trivial solution
        // Matrix consists of main diagonal only
        for (unsigned int i = 0; i < dimension; ++i)
                y.setval(i) *= pivots.get(i);
        return y;
}

NAMESPACE_END
                
#endif  /* TBCI_SOLVER_DILUPRECOND_H */