amgcl_solver.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                     Kratos default license: kratos/license.txt
//
//  Main authors:    Denis Demidov
//                   Riccardo Rossi
//
 
#if !defined(KRATOS_AMGCL_SOLVER )
#define  KRATOS_AMGCL_SOLVER
 
#ifndef AMGCL_PARAM_UNKNOWN
#include "input_output/logger.h"
#  define AMGCL_PARAM_UNKNOWN(name)                                            \
    Kratos::Logger("AMGCL") << KRATOS_CODE_LOCATION << Kratos::Logger::Severity::WARNING << "Unknown parameter " << name << std::endl
#endif
 
// System includes
#include <iostream>
#include <fstream>
#include <utility>
 
// External includes
/* BOOST */
#include <boost/range/iterator_range.hpp>
#include <boost/property_tree/json_parser.hpp>
 
// Project includes
#include "includes/define.h"
#include "includes/kratos_parameters.h"
#include "linear_solvers/iterative_solver.h"
#include "includes/ublas_interface.h"
#include "spaces/ublas_space.h"
 
#include <amgcl/coarsening/rigid_body_modes.hpp>
 
namespace Kratos
{
 
 
enum AMGCLSmoother
{
    SPAI0,SPAI1,ILU0,DAMPED_JACOBI,GAUSS_SEIDEL,CHEBYSHEV
};
 
enum AMGCLIterativeSolverType
{
    LGMRES,FGMRES,GMRES,BICGSTAB,CG,BICGSTAB_WITH_GMRES_FALLBACK,BICGSTAB2
};
 
enum AMGCLCoarseningType
{
    RUGE_STUBEN,AGGREGATION,SA,SA_EMIN
};
 
 
void KRATOS_API(KRATOS_CORE) AMGCLSolve(
    int block_size,
    TUblasSparseSpace<double>::MatrixType& rA,
    TUblasSparseSpace<double>::VectorType& rX,
    TUblasSparseSpace<double>::VectorType& rB,
    TUblasSparseSpace<double>::IndexType& rIterationNumber,
    double& rResidual,
    boost::property_tree::ptree amgclParams,
    int verbosity_level,
    bool use_gpgpu
    );
 
 
template< class TSparseSpaceType, class TDenseSpaceType,
          class TReordererType = Reorderer<TSparseSpaceType, TDenseSpaceType> >
class AMGCLSolver : public LinearSolver< TSparseSpaceType,
    TDenseSpaceType, TReordererType>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION( AMGCLSolver );
 
    typedef LinearSolver<TSparseSpaceType, TDenseSpaceType, TReordererType> BaseType;
 
    typedef typename TSparseSpaceType::MatrixType SparseMatrixType;
 
    typedef typename TSparseSpaceType::VectorType VectorType;
 
    typedef typename TDenseSpaceType::MatrixType DenseMatrixType;
 
    typedef ModelPart::DofsArrayType DofsArrayType;
 
    typedef typename TSparseSpaceType::IndexType IndexType;
 
