diagprecond.h

Go to the documentation of this file.

#ifndef TBCI_SOLVER_DIAGPRECOND_H
#define TBCI_SOLVER_DIAGPRECOND_H  

#include "tbci/solver/precond.h"

NAMESPACE_TBCI

template <typename T, typename MatrixType = BdMatrix<T> >
class DiagPreconditioner : public Preconditioner_Sig<T, MatrixType>
{
 private:
        BVector<T> invDiag;
        bool checkdiag;
 public:
        friend class SystemMatrix;
        
        DiagPreconditioner()
        {}
        DiagPreconditioner(const MatrixType &A)
        { update (A); }
        ~DiagPreconditioner()  {}
        
        // update-Funktion
        void update(const MatrixType &A)
        {
                const unsigned int rw = A.rows();
                if (invDiag.size() != A.rows()) 
                        invDiag.resize(A.rows());
                BCHK(invDiag.size() == 0, VecErr, "DiagPrecond::update: NULL ptr", 0,);
                if (do_expcheck()) {
                        /* Check for diagonal zeros separately */
                        for (unsigned int j = 0; j < rw; ++j)
                                EXPCHK(A.get(j,j) == T(0.0), VecErr, "DiagPrecond::update: Diagonal ZERO", j,);
                }
                for (unsigned int i = 0; i < rw; ++i)
                        invDiag.set(i) = T(1.0) / A.get(i,i);
        }
        
        
        // Solve-Funktionen
        /* 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> solve(TVector<T> tv) const
        {
                const unsigned long sz = invDiag.size();
                //const unsigned long sz = tv.size();
                // FIXME: This is do_vec_mul_vec!
                do_vec_mul_vec<T> (sz, &tv.setval(0), &invDiag.getcref(0));
                //for (unsigned int i = 0; i < sz; ++i)
                //      tv.set(tv.get(i) * invDiag.get(i), i);
                return 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); } 

};

NAMESPACE_END

#endif  /* TBCI_SOLVER_DIAGPRECOND_H */