ilu0precond.h

Go to the documentation of this file.

#ifndef TBCI_SOLVER_ILU0PRECOND_H
#define TBCI_SOLVER_ILU0PRECOND_H  

#include "tbci/solver/precond.h"

NAMESPACE_TBCI

template <typename T> class Symm_BdMatrix;

template <typename T>
class ILU0_Symm_BdMatrixPreconditioner : public Preconditioner_Sig<T, Symm_BdMatrix<T> >
{
    public:
        ILU0_Symm_BdMatrixPreconditioner () : matrix(0)  {}
        ILU0_Symm_BdMatrixPreconditioner (const Symm_BdMatrix<T> &A)  
                { update (A); }
        ~ILU0_Symm_BdMatrixPreconditioner ()  {}
        
        
        // update function
        void update (const Symm_BdMatrix<T> &A);
        //void update (const MatrixType &A);
        
        // solve functions
        TVector<T> solve (TVector<T> x) const;
        
        /* template <typename T> */
        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;
        const Symm_BdMatrix<T> *matrix;

        // copy Matrix varibales to local vars
        BVector<unsigned int> conf;
        unsigned int dimension;
        BVector<T*> rowPtr;
};

template<typename T>
void ILU0_Symm_BdMatrixPreconditioner<T>::update (const Symm_BdMatrix<T> &A)
{
        matrix = &A;
        if (pivots.size() != matrix->rows()) 
                pivots.resize(matrix->rows());  
        if (UNLIKELY(pivots.size() == 0))
                STD__ cerr << "\n NULL pointer encountered! \n" << STD__ flush;
        for (unsigned int i=0; i<matrix->rows(); ++i)
                pivots(i) = T(1.0)/(*matrix)(i,i);
        // copy Matrix varibales to local vars
        conf.resize(matrix->conf);
        dimension = matrix->dimension;
        rowPtr.resize(matrix->rowPtr);
}


template<typename T>
/*inline*/ TVector<T> ILU0_Symm_BdMatrixPreconditioner<T>::solve (TVector<T> y) const
{
        unsigned int* confPtr;
        T* elementPtr;
        ALIGN2(REGISTER T sum, MIN_ALIGN2);

        const unsigned int csize = conf.size();
        // check for trivial solution
        // Matrix consists of just the main diagonal
        if (LIKELY(csize == 0)) {
                for (unsigned int i=0; i<dimension; ++i)
                        y.setval(i) = pivots(i)*y.get(i);
                return y;
        }               
        
        
        // Otherwise solve via incompl. LU decomp.:
        // LUy=x, Uy=z, L=D+L, U=I+D^(-1)*U
        // We only generate one TVector, which will be written to twice

        // Lz=y-> (D + L)*z = y
        // write to y instead of temporary vector (z)
        for (int g = 0; g < (int)dimension; ++g) {
                sum = (T)0.0; 
                int cnfi;
                for (unsigned int i=0; (i<csize) && ((cnfi = conf(i)) <= g);) {
                        elementPtr = rowPtr(g - cnfi) + ++i;
                        sum += (*elementPtr) * y.get(g - cnfi);
                }

                y.setval(g) = pivots(g) * (y.get(g) - sum);
        }
                        

        // Ux=z -> (I + D^(-1)*U)*x = z
        // take the values obtained above for the RHS (z) 
        // and overwrite step by step 
        for (int h = dimension-2; h >= 0; --h)
        {
                confPtr = conf.vecptr();
                elementPtr = rowPtr(h) +1;
                                
                sum = (T)0.0;
                while (elementPtr != rowPtr(h+1)) {
                        sum += (*elementPtr) * y.get(h + (*confPtr));
                        confPtr++;
                        elementPtr++;
                }
        
                y.setval(h) -= pivots(h) * sum;
        }

        // presto!
        return y;
}


// ********************************************************************************
// Now for general BdMatrices
// ********************************************************************************


template <typename T> class BdMatrix;

template <typename T>
class ILU0_BdMatrixPreconditioner : public Preconditioner_Sig<T, BdMatrix<T> >
{
    public:
        ILU0_BdMatrixPreconditioner ()  {}
        ILU0_BdMatrixPreconditioner (const BdMatrix<T> &A)
                { update(A); }
        ~ILU0_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;

        BVector<T*> ldiag, udiag;
        BVector<unsigned int> bdconf;

        unsigned int dimension;
};

template <typename T>
void ILU0_BdMatrixPreconditioner<T>::update (const BdMatrix<T> &A) 
{
        //copy matrix variables to private vars
        dimension = A.dim;
        udiag.resize(A.adiag); ldiag.resize(A.bdiag);

        // calculate Pivots
        if (pivots.size() != dimension) 
                pivots.resize(dimension);
        BCHK(pivots.size() == 0, VecErr, "ILU0_BdMPrecond::update: Null ptr", 0,);
        for (unsigned i = 0; i < dimension; ++i)
                pivots.set(i) = T(1.0) / A.get(i,i);

        // fill in configuration vectors for BandMatrix
        bdconf.resize(A.diagconf);
}

template <typename T>
TVector<T> ILU0_BdMatrixPreconditioner<T>::solve(TVector<T> y) const
{
        ALIGN2(REGISTER T sum, MIN_ALIGN2);

        // check dimensions
        BCHK(dimension != y.size(), VecErr, "ILU0_BdMPrecond::solve: Dims don't match", dimension, y);
        // check trivial solution
        // Matrix consists of main diagonal only
        if (bdconf.size() == 0) {
                for (unsigned int i=0; i<dimension; ++i)
                        y.setval(i) *= pivots.get(i);
                return y;
        }

        // Otherwise solve via incompl. LU decomposition:
        // LUy=x, Uy=z, L=D+L, U=I+D^(-1)*U
        // We only use the passed TVector, which will be written to twice

        // Lz=y-> (D + L)*z = y
        // write to y instead of temporary Vector (z)
        y.setval(0) = pivots.get(0) * (y.get(0));
        int h;
        const unsigned int bsize = bdconf.size();
        for (h = 1; h < (int)dimension; ++h) {
                sum = (T)0; int bdi;
                for (unsigned int i=0; i < bsize && (bdi = bdconf.get(i)) <= h; ++i)
                        sum += (*(ldiag.get(bdi) + h-bdi)) * y.get (h-bdi);
                y.setval(h) = pivots.get(h) * (y.get(h) - sum);
        }

        // Ux=z -> (I + D^(-1)*U)*x = z
        // for the RHS (z), take the values obtained above
        // and overwrite y piecewise

        //y.setval(dimension-1) = y.get(dimension-1);
        for (h = dimension-2; h >= 0; --h) {
                sum = (T)0; int bdi;
                for (unsigned int i=0; i < bsize && (bdi = bdconf.get(i))+h < (int)dimension; ++i)
                        sum += (*(udiag.get(bdi) + h)) * y.get(bdi + h);
                y.setval(h) -= pivots.get(h) * sum;
        }

        // presto!
        return y;
}

NAMESPACE_END
                
#endif  /* TBCI_SOLVER_ILU0PRECOND_H */