    typedef std::size_t SizeType;
 
 
    AMGCLSolver(Parameters ThisParameters = Parameters(R"({})"))
    {
        Parameters default_parameters( R"(
        {
            "preconditioner_type"            : "amg",
            "solver_type"                    : "AMGCL",
            "smoother_type"                  : "ilu0",
            "krylov_type"                    : "gmres",
            "coarsening_type"                : "aggregation",
            "max_iteration"                  : 100,
            "provide_coordinates"            : false,
            "gmres_krylov_space_dimension"   : 100,
            "verbosity"                      : 1,
            "tolerance"                      : 1e-6,
            "scaling"                        : false,
            "block_size"                     : 1,
            "use_block_matrices_if_possible" : true,
            "coarse_enough"                  : 1000,
            "max_levels"                     : -1,
            "pre_sweeps"                     : 1,
            "post_sweeps"                    : 1,
            "use_gpgpu"                      : false
        }  )" );
 
        // Now validate agains defaults -- this also ensures no type mismatch
        ThisParameters.ValidateAndAssignDefaults(default_parameters);
 
        // specify available options
        std::set<std::string> available_smoothers = {"spai0","spai1","ilu0","ilut","iluk","damped_jacobi","gauss_seidel","chebyshev"};
        std::set<std::string> available_solvers = {"gmres","bicgstab","cg","bicgstabl","lgmres","fgmres", "bicgstab_with_gmres_fallback","idrs"};
        std::set<std::string> available_coarsening = {"ruge_stuben","aggregation","smoothed_aggregation","smoothed_aggr_emin"};
        std::set<std::string> available_preconditioner = {"amg","relaxation","dummy"};
 
        // Validate if values are admissible
        CheckIfSelectedOptionIsAvailable(ThisParameters, "smoother_type",       available_smoothers);
        CheckIfSelectedOptionIsAvailable(ThisParameters, "krylov_type",         available_solvers);
        CheckIfSelectedOptionIsAvailable(ThisParameters, "coarsening_type",     available_coarsening);
        CheckIfSelectedOptionIsAvailable(ThisParameters, "preconditioner_type", available_preconditioner);
 
        //selecting preconditioner type - default is AMG
        mAMGCLParameters.put("precond.class", ThisParameters["preconditioner_type"].GetString());
        if(ThisParameters["preconditioner_type"].GetString() != "amg"){
            mUseAMGPreconditioning = false;
        }
 
        if(ThisParameters["preconditioner_type"].GetString() == "relaxation") //this implies not using. Use a relaxation sweep as preconditioning. Relaxation type is taken from smoother_type
        {
            mAMGCLParameters.put("precond.type", ThisParameters["smoother_type"].GetString());
        }
 
        mProvideCoordinates = ThisParameters["provide_coordinates"].GetBool();
        mCoarseEnough = ThisParameters["coarse_enough"].GetInt();
 
        mBlockSize = ThisParameters["block_size"].GetInt(); //set the mndof to an inital number
        mTolerance = ThisParameters["tolerance"].GetDouble();
        mMaxIterationsNumber = ThisParameters["max_iteration"].GetInt();
        mVerbosity=ThisParameters["verbosity"].GetInt();
        mGMRESSize = ThisParameters["gmres_krylov_space_dimension"].GetInt();
 
        const std::string& solver_type = ThisParameters["krylov_type"].GetString();
        mAMGCLParameters.put("solver.type", solver_type);
        mFallbackToGMRES = false;
 
        if(solver_type == "bicgstab_with_gmres_fallback") {
            mFallbackToGMRES = true;
            mAMGCLParameters.put("solver.type", "bicgstab");
        }
 
        //settings only needed if full AMG is used
        if(mUseAMGPreconditioning)
        {
            mAMGCLParameters.put("precond.relax.type", ThisParameters["smoother_type"].GetString());
            mAMGCLParameters.put("precond.coarsening.type",  ThisParameters["coarsening_type"].GetString());
 
            int max_levels = ThisParameters["max_levels"].GetInt();
            if(max_levels >= 0)
                mAMGCLParameters.put("precond.max_levels",  max_levels);
 
            mAMGCLParameters.put("precond.npre",  ThisParameters["pre_sweeps"].GetInt());
            mAMGCLParameters.put("precond.npost",  ThisParameters["post_sweeps"].GetInt());
        }
 
        mUseBlockMatricesIfPossible = ThisParameters["use_block_matrices_if_possible"].GetBool();
 
        mUseGPGPU = ThisParameters["use_gpgpu"].GetBool();
    }
 
    AMGCLSolver(
        AMGCLSmoother Smoother,
        AMGCLIterativeSolverType Solver,
        double Tolerance,
        int MaxIterationsNumber,
        int Verbosity,
        int GMRESSize = 50
        ) : mTolerance(Tolerance),
            mMaxIterationsNumber(MaxIterationsNumber),
            mVerbosity(Verbosity),
            mBlockSize(1),
            mGMRESSize(GMRESSize),
            mCoarseEnough(1000),
            mFallbackToGMRES(false),
            mProvideCoordinates(false),
            mUseBlockMatricesIfPossible(false)
    {
        KRATOS_INFO_IF("AMGCL Linear Solver", mVerbosity > 0) << "Setting up AMGCL for iterative solve " << std::endl;
 
        // Choose smoother in the list "gauss_seidel, multicolor_gauss_seidel, ilu0, parallel_ilu0, ilut, damped_jacobi, spai0, chebyshev"
        SetSmootherType(Smoother);
        // Setting iterative solver
        SetIterativeSolverType(Solver);
 
        mAMGCLParameters.put("precond.coarsening.type", "aggregation");
    }
 
