amesos_solver.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//  Collaborator:    Philipp Bucher
//
 
#if !defined (KRATOS_AMESOS_SOLVER_H_INCLUDED)
#define KRATOS_AMESOS_SOLVER_H_INCLUDED
 
// External includes
 
// Project includes
#include "includes/define.h"
#include "linear_solvers/linear_solver.h"
 
// Amesos solver includes
#include "Amesos.h"
#include "Epetra_LinearProblem.h"
 
namespace Kratos
{
 
 
template< class TSparseSpaceType, class TDenseSpaceType,
          class TReordererType = Reorderer<TSparseSpaceType, TDenseSpaceType> >
class AmesosSolver : public LinearSolver< TSparseSpaceType,
    TDenseSpaceType, TReordererType>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(AmesosSolver);
 
    typedef typename TSparseSpaceType::MatrixType SparseMatrixType;
 
    typedef typename TSparseSpaceType::VectorType VectorType;
 
    typedef typename TDenseSpaceType::MatrixType DenseMatrixType;
 
 
    AmesosSolver(Parameters settings)
    {
        Parameters default_settings( R"({
            "solver_type"                    : "amesos",
            "amesos_solver_type"             : "Amesos_Klu",
            "trilinos_amesos_parameter_list" : { }
        }  )" );
 
        // choose solver-type
        if (settings["solver_type"].GetString() == "klu") {
            if (settings.Has("amesos_solver_type")) {
                KRATOS_INFO("Amesos-Solver") << "Ignoring setting \"amesos_solver_type\"" << std::endl;
            } else {
                settings.AddEmptyValue("amesos_solver_type");
            }
            settings["amesos_solver_type"].SetString("Amesos_Klu");
        }
        else if (settings["solver_type"].GetString() == "super_lu_dist") {
            if (settings.Has("amesos_solver_type")) {
                KRATOS_INFO("Amesos-Solver") << "Ignoring setting \"amesos_solver_type\"" << std::endl;
            }
            else {
                settings.AddEmptyValue("amesos_solver_type");
            }
            settings["amesos_solver_type"].SetString("Amesos_Superludist");
        }
        else if (settings["solver_type"].GetString() == "mumps") {
            if (settings.Has("amesos_solver_type")) {
                KRATOS_INFO("Amesos-Solver") << "Ignoring setting \"amesos_solver_type\"" << std::endl;
            }
            else {
                settings.AddEmptyValue("amesos_solver_type");
            }
            settings["amesos_solver_type"].SetString("Amesos_Mumps");
        }
        else if (settings["solver_type"].GetString() == "amesos") {
            // do nothing here. Leave full control to the user through the "trilinos_amesos_parameter_list"
            // and the "amesos_solver_type"
        }
        else {
            KRATOS_ERROR << "The solver type specified: \"" << settings["solver_type"].GetString() << "\" is not supported";
        }
 
        settings.ValidateAndAssignDefaults(default_settings);
 
        //assign the amesos parameter list, which may contain parameters IN TRILINOS INTERNAL FORMAT to mParameterList
        mParameterList = Teuchos::ParameterList();
        for(auto it = settings["trilinos_amesos_parameter_list"].begin(); it != settings["trilinos_amesos_parameter_list"].end(); it++) {
            if(it->IsString()) mParameterList.set(it.name(), it->GetString());
            else if(it->IsInt()) mParameterList.set(it.name(), it->GetInt());
            else if(it->IsBool()) mParameterList.set(it.name(), it->GetBool());
            else if(it->IsDouble()) mParameterList.set(it.name(), it->GetDouble());
        }
 
        mSolverName = settings["amesos_solver_type"].GetString();
 
        KRATOS_ERROR_IF_NOT(HasSolver(mSolverName)) << "attempting to use Amesos solver \"" << mSolverName
            << "\" unfortunately the current compilation of Trilinos does not include it" << std::endl;
    }
 
    AmesosSolver(const std::string& SolverName, Teuchos::ParameterList& rParameterList)
    {
        mParameterList = rParameterList;
        mSolverName = SolverName;
 
        KRATOS_ERROR_IF_NOT(HasSolver(mSolverName)) << "attempting to use Amesos solver \"" << mSolverName
            << "\" unfortunately the current compilation of Trilinos does not include it" << std::endl;
    }
 
    AmesosSolver(const AmesosSolver& Other) = delete;
 
    ~AmesosSolver() override = default;
 
 
    AmesosSolver& operator=(const AmesosSolver& Other) = delete;
 
 
    bool Solve(SparseMatrixType& rA, VectorType& rX, VectorType& rB) override
    {
        KRATOS_TRY
        rA.Comm().Barrier();
        Epetra_LinearProblem linear_problem(&rA,&rX,&rB);
        Amesos_BaseSolver* p_amesos_solver;
        Amesos amesos_factory;
        p_amesos_solver = amesos_factory.Create(mSolverName, linear_problem); // that the solver exists is checked in the constructor
 
        p_amesos_solver->SetParameters( mParameterList );
 
        p_amesos_solver->SymbolicFactorization();
        p_amesos_solver->NumericFactorization();
        p_amesos_solver->Solve();
 
        delete p_amesos_solver;
 
        rA.Comm().Barrier();
 
        return true;
        KRATOS_CATCH("");
    }
 
    bool Solve(SparseMatrixType& rA, DenseMatrixType& rX, DenseMatrixType& rB) override
    {
        return false;
    }
 
 
    static bool HasSolver(const std::string& rAmesosSolverName)
    {
        Amesos amesos_factory;
        return amesos_factory.Query(rAmesosSolverName);
    }
 
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "Trilinos Amesos-Solver";
    }
 
 
private:
 
    Teuchos::ParameterList mParameterList;
    std::string mSolverName;
 
 
}; // Class AmesosSolver
 
template<class TSparseSpaceType, class TDenseSpaceType, class TReordererType>
inline std::ostream& operator << (std::ostream& rOStream,
                                  const AmesosSolver<TSparseSpaceType,
                                  TDenseSpaceType, TReordererType>& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
}  // namespace Kratos.
 
#endif // KRATOS_AMESOS_SOLVER_H_INCLUDED defined