23#include <ginkgo/ginkgo.hpp>
32#define MFEM_GINKGO_VERSION \
33 ((GKO_VERSION_MAJOR*100 + GKO_VERSION_MINOR)*100 + GKO_VERSION_PATCH)
40template <
typename T>
using gko_array = gko::array<T>;
70 gko::size_type size,
Vector *mfem_vec,
71 bool ownership =
false)
72 : gko::matrix::Dense<
real_t>(
74 gko::
dim<2> {size, 1},
78 exec != exec->get_master() ?
true : false)),
89 wrapped_vec = std::unique_ptr<
Vector,
90 std::function<void(
Vector *)>>(
95 using deleter = gko::null_deleter<Vector>;
96 wrapped_vec = std::unique_ptr<Vector,
97 std::function<void(Vector *)>>(
102 static std::unique_ptr<VectorWrapper>
create(
103 std::shared_ptr<const gko::Executor> exec,
106 bool ownership =
false)
108 return std::unique_ptr<VectorWrapper>(
120 std::unique_ptr<gko::matrix::Dense<real_t>>
125 this->wrapped_vec.get()->GetMemory().GetMemoryType());
127 mfem_vec->
UseDevice(this->wrapped_vec.get()->UseDevice());
143 std::shared_ptr<const gko::Executor> exec,
144 const gko::dim<2> &size,
145 gko::size_type stride)
const override
152 "VectorWrapper cannot be created with stride > 1");
155 gko::size_type total_size = size[0]*size[1];
158 this->wrapped_vec.get()->GetMemory().GetMemoryType());
160 mfem_vec->
UseDevice(this->wrapped_vec.get()->UseDevice());
166 this->get_executor(), total_size, mfem_vec,
173 const gko::span &rows,
174 const gko::span &columns,
175 const gko::size_type stride)
override
178 gko::size_type num_rows = rows.end - rows.begin;
179 gko::size_type num_cols = columns.end - columns.begin;
183 if (num_cols > 1 || stride > 1)
185 throw gko::BadDimension(
186 __FILE__, __LINE__, __func__,
"new_submatrix", num_rows,
188 "VectorWrapper submatrix must have one column and stride = 1");
190 int data_size =
static_cast<int>(num_rows * num_cols);
191 int start =
static_cast<int>(rows.begin);
194 mfem_vec->
MakeRef(*(this->wrapped_vec.get()), start, data_size);
195 mfem_vec->
UseDevice(this->wrapped_vec.get()->UseDevice());
202 this->get_executor(), data_size, mfem_vec,
207 std::unique_ptr<
Vector, std::function<void(
Vector *)>> wrapped_vec;
219 :
public gko::EnableLinOp<OperatorWrapper>,
220 public gko::EnableCreateMethod<OperatorWrapper>
224 gko::size_type size = 0,
227 gko::EnableCreateMethod<OperatorWrapper>()
229 this->wrapped_oper = oper;
233 void apply_impl(
const gko::LinOp *
b, gko::LinOp *x)
const override;
235 const gko::LinOp *
beta, gko::LinOp *x)
const override;
243template <
typename ValueType=real_t>
247 auto cpu_norm = clone(norm->get_executor()->get_master(), norm);
249 return cpu_norm->at(0, 0);
254template <
typename ValueType=real_t>
258 auto exec =
b->get_executor();
260 auto b_norm = gko::initialize<gko::matrix::Dense<ValueType>>({0.0}, exec);
262 b->compute_norm2(b_norm);
264 return std::pow(
get_norm(b_norm.get()),2);
278template <
typename ValueType=real_t>
285 if (compute_real_residual)
287 mfem::out <<
"Iteration log with real residual norms:" << std::endl;
291 mfem::out <<
"Iteration log with residual norms:" << std::endl;
293 mfem::out <<
'|' << std::setw(10) <<
"Iteration" <<
'|' << std::setw(25)
294 <<
"Residual Norm" <<
'|' << std::endl;
297 mfem::out <<
'|' << std::setfill(
'-') << std::setw(11) <<
'|' <<
298 std::setw(26) <<
'|' << std::setfill(
' ') << std::endl;
301 for (std::size_t i = 0; i < iterations.size(); i++)
303 mfem::out <<
'|' << std::setw(10) << iterations[i] <<
'|'
304 << std::setw(25) << residual_norms[i] <<
'|' << std::endl;
308 mfem::out.unsetf(std::ios_base::floatfield);
310 mfem::out <<
'|' << std::setfill(
'-') << std::setw(11) <<
'|' <<
311 std::setw(26) <<
'|' << std::setfill(
' ') << std::endl;
318 const gko::size_type &iteration,
319 const gko::LinOp *residual,
320 const gko::LinOp *solution,
321 const gko::LinOp *residual_norm)
const override
323 iteration_complete_core(iteration, residual, solution, residual_norm,
328 const gko::size_type &iteration,
329 const gko::LinOp *residual,
330 const gko::LinOp *solution,
331 const gko::LinOp *residual_norm,
332 const gko::LinOp *implicit_sq_residual_norm)
const override
334 iteration_complete_core(iteration, residual, solution, residual_norm,
335 implicit_sq_residual_norm);
339 const gko::LinOp *rhs,
340 const gko::LinOp *solution,
341 const gko::size_type &iteration,
342 const gko::LinOp *residual,
343 const gko::LinOp *residual_norm,
344 const gko::LinOp *implicit_sq_residual_norm,
345 const gko::array<gko::stopping_status>* status,
346 bool stopped)
const override
348 iteration_complete_core(iteration, residual, solution, residual_norm,
349 implicit_sq_residual_norm);
354 const gko::LinOp *matrix,
const gko_dense *b,
355 bool compute_real_residual=
false)
357 gko::log::Logger(gko::log::Logger::iteration_complete_mask),
360 compute_real_residual{compute_real_residual}
362 if (compute_real_residual ==
true)
367 res = std::move(gko_dense::create_with_config_of(b_cast).release());
371 res = std::move(gko::clone(
b).release());
379 void iteration_complete_core(
const gko::size_type &iteration,
380 const gko::LinOp *residual,
381 const gko::LinOp *solution,
382 const gko::LinOp *residual_norm,
383 const gko::LinOp *implicit_sq_residual_norm)
const
387 if (solution && compute_real_residual)
390 auto exec = matrix->get_executor();
393 matrix->apply(solution, res);
398 const VectorWrapper *b_cast = gko::as<const VectorWrapper>(b);
406 auto neg_one = gko::initialize<gko_dense>({-1.0}, exec);
407 res->add_scaled(neg_one, b);
417 if (implicit_sq_residual_norm)
419 auto dense_norm = gko::as<gko_dense>(implicit_sq_residual_norm);
421 residual_norms.push_back(
get_norm(dense_norm));
424 else if (residual_norm)
426 auto dense_norm = gko::as<gko_dense>(residual_norm);
428 residual_norms.push_back(
get_norm(dense_norm));
433 auto dense_residual = gko::as<gko_dense>(residual);
437 residual_norms.push_back(norm);
441 iterations.push_back(iteration);
445 const gko::LinOp *matrix;
451 mutable std::vector<ValueType> residual_norms{};
453 mutable std::vector<std::size_t> iterations{};
456 const bool compute_real_residual;
526 return this->executor;
530 std::shared_ptr<gko::Executor> executor;
733 std::shared_ptr<gko::stop::ResidualNorm<>::Factory>
740 std::shared_ptr<gko::stop::ResidualNorm<>::Factory>
747 std::shared_ptr<gko::stop::ImplicitResidualNorm<>::Factory>
754 std::shared_ptr<gko::stop::ImplicitResidualNorm<>::Factory>
805 initialize_ginkgo_log(gko::matrix::Dense<real_t>*
b)
const;
811 std::shared_ptr<gko::LinOp> system_oper;
821template<
typename SolverType>
832 auto current_params = gko::as<typename SolverType::Factory>
846 auto current_params = gko::as<typename SolverType::Factory>
860 auto current_params = gko::as<typename SolverType::Factory>
1026using gko::solver::cb_gmres::storage_precision;
1057 storage_precision prec = storage_precision::reduce1);
1076 storage_precision prec = storage_precision::reduce1);
1137 const std::string &storage_opt =
"none",
1138 const real_t accuracy = 1.e-1,
1139 const int max_block_size = 32
1171 const std::string &factorization_type =
"exact",
1172 const int sweeps = 0,
1173 const bool skip_sort =
false
1213 const std::string &factorization_type =
"exact",
1214 const int sweeps = 0,
1215 const int sparsity_power = 1,
1216 const bool skip_sort =
false
1248 const std::string &factorization_type =
"exact",
1249 const int sweeps = 0,
1250 const bool skip_sort =
false
1290 const std::string &factorization_type =
"exact",
1291 const int sweeps = 0,
1292 const int sparsity_power = 1,
1293 const bool skip_sort =
false
1313 const Solver &mfem_precond
1323 MFEM_ABORT(
"Ginkgo::MFEMPreconditioner must be constructed "
1324 "with the MFEM Operator that it will wrap as an argument;\n"
1325 "calling SetOperator() is not allowed.");
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
BICGSTABSolver(GinkgoExecutor &exec)
CBGMRESSolver(GinkgoExecutor &exec, int dim=0, storage_precision prec=storage_precision::reduce1)
CGSSolver(GinkgoExecutor &exec)
CGSolver(GinkgoExecutor &exec)
void SetRelTol(real_t rtol)
void SetAbsTol(real_t atol)
void SetMaxIter(int max_it)
EnableGinkgoSolver(GinkgoExecutor &exec, bool use_implicit_res_norm)
FCGSolver(GinkgoExecutor &exec)
GMRESSolver(GinkgoExecutor &exec, int dim=0)
std::shared_ptr< gko::Executor > GetExecutor() const
GinkgoExecutor(ExecType exec_type)
virtual ~GinkgoExecutor()=default
@ REFERENCE
Reference CPU Executor.
@ OMP
OpenMP CPU Executor.
void SetOperator(const Operator &op) override
bool sub_op_needs_wrapped_vecs
bool UsesVectorWrappers() const
std::shared_ptr< gko::stop::ImplicitResidualNorm<>::Factory > imp_abs_criterion
void Mult(const Vector &x, Vector &y) const override
std::shared_ptr< ResidualLogger<> > residual_logger
void update_stop_factory()
bool use_implicit_res_norm
std::shared_ptr< gko::stop::ResidualNorm<>::Factory > abs_criterion
std::shared_ptr< gko::stop::Combined::Factory > combined_factory
real_t GetFinalNorm() const
int GetNumIterations() const
const std::shared_ptr< gko::LinOpFactory > GetFactory() const
std::shared_ptr< gko::log::Convergence<> > convergence_logger
GinkgoIterativeSolver(GinkgoExecutor &exec, bool use_implicit_res_norm)
void SetPrintLevel(int print_lvl)
std::shared_ptr< gko::LinOpFactory > solver_gen
std::shared_ptr< gko::stop::ResidualNorm<>::Factory > rel_criterion
virtual ~GinkgoIterativeSolver()=default
std::shared_ptr< gko::stop::ImplicitResidualNorm<>::Factory > imp_rel_criterion
std::shared_ptr< gko::Executor > executor
std::shared_ptr< gko::LinOp > solver
void SetOperator(const Operator &op) override
std::shared_ptr< gko::Executor > executor
const std::shared_ptr< gko::LinOpFactory > GetFactory() const
void Mult(const Vector &x, Vector &y) const override
virtual ~GinkgoPreconditioner()=default
std::shared_ptr< gko::LinOp > generated_precond
bool HasGeneratedPreconditioner() const
bool has_generated_precond
GinkgoPreconditioner(GinkgoExecutor &exec)
std::shared_ptr< gko::LinOpFactory > precond_gen
const std::shared_ptr< gko::LinOp > GetGeneratedPreconditioner() const
IRSolver(GinkgoExecutor &exec)
IcIsaiPreconditioner(GinkgoExecutor &exec, const std::string &factorization_type="exact", const int sweeps=0, const int sparsity_power=1, const bool skip_sort=false)
IcPreconditioner(GinkgoExecutor &exec, const std::string &factorization_type="exact", const int sweeps=0, const bool skip_sort=false)
IluIsaiPreconditioner(GinkgoExecutor &exec, const std::string &factorization_type="exact", const int sweeps=0, const int sparsity_power=1, const bool skip_sort=false)
IluPreconditioner(GinkgoExecutor &exec, const std::string &factorization_type="exact", const int sweeps=0, const bool skip_sort=false)
JacobiPreconditioner(GinkgoExecutor &exec, const std::string &storage_opt="none", const real_t accuracy=1.e-1, const int max_block_size=32)
MFEMPreconditioner(GinkgoExecutor &exec, const Solver &mfem_precond)
void SetOperator(const Operator &op) override
OperatorWrapper(std::shared_ptr< const gko::Executor > exec, gko::size_type size=0, const Operator *oper=NULL)
void apply_impl(const gko::LinOp *b, gko::LinOp *x) const override
static std::unique_ptr< VectorWrapper > create(std::shared_ptr< const gko::Executor > exec, gko::size_type size, Vector *mfem_vec, bool ownership=false)
Vector & get_mfem_vec_ref()
std::unique_ptr< gko::matrix::Dense< real_t > > create_with_same_config() const override
std::unique_ptr< gko::matrix::Dense< real_t > > create_submatrix_impl(const gko::span &rows, const gko::span &columns, const gko::size_type stride) override
std::unique_ptr< gko::matrix::Dense< real_t > > create_with_type_of_impl(std::shared_ptr< const gko::Executor > exec, const gko::dim< 2 > &size, gko::size_type stride) const override
VectorWrapper(std::shared_ptr< const gko::Executor > exec, gko::size_type size, Vector *mfem_vec, bool ownership=false)
const Vector & get_mfem_vec_const_ref() const
void operator()(pointer ptr) const noexcept
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
real_t get_norm(const gko::matrix::Dense< ValueType > *norm)
gko::array< T > gko_array
real_t compute_norm(const gko::matrix::Dense< ValueType > *b)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
void on_iteration_complete(const gko::LinOp *op, const gko::size_type &iteration, const gko::LinOp *residual, const gko::LinOp *solution, const gko::LinOp *residual_norm) const override
void on_iteration_complete(const gko::LinOp *op, const gko::LinOp *rhs, const gko::LinOp *solution, const gko::size_type &iteration, const gko::LinOp *residual, const gko::LinOp *residual_norm, const gko::LinOp *implicit_sq_residual_norm, const gko::array< gko::stopping_status > *status, bool stopped) const override
ResidualLogger(std::shared_ptr< const gko::Executor > exec, const gko::LinOp *matrix, const gko_dense *b, bool compute_real_residual=false)
gko::matrix::Dense< ValueType > gko_dense
void on_iteration_complete(const gko::LinOp *op, const gko::size_type &iteration, const gko::LinOp *residual, const gko::LinOp *solution, const gko::LinOp *residual_norm, const gko::LinOp *implicit_sq_residual_norm) const override