    AMGCLSolver(
        AMGCLSmoother Smoother,
        AMGCLIterativeSolverType Solver,
        AMGCLCoarseningType Coarsening,
        double Tolerance,
        int MaxIterationsNumber,
        int Verbosity,
        int GMRESSize = 50,
        bool ProvideCoordinates = false
        ) : mTolerance(Tolerance),
            mMaxIterationsNumber(MaxIterationsNumber),
            mVerbosity(Verbosity),
            mBlockSize(1),
            mGMRESSize(GMRESSize),
            mCoarseEnough(1000),
            mFallbackToGMRES(false),
            mProvideCoordinates(ProvideCoordinates),
            mUseBlockMatricesIfPossible(false)
    {
        KRATOS_INFO_IF("AMGCL Linear Solver", mVerbosity > 0) << "Setting up AMGCL for iterative solve " << std::endl;
 
        // Choose smoother in the list "gauss_seidel, multicolor_gauss_seidel, ilu0, parallel_ilu0, ilut, damped_jacobi, spai0, chebyshev"
        SetSmootherType(Smoother);
        // Setting iterative solver
        SetIterativeSolverType(Solver);
        // Setting coarsening type
        SetCoarseningType(Coarsening);
    }
 
    ~AMGCLSolver() override {};
 
 
 
