Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Register Single Precision Linear Solvers #13087

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
+250 −105
Open

Register Single Precision Linear Solvers #13087

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
40c3bfb
add support for reading/writing MatrixMarket in single precision
matekelemen Feb 4, 2025
3a44c9c
fix templated core solvers' value types
matekelemen Feb 4, 2025
59e48e9
register single precision versions of templated core linear solvers
matekelemen Feb 4, 2025
42aafd3
fix more sparse space value type inconsistencies
matekelemen Feb 4, 2025
dca7e3e
extend serializer for std::complex<float>
matekelemen Feb 4, 2025
83c3e3c
expose single precision linear solvers in KratosComponents
matekelemen Feb 15, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
124 changes: 74 additions & 50 deletions kratos/factories/standard_linear_solver_factory.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,7 @@
// External includes

// Project includes
#include "includes/define.h"
#include "factories/standard_linear_solver_factory.h"
#include "linear_solvers/linear_solver.h"
#include "linear_solvers/cg_solver.h"
#include "linear_solvers/deflated_cg_solver.h"
#include "linear_solvers/bicgstab_solver.h"
Expand All @@ -31,55 +29,81 @@
#include "linear_solvers/skyline_lu_custom_scalar_solver.h"
#include "spaces/ublas_space.h"

namespace Kratos
namespace Kratos {


namespace detail {


template <class TSparseDataType,
class TDenseDataType>
void RegisterLinearSolvers()
{
void RegisterLinearSolvers()
{
using SpaceType = TUblasSparseSpace<double>;
using LocalSpaceType = TUblasDenseSpace<double>;
using ComplexSpaceType = TUblasSparseSpace<std::complex<double>>;
using ComplexLocalSpaceType = TUblasDenseSpace<std::complex<double>>;
// using LinearSolverType = LinearSolver<SpaceType, LocalSpaceType>;
// using IterativeSolverType = IterativeSolver<SpaceType, LocalSpaceType>;
using CGSolverType = CGSolver<SpaceType, LocalSpaceType>;
using SpaceType = TUblasSparseSpace<TSparseDataType>;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

BTW, we should eventually replace the linear solver factory with the register @KratosMultiphysics/technical-committee I don't know if a good moment to do that

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was thinking about that too, but other parts of Kratos still expect linear solvers to be in KratosComponents and I didn't want to change too many unrelated things at once.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A priori we can define in both database and add a warning in the factory, opinion @KratosMultiphysics/technical-committee ?

using LocalSpaceType = TUblasDenseSpace<TDenseDataType>;
using ComplexSpaceType = TUblasSparseSpace<std::complex<TSparseDataType>>;
using ComplexLocalSpaceType = TUblasDenseSpace<std::complex<TDenseDataType>>;

using CGSolverType = CGSolver<SpaceType, LocalSpaceType>;
static auto CGSolverFactory = StandardLinearSolverFactory<SpaceType,LocalSpaceType,CGSolverType>();
KratosComponents<LinearSolverFactory<SpaceType,LocalSpaceType>>::Add("cg", CGSolverFactory);

using BICGSTABSolverType = BICGSTABSolver<SpaceType, LocalSpaceType>;
static auto BICGSTABSolverFactory = StandardLinearSolverFactory<SpaceType,LocalSpaceType,BICGSTABSolverType>();
KratosComponents<LinearSolverFactory<SpaceType,LocalSpaceType>>::Add("bicgstab", BICGSTABSolverFactory);

using SkylineLUFactorizationSolverType = SkylineLUFactorizationSolver<SpaceType, LocalSpaceType>;
static auto SkylineLUFactorizationSolverFactory= StandardLinearSolverFactory<SpaceType,LocalSpaceType,SkylineLUFactorizationSolverType>();
KratosComponents<LinearSolverFactory<SpaceType,LocalSpaceType>>::Add("skyline_lu_factorization",SkylineLUFactorizationSolverFactory );

using TFQMRSolverType = TFQMRSolver<SpaceType, LocalSpaceType>;
static auto TFQMRSolverFactory = StandardLinearSolverFactory<SpaceType,LocalSpaceType,TFQMRSolverType>();
KratosComponents<LinearSolverFactory<SpaceType,LocalSpaceType>>::Add("tfqmr", TFQMRSolverFactory);

using AMGCLSolverType = AMGCLSolver<SpaceType, LocalSpaceType>;
static auto AMGCLSolverFactory= StandardLinearSolverFactory<SpaceType,LocalSpaceType,AMGCLSolverType>();
KratosComponents<LinearSolverFactory<SpaceType,LocalSpaceType>>::Add("amgcl", AMGCLSolverFactory);

using AMGCL_NS_SolverType = AMGCL_NS_Solver<SpaceType, LocalSpaceType>;
static auto AMGCL_NS_SolverFactory= StandardLinearSolverFactory<SpaceType,LocalSpaceType,AMGCL_NS_SolverType>();
KratosComponents<LinearSolverFactory<SpaceType,LocalSpaceType>>::Add("amgcl_ns",AMGCL_NS_SolverFactory );

using ScalingSolverType = ScalingSolver<SpaceType, LocalSpaceType>;
static auto ScalingSolverFactory= StandardLinearSolverFactory<SpaceType,LocalSpaceType,ScalingSolverType>();
KratosComponents<LinearSolverFactory<SpaceType,LocalSpaceType>>::Add("scaling",ScalingSolverFactory );

using FallbackLinearSolverType = FallbackLinearSolver<SpaceType, LocalSpaceType>;
static auto FallbackLinearSolverFactory = StandardLinearSolverFactory<SpaceType, LocalSpaceType, FallbackLinearSolverType>();
KratosComponents<LinearSolverFactory<SpaceType,LocalSpaceType>>::Add("fallback_linear_solver", FallbackLinearSolverFactory);

using MonotonicityPreservingSolverType = MonotonicityPreservingSolver<SpaceType, LocalSpaceType>;
static auto MonotonicityPreservingSolverFactory= StandardLinearSolverFactory<SpaceType,LocalSpaceType,MonotonicityPreservingSolverType>();
KratosComponents<LinearSolverFactory<SpaceType,LocalSpaceType>>::Add("monotonicity_preserving",MonotonicityPreservingSolverFactory );

using SkylineLUComplexSolverType = SkylineLUCustomScalarSolver<ComplexSpaceType, ComplexLocalSpaceType>;
static auto SkylineLUComplexSolverFactory = StandardLinearSolverFactory<ComplexSpaceType, ComplexLocalSpaceType, SkylineLUComplexSolverType>();
KratosComponents<LinearSolverFactory<ComplexSpaceType,ComplexLocalSpaceType>>::Add("skyline_lu_complex", SkylineLUComplexSolverFactory);

// using LinearSolverType = LinearSolver<SpaceType, LocalSpaceType>;
// using IterativeSolverType = IterativeSolver<SpaceType, LocalSpaceType>;
// KRATOS_REGISTER_LINEAR_SOLVER("LinearSolver", StandardLinearSolverFactory<SpaceType,LocalSpaceType,LinearSolverType>());

// DeflatedCGSolver's utilities are not templated on the space type (and thus the floating point type),
// so they can only be defined for double precision floats.
if constexpr (std::is_same_v<TSparseDataType,double>) {
using DeflatedCGSolverType = DeflatedCGSolver<SpaceType, LocalSpaceType>;
using BICGSTABSolverType = BICGSTABSolver<SpaceType, LocalSpaceType>;
using TFQMRSolverType = TFQMRSolver<SpaceType, LocalSpaceType>;
using SkylineLUFactorizationSolverType = SkylineLUFactorizationSolver<SpaceType, LocalSpaceType>;
using AMGCLSolverType = AMGCLSolver<SpaceType, LocalSpaceType>;
using AMGCL_NS_SolverType = AMGCL_NS_Solver<SpaceType, LocalSpaceType>;
using SkylineLUComplexSolverType = SkylineLUCustomScalarSolver<ComplexSpaceType, ComplexLocalSpaceType>;

using ScalingSolverType = ScalingSolver<SpaceType, LocalSpaceType>;
using FallbackLinearSolverType = FallbackLinearSolver<SpaceType, LocalSpaceType>;
using MonotonicityPreservingSolverType = MonotonicityPreservingSolver<SpaceType, LocalSpaceType>;

//NOTE: here we must create persisting objects for the linear solvers
static auto CGSolverFactory = StandardLinearSolverFactory<SpaceType,LocalSpaceType,CGSolverType>();
static auto BICGSTABSolverFactory = StandardLinearSolverFactory<SpaceType,LocalSpaceType,BICGSTABSolverType>();
static auto DeflatedCGSolverFactory= StandardLinearSolverFactory<SpaceType,LocalSpaceType,DeflatedCGSolverType>();
static auto SkylineLUFactorizationSolverFactory= StandardLinearSolverFactory<SpaceType,LocalSpaceType,SkylineLUFactorizationSolverType>();
static auto TFQMRSolverFactory = StandardLinearSolverFactory<SpaceType,LocalSpaceType,TFQMRSolverType>();
static auto AMGCLSolverFactory= StandardLinearSolverFactory<SpaceType,LocalSpaceType,AMGCLSolverType>();
static auto AMGCL_NS_SolverFactory= StandardLinearSolverFactory<SpaceType,LocalSpaceType,AMGCL_NS_SolverType>();
static auto ScalingSolverFactory= StandardLinearSolverFactory<SpaceType,LocalSpaceType,ScalingSolverType>();
static auto FallbackLinearSolverFactory = StandardLinearSolverFactory<SpaceType, LocalSpaceType, FallbackLinearSolverType>();
static auto MonotonicityPreservingSolverFactory= StandardLinearSolverFactory<SpaceType,LocalSpaceType,MonotonicityPreservingSolverType>();
static auto SkylineLUComplexSolverFactory = StandardLinearSolverFactory<ComplexSpaceType, ComplexLocalSpaceType, SkylineLUComplexSolverType>();

//registration of linear solvers
// KRATOS_REGISTER_LINEAR_SOLVER("LinearSolver", StandardLinearSolverFactory<SpaceType,LocalSpaceType,LinearSolverType>());
KRATOS_REGISTER_LINEAR_SOLVER("cg", CGSolverFactory);
KRATOS_REGISTER_LINEAR_SOLVER("bicgstab", BICGSTABSolverFactory);
KRATOS_REGISTER_LINEAR_SOLVER("deflated_cg", DeflatedCGSolverFactory);
KRATOS_REGISTER_LINEAR_SOLVER("tfqmr", TFQMRSolverFactory);
KRATOS_REGISTER_LINEAR_SOLVER("skyline_lu_factorization",SkylineLUFactorizationSolverFactory );
KRATOS_REGISTER_LINEAR_SOLVER("amgcl", AMGCLSolverFactory);
KRATOS_REGISTER_LINEAR_SOLVER("amgcl_ns",AMGCL_NS_SolverFactory );
KRATOS_REGISTER_LINEAR_SOLVER("scaling",ScalingSolverFactory );
KRATOS_REGISTER_LINEAR_SOLVER("fallback_linear_solver", FallbackLinearSolverFactory);
KRATOS_REGISTER_LINEAR_SOLVER("monotonicity_preserving",MonotonicityPreservingSolverFactory );
KRATOS_REGISTER_COMPLEX_LINEAR_SOLVER("skyline_lu_complex", SkylineLUComplexSolverFactory);

};
KratosComponents<LinearSolverFactory<SpaceType,LocalSpaceType>>::Add("deflated_cg", DeflatedCGSolverFactory);
}
}


} // namespace detail


void RegisterLinearSolvers()
{
detail::RegisterLinearSolvers</*TSparseDataType=*/double,/*TDenseDataType*/double>();
detail::RegisterLinearSolvers</*TSparseDataType=*/float,/*TDenseDataType*/double>();
};
} // Namespace Kratos
2 changes: 2 additions & 0 deletions kratos/includes/serializer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -541,6 +541,7 @@ class KRATOS_API(KRATOS_CORE) Serializer
KRATOS_SERIALIZATION_DIRECT_LOAD(std::size_t)
#endif
KRATOS_SERIALIZATION_DIRECT_LOAD(std::complex<double>)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should remove those evil tabs

