diagonal_preconditioner.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//
//
 
#if !defined(KRATOS_DIAGONAL_PRECONDITIONER_H_INCLUDED )
#define  KRATOS_DIAGONAL_PRECONDITIONER_H_INCLUDED
 
// System includes
 
 
// External includes
#include <boost/numeric/ublas/vector.hpp>
 
 
// Project includes
#include "includes/define.h"
 
 
namespace Kratos
{
 
 
 
 
 
 
 
 
template<class TSparseSpaceType, class TDenseSpaceType>
class DiagonalPreconditioner : public Preconditioner<TSparseSpaceType, TDenseSpaceType>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(DiagonalPreconditioner);
 
    typedef  Preconditioner<TSparseSpaceType, TDenseSpaceType> BaseType;
 
    typedef typename TSparseSpaceType::DataType DataType;
 
    typedef typename TSparseSpaceType::MatrixType SparseMatrixType;
 
    typedef typename TSparseSpaceType::VectorType VectorType;
 
    typedef typename TDenseSpaceType::MatrixType DenseMatrixType;
 
 
    DiagonalPreconditioner() {}
 
    DiagonalPreconditioner(const DiagonalPreconditioner& Other)
        : BaseType(Other), mDiagonal(Other.mDiagonal), mTemp(Other.mTemp) {}
 
    ~DiagonalPreconditioner() override {}
 
 
 
    DiagonalPreconditioner& operator=(const DiagonalPreconditioner& Other)
    {
        BaseType::operator=(Other);
        mDiagonal = Other.mDiagonal;
        mTemp = Other.mTemp;
        return *this;
    }
 
 
 
 
    void Initialize(SparseMatrixType& rA, VectorType& rX, VectorType& rB) override
    {
        mDiagonal.resize(TSparseSpaceType::Size(rX));
        mTemp.resize(TSparseSpaceType::Size(rX));
 
        const DataType zero = DataType();
 
        IndexPartition<std::size_t>(rA.size1()).for_each([&](std::size_t Index){
            double diag_Aii = rA(Index,Index);
            if(diag_Aii != zero)
                mDiagonal[Index] = 1.00 / std::sqrt(std::abs(diag_Aii));
            else
                KRATOS_THROW_ERROR(std::logic_error,"zero found in the diagonal. Diagonal preconditioner can not be used","");
        });
 
    }
 
    void Initialize(SparseMatrixType& rA, DenseMatrixType& rX, DenseMatrixType& rB) override
    {
        BaseType::Initialize(rA, rX, rB);
    }
 
    void Mult(SparseMatrixType& rA, VectorType& rX, VectorType& rY) override
    {
        IndexPartition<std::size_t>(TSparseSpaceType::Size(rX)).for_each([&](std::size_t Index){
            mTemp[Index] = rX[Index] * mDiagonal[Index];
        });
 
        TSparseSpaceType::Mult(rA,mTemp, rY);
        ApplyLeft(rY);
    }
 
    void TransposeMult(SparseMatrixType& rA, VectorType& rX, VectorType& rY) override
    {
        IndexPartition<std::size_t>(TSparseSpaceType::Size(rX)).for_each([&](std::size_t Index){
            mTemp[Index] = rX[Index] * mDiagonal[Index];
        });
 
        TSparseSpaceType::TransposeMult(rA,mTemp, rY);
        ApplyRight(rY);
    }
 
    VectorType& ApplyLeft(VectorType& rX) override
    {
        IndexPartition<std::size_t>(TSparseSpaceType::Size(rX)).for_each([&](std::size_t Index){
            rX[Index] *= mDiagonal[Index];
        });
 
        return rX;
    }
 
    VectorType& ApplyRight(VectorType& rX) override
    {
        IndexPartition<std::size_t>(TSparseSpaceType::Size(rX)).for_each([&](std::size_t Index){
            rX[Index] *= mDiagonal[Index];
        });
 
        return rX;
    }
 
    VectorType& ApplyTransposeLeft(VectorType& rX) override
    {
        IndexPartition<std::size_t>(TSparseSpaceType::Size(rX)).for_each([&](std::size_t Index){
            rX[Index] *= mDiagonal[Index];
        });
 
        return rX;
    }
 
    VectorType& ApplyTransposeRight(VectorType& rX) override
    {
        IndexPartition<std::size_t>(TSparseSpaceType::Size(rX)).for_each([&](std::size_t Index){
            rX[Index] *= mDiagonal[Index];
        });
 
        return rX;
    }
 
    VectorType& ApplyInverseRight(VectorType& rX) override
    {
        IndexPartition<std::size_t>(TSparseSpaceType::Size(rX)).for_each([&](std::size_t Index){
            rX[Index] /= mDiagonal[Index];
        });
 
        return rX;
    }
 
    VectorType& Finalize(VectorType& rX) override
    {
        IndexPartition<std::size_t>(TSparseSpaceType::Size(rX)).for_each([&](std::size_t Index){
            rX[Index] *= mDiagonal[Index];
        });
 
        return rX;
    }
 
 
 
 
 
    std::string  Info() const override
    {
        return "Diagonal preconditioner";
    }
 
    void  PrintInfo(std::ostream& OStream) const override
    {
        OStream << "Diagonal preconditioner";
    }
 
 
 
 
 
protected:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
private:
 
 
 
    VectorType mDiagonal;
 
    VectorType mTemp;
 
 
 
 
 
 
 
 
 
 
 
 
 
}; // Class DiagonalPreconditioner
 
 
 
 
 
 
 
template<class TSparseSpaceType, class TDenseSpaceType>
inline std::istream& operator >> (std::istream& IStream,
                                  DiagonalPreconditioner<TSparseSpaceType, TDenseSpaceType>& rThis)
{
    return IStream;
}
 
template<class TSparseSpaceType, class TDenseSpaceType>
inline std::ostream& operator << (std::ostream& OStream,
                                  const DiagonalPreconditioner<TSparseSpaceType, TDenseSpaceType>& rThis)
{
    rThis.PrintInfo(OStream);
    OStream << std::endl;
    rThis.PrintData(OStream);
 
    return OStream;
}
 
 
}  // namespace Kratos.
 
#endif // KRATOS_DIAGONAL_PRECONDITIONER_H_INCLUDED  defined