    bool Solve(
        SparseMatrixType& rA,
        VectorType& rX,
        VectorType& rB
        ) override
    {
        // Initial checks
        KRATOS_ERROR_IF(TSparseSpaceType::Size1(rA) != TSparseSpaceType::Size2(rA) ) << "matrix A is not square! sizes are "
            << TSparseSpaceType::Size1(rA) << " and " << TSparseSpaceType::Size2(rA) << std::endl;
        KRATOS_ERROR_IF(TSparseSpaceType::Size(rX)  != TSparseSpaceType::Size1(rA)) << "size of x does not match the size of A. x size is " << TSparseSpaceType::Size(rX)
            << " matrix size is " << TSparseSpaceType::Size1(rA) << std::endl;
        KRATOS_ERROR_IF(TSparseSpaceType::Size(rB) != TSparseSpaceType::Size1(rA)) << "size of b does not match the size of A. b size is " << TSparseSpaceType::Size(rB)
            << " matrix size is " << TSparseSpaceType::Size1(rA) << std::endl;
 
        mAMGCLParameters.put("solver.tol", mTolerance);
        mAMGCLParameters.put("solver.maxiter", mMaxIterationsNumber);
 
        if(mUseAMGPreconditioning)
            mAMGCLParameters.put("precond.coarse_enough",mCoarseEnough/mBlockSize);
 
        // Use rigid body modes or set block size
        int static_block_size = mUseBlockMatricesIfPossible ? mBlockSize : 1;
        std::vector<double> B;
        if(mUseAMGPreconditioning && mProvideCoordinates && (mBlockSize == 2 || mBlockSize == 3)) {
            int nmodes = amgcl::coarsening::rigid_body_modes(mBlockSize,
                    boost::make_iterator_range(
                        &mCoordinates[0][0],
                        &mCoordinates[0][0] + TSparseSpaceType::Size1(rA)),
                    B);
 
            if (static_block_size != 1 && static_block_size != 3) {
                KRATOS_WARNING("AMGCL Linear Solver") << "Can only combine block matrices with coordinates in 3D. Falling back to scalar matrices" << std::endl;
                static_block_size = 1;
            }
 
            mAMGCLParameters.put("precond.coarsening.aggr.eps_strong", 0.0);
            mAMGCLParameters.put("precond.coarsening.aggr.block_size", 1);
            mAMGCLParameters.put("precond.coarsening.nullspace.cols",  nmodes);
            mAMGCLParameters.put("precond.coarsening.nullspace.rows",  TSparseSpaceType::Size1(rA));
            mAMGCLParameters.put("precond.coarsening.nullspace.B",     &B[0]);
        } else if(mUseAMGPreconditioning && mAMGCLParameters.get<std::string>("precond.coarsening.type") != std::string("ruge_stuben")) {
            mAMGCLParameters.put("precond.coarsening.aggr.eps_strong", 0.0);
            mAMGCLParameters.put("precond.coarsening.aggr.block_size", mBlockSize);
        }
 
        if (mVerbosity > 2) {
            write_json(std::cout, mAMGCLParameters);
        }
 
        if(mVerbosity == 4) {
            //output to matrix market
            std::stringstream matrix_market_name;
            matrix_market_name << "A" <<  ".mm";
            TSparseSpaceType::WriteMatrixMarketMatrix((char*) (matrix_market_name.str()).c_str(), rA, false);
 
            std::stringstream matrix_market_vectname;
            matrix_market_vectname << "b" << ".mm.rhs";
            TSparseSpaceType::WriteMatrixMarketVector((char*) (matrix_market_vectname.str()).c_str(), rB);
 
            if(mProvideCoordinates) {
                //output of coordinates
                std::ofstream coordsfile;
                coordsfile.open ("coordinates.txt");
                for(unsigned int i=0; i<mCoordinates.size(); i++) {
                    coordsfile << mCoordinates[i][0] << " " << mCoordinates[i][1] << " " << mCoordinates[i][2] << "\n";
                }
                coordsfile.close();
            }
 
            KRATOS_ERROR << " Verbosity = 4 prints the matrix and exits" << std::endl;
        }
 
        size_t iters;
        double resid;
        {
            if(mFallbackToGMRES) mAMGCLParameters.put("solver.type", "bicgstab"); //first we need to try with bicgstab
 
            if(mAMGCLParameters.get<std::string>("solver.type") == "gmres" ||
                mAMGCLParameters.get<std::string>("solver.type") == "lgmres" ||
                mAMGCLParameters.get<std::string>("solver.type") == "fgmres" )
                mAMGCLParameters.put("solver.M",  mGMRESSize);
            else
                mAMGCLParameters.erase("solver.M");
 
            if(mUseBlockMatricesIfPossible) {
                KRATOS_ERROR_IF(TSparseSpaceType::Size1(rA)%mBlockSize != 0) << "The block size employed " << mBlockSize << " is not an exact multiple of the matrix size "
                    << TSparseSpaceType::Size1(rA) << std::endl;
            }
            AMGCLSolve(static_block_size, rA,rX,rB, iters, resid, mAMGCLParameters, mVerbosity, mUseGPGPU);
        } //please do not remove this parenthesis!
 
        if(mFallbackToGMRES && resid > mTolerance ) {
            mAMGCLParameters.put("solver.type", "gmres");
            mAMGCLParameters.put("solver.M",  mGMRESSize);
            AMGCLSolve(1, rA,rX,rB, iters, resid, mAMGCLParameters, mVerbosity, mUseGPGPU);
        }
 
        KRATOS_WARNING_IF("AMGCL Linear Solver", mTolerance < resid)<<"Non converged linear solution. ["<< resid << " > "<< mTolerance << "]" << std::endl;
 
        KRATOS_INFO_IF("AMGCL Linear Solver", mVerbosity > 1)
            << "Iterations: " << iters << std::endl
            << "Error: "      << resid << std::endl;
 
        // Setting values
        SetResidualNorm(resid);
        SetIterationsNumber(iters);
 
        // We check the convergence
        if(resid > mTolerance)
            return false;
 
        return true;
    }
 
    virtual void SetIterationsNumber(const IndexType IterationsNumber)
    {
        mIterationsNumber = IterationsNumber;
    }
 
    virtual void SetResidualNorm(const double ResidualNorm)
    {
        mResidualNorm = ResidualNorm;
    }
 
    IndexType GetIterationsNumber() override
    {
        return mIterationsNumber;
    }
 
    virtual double GetResidualNorm()
    {
        return mResidualNorm;
    }
 
    bool Solve(SparseMatrixType& rA, DenseMatrixType& rX, DenseMatrixType& rB) override
    {
        return false;
    }
 
    bool AdditionalPhysicalDataIsNeeded() override
    {
        return true;
    }
 