KRATOS_SERIALIZATION_DIRECT_LOAD(std::complex<float>)

template<class TDataType, std::size_t TDataSize>
void save(std::string const & rTag, std::array<TDataType, TDataSize> const& rObject)
Expand Down Expand Up @@ -787,6 +788,7 @@ class KRATOS_API(KRATOS_CORE) Serializer
KRATOS_SERIALIZATION_DIRECT_SAVE(std::size_t)
#endif
KRATOS_SERIALIZATION_DIRECT_SAVE(std::complex<double>)
KRATOS_SERIALIZATION_DIRECT_SAVE(std::complex<float>)


template<class TDataType>
Expand Down
22 changes: 20 additions & 2 deletions kratos/linear_solvers/amgcl_solver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,21 @@ void KRATOS_API(KRATOS_CORE) AMGCLSolve(
bool use_gpgpu
);


/// @copdoc AMGCLSolve(int,TUblasSparseSpace<float>::MatrixType&,TUblasSparseSpace<float>::VectorType&,TUblasSparseSpace<float>::VectorType&,TUblasSparseSpace<float>::IndexType&,float&,boost::property_tree::ptree,int,bool)
void KRATOS_API(KRATOS_CORE) AMGCLSolve(
int block_size,
TUblasSparseSpace<float>::MatrixType& rA,
TUblasSparseSpace<float>::VectorType& rX,
TUblasSparseSpace<float>::VectorType& rB,
TUblasSparseSpace<float>::IndexType& rIterationNumber,
float& rResidual,
boost::property_tree::ptree amgclParams,
int verbosity_level,
bool use_gpgpu
);


