eigen_direct_solver.h

Go to the documentation of this file.

/* KRATOS  _     _                       ____        _
//        | |   (_)_ __   ___  __ _ _ __/ ___|  ___ | |_   _____ _ __ ___
//        | |   | | '_ \ / _ \/ _` | '__\___ \ / _ \| \ \ / / _ \ '__/ __|
//        | |___| | | | |  __/ (_| | |   ___) | (_) | |\ V /  __/ |  \__ |
//        |_____|_|_| |_|\___|\__,_|_|  |____/ \___/|_| \_/ \___|_|  |___/ Application
//
//  Author: Thomas Oberbichler
*/
 
#if !defined(KRATOS_EIGEN_SOLVER_H_INCLUDED)
#define KRATOS_EIGEN_SOLVER_H_INCLUDED
 
// External includes
#include <Eigen/Core>
#include <Eigen/Sparse>
 
// Project includes
#include "includes/define.h"
#include "custom_utilities/ublas_wrapper.h"
#include "factories/standard_linear_solver_factory.h"
#include "includes/ublas_complex_interface.h"
#include "includes/ublas_interface.h"
#include "linear_solvers/direct_solver.h"
#include "spaces/ublas_space.h"
 
namespace Kratos {
 
template <typename Scalar>
struct SpaceType;
 
template <>
struct SpaceType<double>
{
    using Global = UblasSpace<double, CompressedMatrix, Vector>;
    using Local = UblasSpace<double, Matrix, Vector>;
};
 
template <>
struct SpaceType<std::complex<double>>
{
    using Global = UblasSpace<std::complex<double>, ComplexCompressedMatrix, ComplexVector>;
    using Local = UblasSpace<std::complex<double>, ComplexMatrix, ComplexVector>;
};
 
template <
    class TSolverType,
    class TSparseSpaceType = typename SpaceType<typename TSolverType::Scalar>::Global,
    class TDenseSpaceType = typename SpaceType<typename TSolverType::Scalar>::Local,
    class TReordererType = Reorderer<TSparseSpaceType, TDenseSpaceType>>
class EigenDirectSolver
    : public DirectSolver<TSparseSpaceType, TDenseSpaceType, TReordererType>
{
private:
    TSolverType m_solver;
 
    EigenDirectSolver &operator=(const EigenDirectSolver &Other);
 
    EigenDirectSolver(const EigenDirectSolver &Other);
 
    UblasWrapper<typename TSolverType::Scalar> m_a_wrapper;
 
public:
    KRATOS_CLASS_POINTER_DEFINITION(EigenDirectSolver);
 
    typedef DirectSolver<TSparseSpaceType, TDenseSpaceType, TReordererType> BaseType;
 
    typedef typename TSparseSpaceType::MatrixType SparseMatrixType;
 
    typedef typename TSparseSpaceType::VectorType VectorType;
 
    typedef typename TSparseSpaceType::DataType DataType;
 
    typedef typename TDenseSpaceType::MatrixType DenseMatrixType;
 
    typedef StandardLinearSolverFactory<TSparseSpaceType, TDenseSpaceType, EigenDirectSolver> FactoryType;
 
    EigenDirectSolver() {}
 
    EigenDirectSolver(Parameters settings) : BaseType(settings)
    {
        m_solver.Initialize(settings);
    }
 
    ~EigenDirectSolver() override {}
 
    void InitializeSolutionStep(SparseMatrixType& rA, VectorType& rX, VectorType& rB) override
    {
        m_a_wrapper = UblasWrapper<DataType>(rA);
 
        const auto& a = m_a_wrapper.matrix();
 
        const bool success = m_solver.Compute(a);
 
        KRATOS_ERROR_IF(!success) << "Decomposition failed!" << std::endl;
    }
 
    void PerformSolutionStep(SparseMatrixType& rA, VectorType& rX, VectorType& rB) override
    {
        Eigen::Map<Kratos::EigenDynamicVector<DataType>> x(rX.data().begin(), rX.size());
        Eigen::Map<Kratos::EigenDynamicVector<DataType>> b(rB.data().begin(), rB.size());
 
        const bool success = m_solver.Solve(b, x);
 
        KRATOS_ERROR_IF(!success) << "Solving failed!\n" << m_solver.GetSolverErrorMessages() << std::endl;
    }
 
    bool Solve(SparseMatrixType &rA, VectorType &rX, VectorType &rB) override
    {
        InitializeSolutionStep(rA, rX, rB);
        PerformSolutionStep(rA, rX, rB);
 
        return true;
    }
 
    bool Solve(SparseMatrixType &rA, DenseMatrixType &rX, DenseMatrixType &rB) override
    {
        return false;
    }
 
    void PrintInfo(std::ostream &rOStream) const override
    {
        m_solver.PrintInfo(rOStream);
    }
 
    void PrintData(std::ostream &rOStream) const override
    {
    }
 
    static StandardLinearSolverFactory<TSparseSpaceType, TDenseSpaceType, EigenDirectSolver> Factory()
    {
        return StandardLinearSolverFactory<TSparseSpaceType, TDenseSpaceType, EigenDirectSolver>();
    }
}; // class EigenDirectSolver
 
template<
    class TSolverType,
    class TSparseSpaceType,
    class TDenseSpaceType,
    class TReordererType>
inline std::istream &operator>>(
    std::istream &rIStream,
    EigenDirectSolver<TSolverType, TSparseSpaceType, TDenseSpaceType, TReordererType> &rThis)
{
    return rIStream;
}
 
template<
    class TSolverType,
    class TSparseSpaceType,
    class TDenseSpaceType,
    class TReordererType
    >
inline std::ostream &operator<<(
    std::ostream &rOStream,
    const EigenDirectSolver<TSolverType, TSparseSpaceType, TDenseSpaceType, TReordererType> &rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
} // namespace Kratos
 
#endif // defined(KRATOS_EIGEN_SOLVER_H_INCLUDED)