    void ProvideAdditionalData (
        SparseMatrixType& rA,
        VectorType& rX,
        VectorType& rB,
        DofsArrayType& rDofSet,
        ModelPart& rModelPart
        ) override
    {
        int old_ndof = -1;
        int ndof=0;
 
        if (!rModelPart.IsDistributed())
        {
            unsigned int old_node_id = rDofSet.size() ? rDofSet.begin()->Id() : 0;
            for (auto it = rDofSet.begin(); it!=rDofSet.end(); it++) {
                if(it->EquationId() < TSparseSpaceType::Size1(rA) ) {
                    IndexType id = it->Id();
                    if(id != old_node_id) {
                        old_node_id = id;
                        if(old_ndof == -1) old_ndof = ndof;
                        else if(old_ndof != ndof) { //if it is different than the block size is 1
                            old_ndof = -1;
                            break;
                        }
 
                        ndof=1;
                    } else {
                        ndof++;
                    }
                }
            }
 
            if(old_ndof == -1)
                mBlockSize = 1;
            else
                mBlockSize = ndof;
 
        }
        else //distribute
        {
            const std::size_t system_size = TSparseSpaceType::Size1(rA);
            int current_rank = rModelPart.GetCommunicator().GetDataCommunicator().Rank();
            unsigned int old_node_id = rDofSet.size() ? rDofSet.begin()->Id() : 0;
            for (auto it = rDofSet.begin(); it!=rDofSet.end(); it++) {
                if(it->EquationId() < system_size  && it->GetSolutionStepValue(PARTITION_INDEX) == current_rank) {
                    IndexType id = it->Id();
                    if(id != old_node_id) {
                        old_node_id = id;
                        if(old_ndof == -1) old_ndof = ndof;
                        else if(old_ndof != ndof) { //if it is different than the block size is 1
                            old_ndof = -1;
                            break;
                        }
 
                        ndof=1;
                    } else {
                        ndof++;
                    }
                }
            }
 
            if(old_ndof != -1)
                mBlockSize = ndof;
 
            int max_block_size = rModelPart.GetCommunicator().GetDataCommunicator().MaxAll(mBlockSize);
 
            if( old_ndof == -1) {
                mBlockSize = max_block_size;
            }
 
            KRATOS_ERROR_IF(mBlockSize != max_block_size) << "Block size is not consistent. Local: " << mBlockSize  << " Max: " << max_block_size << std::endl;
        }
 
        KRATOS_INFO_IF("AMGCL Linear Solver", mVerbosity > 1) << "mndof: " << mBlockSize << std::endl;
 
        if(mProvideCoordinates) {
            mCoordinates.resize(TSparseSpaceType::Size1(rA)/mBlockSize);
            unsigned int i=0;
            for (auto it_dof = rDofSet.begin(); it_dof!=rDofSet.end(); it_dof+=mBlockSize) {
                if(it_dof->EquationId() < TSparseSpaceType::Size1(rA) ) {
                    auto it_node = rModelPart.Nodes().find(it_dof->Id());
                    mCoordinates[ i ] = it_node->Coordinates();
                    i++;
                }
            }
        }
    }
 
 
 
 
 
 
    void  PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << "AMGCL solver:";
    }
 
    void  PrintData(std::ostream& rOStream) const override
    {
        rOStream << "Settings: ";
        write_json(rOStream, mAMGCLParameters);
    }
 
 
 
 