///@}
///@name Kratos Classes
///@{
Expand Down Expand Up @@ -401,7 +416,7 @@ class AMGCLSolver : public LinearSolver< TSparseSpaceType,
}

size_t iters;
double resid;
typename TSparseSpaceType::DataType resid;
{
if(mFallbackToGMRES) mAMGCLParameters.put("solver.type", "bicgstab"); //first we need to try with bicgstab

Expand All @@ -416,7 +431,10 @@ class AMGCLSolver : public LinearSolver< TSparseSpaceType,
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);
AMGCLSolve(static_block_size,
rA, rX, rB,
iters, resid,
mAMGCLParameters, mVerbosity, mUseGPGPU);
} //please do not remove this parenthesis!

if(mFallbackToGMRES && resid > mTolerance ) {
Expand Down
112 changes: 80 additions & 32 deletions kratos/linear_solvers/amgcl_solver_impl.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,19 +33,20 @@ vex::Context& vexcl_context() {
return ctx;
}

template <int TBlockSize>
template <class TValue, int TBlockSize>
void register_vexcl_static_matrix_type() {
static vex::scoped_program_header header(vexcl_context(),
amgcl::backend::vexcl_static_matrix_declaration<double,TBlockSize>());
amgcl::backend::vexcl_static_matrix_declaration<TValue,TBlockSize>());
}
#endif