protected:
 
 
 
 
    // Helper function for checking if a selected option is available
    // and printing the available options
    void CheckIfSelectedOptionIsAvailable(
        const Parameters ThisParameters,
        const std::string& rOptionName,
        const std::set<std::string>& rAvailableOptions)
    {
        if (rAvailableOptions.find(ThisParameters[rOptionName].GetString()) == rAvailableOptions.end()) {
            std::stringstream msg;
            msg << "Currently prescribed " << rOptionName << " : " << ThisParameters[rOptionName].GetString() << std::endl;
            msg << "Admissible values are :";
            for (const auto& r_name : rAvailableOptions) {
                msg << std::endl << "    " << r_name;
            }
            KRATOS_ERROR << "AMGCL Linear Solver : " << rOptionName << " is invalid!" << std::endl << msg.str() << std::endl;
        }
    }
 
 
 
 
 
 
 
    double mTolerance;                
    IndexType mMaxIterationsNumber;   
    int mVerbosity;                   
    int mBlockSize;                   
    SizeType mGMRESSize;              
    SizeType mCoarseEnough;           
    bool mFallbackToGMRES;            
    bool mProvideCoordinates;         
    bool mUseBlockMatricesIfPossible; 
    bool mUseGPGPU;                   
 
    std::vector<array_1d<double,3> > mCoordinates; 
 
    boost::property_tree::ptree mAMGCLParameters;  
 
    double mResidualNorm = 0.0;      
    IndexType mIterationsNumber = 0; 
    bool mUseAMGPreconditioning = true; 
 
 
    void SetSmootherType(const AMGCLSmoother SmootherType)
    {
        switch(SmootherType)
        {
            case SPAI0:
            {
                mAMGCLParameters.put("precond.relax.type","spai0");
                break;
            }
            case SPAI1:
            {
                mAMGCLParameters.put("precond.relax.type","spai1");
                break;
            }
            case ILU0:
            {
                mAMGCLParameters.put("precond.relax.type","ilu0");
                break;
            }
            case DAMPED_JACOBI:
            {
                mAMGCLParameters.put("precond.relax.type","damped_jacobi");
                break;
            }
            case GAUSS_SEIDEL:
            {
                mAMGCLParameters.put("precond.relax.type","gauss_seidel");
                break;
            }
            case CHEBYSHEV:
            {
                mAMGCLParameters.put("precond.relax.type","chebyshev");
                break;
            }
        };
    }
 
    void SetIterativeSolverType(const AMGCLIterativeSolverType SolverType)
    {
        switch(SolverType)
        {
            case GMRES:
            {
                mAMGCLParameters.put("solver.M",  mGMRESSize);
                mAMGCLParameters.put("solver.type", "gmres");
                break;
            }
            case FGMRES:
            {
                mAMGCLParameters.put("solver.M",  mGMRESSize);
                mAMGCLParameters.put("solver.type", "fgmres");
                break;
            }
            case LGMRES:
            {
                mAMGCLParameters.put("solver.M",  mGMRESSize);
                mAMGCLParameters.put("solver.type", "lgmres");
                break;
            }
            case BICGSTAB:
            {
                mAMGCLParameters.put("solver.type", "bicgstab");
                break;
            }
            case CG:
            {
                mAMGCLParameters.put("solver.type", "cg");
                break;
            }
            case BICGSTAB2:
            {
                mAMGCLParameters.put("solver.type", "bicgstabl");
                break;
            }
            case BICGSTAB_WITH_GMRES_FALLBACK:
            {
                mAMGCLParameters.put("solver.M",  mGMRESSize);
                mAMGCLParameters.put("solver.type", "bicgstab");
                mFallbackToGMRES=true;
                break;
            }
        };
    }
 
    void SetCoarseningType(const AMGCLCoarseningType CoarseningType)
    {
        switch(CoarseningType)
        {
            case RUGE_STUBEN:
            {
                mAMGCLParameters.put("precond.coarsening.type", "ruge_stuben");
                break;
            }
            case AGGREGATION:
            {
                mAMGCLParameters.put("precond.coarsening.type", "aggregation");
                break;
            }
            case SA:
            {
                mAMGCLParameters.put("precond.coarsening.type", "smoothed_aggregation");
                break;
            }
            case SA_EMIN:
            {
                mAMGCLParameters.put("precond.coarsening.type", "smoothed_aggr_emin");
                break;
            }
        };
    }
 
 
 
 
 
 
 
    AMGCLSolver& operator=(const AMGCLSolver& Other);
 
    AMGCLSolver(const AMGCLSolver& Other);
 
}; // Class AMGCLSolver
 
 
template<class TSparseSpaceType, class TDenseSpaceType,class TReordererType>
inline std::istream& operator >> (std::istream& rIStream, AMGCLSolver< TSparseSpaceType,
                                  TDenseSpaceType, TReordererType>& rThis)
{
    return rIStream;
}
 
template<class TSparseSpaceType, class TDenseSpaceType, class TReordererType>
inline std::ostream& operator << (std::ostream& rOStream,
                                  const AMGCLSolver<TSparseSpaceType,
                                  TDenseSpaceType, TReordererType>& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
}  // namespace Kratos.
 
 
#endif // KRATOS_AMGCL_SOLVER  defined