template <class TValue>
void AMGCLScalarSolve(
TUblasSparseSpace<double>::MatrixType& rA,
TUblasSparseSpace<double>::VectorType& rX,
TUblasSparseSpace<double>::VectorType& rB,
TUblasSparseSpace<double>::IndexType& rIterationNumber,
double& rResidual,
typename TUblasSparseSpace<TValue>::MatrixType& rA,
typename TUblasSparseSpace<TValue>::VectorType& rX,
typename TUblasSparseSpace<TValue>::VectorType& rB,
typename TUblasSparseSpace<TValue>::IndexType& rIterationNumber,
TValue& rResidual,
const boost::property_tree::ptree &amgclParams,
int verbosity_level,
bool use_gpgpu
Expand All @@ -55,24 +56,24 @@ void AMGCLScalarSolve(
if (use_gpgpu && vexcl_context()) {
auto &ctx = vexcl_context();

typedef amgcl::backend::vexcl<double> Backend;
typedef amgcl::backend::vexcl<TValue> Backend;
typedef amgcl::make_solver<
amgcl::runtime::preconditioner<Backend>,
amgcl::runtime::solver::wrapper<Backend>
> Solver;

Backend::params bprm;
typename Backend::params bprm;
bprm.q = ctx;

Solver solve(amgcl::adapter::zero_copy(
TUblasSparseSpace<double>::Size1(rA),
TUblasSparseSpace<TValue>::Size1(rA),
rA.index1_data().begin(),
rA.index2_data().begin(),
rA.value_data().begin()),
amgclParams, bprm);

vex::vector<double> b(ctx, rB.size(), &rB[0]);
vex::vector<double> x(ctx, rX.size(), &rX[0]);
vex::vector<TValue> b(ctx, rB.size(), &rB[0]);
vex::vector<TValue> x(ctx, rX.size(), &rX[0]);

std::tie(rIterationNumber, rResidual) = solve(b, x);

Expand All @@ -83,7 +84,7 @@ void AMGCLScalarSolve(
} else
#endif
{
typedef amgcl::backend::builtin<double> Backend;
typedef amgcl::backend::builtin<TValue> Backend;
typedef amgcl::make_solver<
amgcl::runtime::preconditioner<Backend>,
amgcl::runtime::solver::wrapper<Backend>
Expand All @@ -103,13 +104,13 @@ void AMGCLScalarSolve(
}
}

template <int TBlockSize>
template <class TValue, int TBlockSize>
void AMGCLBlockSolve(
TUblasSparseSpace<double>::MatrixType & rA,
TUblasSparseSpace<double>::VectorType& rX,
TUblasSparseSpace<double>::VectorType& rB,
TUblasSparseSpace<double>::IndexType& rIterationNumber,
double& rResidual,
typename TUblasSparseSpace<TValue>::MatrixType & rA,
typename TUblasSparseSpace<TValue>::VectorType& rX,
typename TUblasSparseSpace<TValue>::VectorType& rB,
typename TUblasSparseSpace<TValue>::IndexType& rIterationNumber,
TValue& rResidual,
boost::property_tree::ptree amgclParams,
int verbosity_level,
bool use_gpgpu
Expand All @@ -120,16 +121,16 @@ void AMGCLBlockSolve(
else
amgclParams.put("precond.coarsening.aggr.block_size",1);

typedef amgcl::static_matrix<double, TBlockSize, TBlockSize> value_type;
typedef amgcl::static_matrix<double, TBlockSize, 1> rhs_type;
typedef amgcl::static_matrix<TValue, TBlockSize, TBlockSize> value_type;
typedef amgcl::static_matrix<TValue, TBlockSize, 1> rhs_type;

std::size_t n = TUblasSparseSpace<double>::Size1(rA);
std::size_t n = TUblasSparseSpace<TValue>::Size1(rA);
std::size_t nb = n / TBlockSize;

#ifdef AMGCL_GPGPU
if (use_gpgpu && vexcl_context()) {
auto &ctx = vexcl_context();
register_vexcl_static_matrix_type<TBlockSize>();
register_vexcl_static_matrix_type<TValue,TBlockSize>();

typedef amgcl::backend::vexcl<value_type> Backend;

Expand Down Expand Up @@ -186,13 +187,18 @@ void AMGCLBlockSolve(
}
}


namespace detail {


template <class TValue>
void AMGCLSolve(
int block_size,
TUblasSparseSpace<double>::MatrixType& rA,
TUblasSparseSpace<double>::VectorType& rX,
TUblasSparseSpace<double>::VectorType& rB,
TUblasSparseSpace<double>::IndexType& rIterationNumber,
double& rResidual,
typename TUblasSparseSpace<TValue>::MatrixType& rA,
typename TUblasSparseSpace<TValue>::VectorType& rX,
typename TUblasSparseSpace<TValue>::VectorType& rB,
typename TUblasSparseSpace<TValue>::IndexType& rIterationNumber,
TValue& rResidual,
boost::property_tree::ptree amgclParams,
int verbosity_level,
bool use_gpgpu
Expand All @@ -210,18 +216,60 @@ void AMGCLSolve(

switch (block_size) {
case 2:
AMGCLBlockSolve<2>(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
AMGCLBlockSolve<TValue,2>(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
return;
case 3:
AMGCLBlockSolve<3>(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
AMGCLBlockSolve<TValue,3>(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
return;
case 4:
AMGCLBlockSolve<4>(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
AMGCLBlockSolve<TValue,4>(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
return;
default:
AMGCLScalarSolve(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
AMGCLScalarSolve<TValue>(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
return;
}
}


} // namespace detail


void 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
)
{
detail::AMGCLSolve<double>(block_size,
rA, rX, rB,
rIterationNumber, rResidual,
amgclParams, verbosity_level, use_gpgpu);
}


void AMGCLSolve(
int block_size,
TUblasSparseSpace<float>::MatrixType& rA,
TUblasSparseSpace<float>::VectorType& rX,
TUblasSparseSpace<float>::VectorType& rB,
TUblasSparseSpace<float>::IndexType& rIterationNumber,
float& rResidual,
boost::property_tree::ptree amgclParams,
int verbosity_level,
bool use_gpgpu
)
{
detail::AMGCLSolve<float>(block_size,
rA, rX, rB,
rIterationNumber, rResidual,
amgclParams, verbosity_level, use_gpgpu);
}


}
Loading
Loading