// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced // at the Lawrence Livermore National Laboratory. All Rights reserved. See files // LICENSE and NOTICE for details. LLNL-CODE-806117. // // This file is part of the MFEM library. For more information and source code // availability visit https://mfem.org. // // MFEM is free software; you can redistribute it and/or modify it under the // terms of the BSD-3 license. We welcome feedback and contributions, see file // CONTRIBUTING.md for details. #ifndef MFEM_GINKGO #define MFEM_GINKGO #include "../config/config.hpp" #ifdef MFEM_USE_GINKGO #include "operator.hpp" #include "sparsemat.hpp" #include "solvers.hpp" #include #include #include #include #include #include #include #define MFEM_GINKGO_VERSION \ ((GKO_VERSION_MAJOR*100 + GKO_VERSION_MINOR)*100 + GKO_VERSION_PATCH) namespace mfem { namespace Ginkgo { /// The alias 'gko_array' refers to 'gko::Array' (Ginkgo < 1.5.0) or /// 'gko::array' (Ginkgo >= 1.5.0). #if MFEM_GINKGO_VERSION < 10500 template using gko_array = gko::Array; #else template using gko_array = gko::array; #endif /** * Helper class for a case where a wrapped MFEM Vector * should be owned by Ginkgo, and deleted when the wrapper * object goes out of scope. */ template class gko_mfem_destroy { public: using pointer = T *; // Destroys an MFEM object. Requires object to have a Destroy() method. void operator()(pointer ptr) const noexcept { ptr->Destroy(); } }; /** * This class wraps an MFEM vector object for Ginkgo's use. It * allows Ginkgo and MFEM to operate directly on the same * data, and is necessary to use MFEM Operators with Ginkgo * solvers. * * @ingroup Ginkgo */ class VectorWrapper : public gko::matrix::Dense { public: VectorWrapper(std::shared_ptr exec, gko::size_type size, Vector *mfem_vec, bool ownership = false) : gko::matrix::Dense( exec, gko::dim<2> {size, 1}, gko_array::view(exec, size, mfem_vec->ReadWrite( exec != exec->get_master() ? true : false)), 1) { // This controls whether or not we want Ginkgo to own its MFEM Vector. // Normally, when we are wrapping an MFEM Vector created outside // Ginkgo, we do not want ownership to be true. However, Ginkgo // creates its own temporary vectors as part of its solvers, and // these will be owned (and deleted) by Ginkgo. if (ownership) { using deleter = gko_mfem_destroy; wrapped_vec = std::unique_ptr>( mfem_vec, deleter{}); } else { using deleter = gko::null_deleter; wrapped_vec = std::unique_ptr>( mfem_vec, deleter{}); } } static std::unique_ptr create( std::shared_ptr exec, gko::size_type size, Vector *mfem_vec, bool ownership = false) { return std::unique_ptr( new VectorWrapper(exec, size, mfem_vec, ownership)); } // Return reference to MFEM Vector object Vector &get_mfem_vec_ref() { return *(this->wrapped_vec.get()); } // Return const reference to MFEM Vector object const Vector &get_mfem_vec_const_ref() const { return *(this->wrapped_vec.get()); } // Override base Dense class implementation for creating new vectors // with same executor and size as self virtual std::unique_ptr> create_with_same_config() const override { Vector *mfem_vec = new Vector( this->get_size()[0], this->wrapped_vec.get()->GetMemory().GetMemoryType()); mfem_vec->UseDevice(this->wrapped_vec.get()->UseDevice()); // If this function is called, Ginkgo is creating this // object and should control the memory, so ownership is // set to true return VectorWrapper::create(this->get_executor(), this->get_size()[0], mfem_vec, true); } // Override base Dense class implementation for creating new vectors // with same executor and type as self, but with a different size. // This function will create "one large VectorWrapper" of size // size[0] * size[1], since MFEM Vectors only have one dimension. virtual std::unique_ptr> create_with_type_of_impl( std::shared_ptr exec, const gko::dim<2> &size, gko::size_type stride) const override { // Only stride of 1 is allowed for VectorWrapper type if (stride > 1) { throw gko::Error( __FILE__, __LINE__, "VectorWrapper cannot be created with stride > 1"); } // Compute total size of new Vector gko::size_type total_size = size[0]*size[1]; Vector *mfem_vec = new Vector( total_size, this->wrapped_vec.get()->GetMemory().GetMemoryType()); mfem_vec->UseDevice(this->wrapped_vec.get()->UseDevice()); // If this function is called, Ginkgo is creating this // object and should control the memory, so ownership is // set to true return VectorWrapper::create( this->get_executor(), total_size, mfem_vec, true); } // Override base Dense class implementation for creating new sub-vectors // from a larger vector. virtual std::unique_ptr> create_submatrix_impl( const gko::span &rows, const gko::span &columns, const gko::size_type stride) override { gko::size_type num_rows = rows.end - rows.begin; gko::size_type num_cols = columns.end - columns.begin; // Data in the Dense matrix will be stored in row-major format. // Check that we only have one column, and that the stride = 1 // (only allowed value for VectorWrappers). if (num_cols > 1 || stride > 1) { throw gko::BadDimension( __FILE__, __LINE__, __func__, "new_submatrix", num_rows, num_cols, "VectorWrapper submatrix must have one column and stride = 1"); } int data_size = static_cast(num_rows * num_cols); int start = static_cast(rows.begin); // Create a new MFEM Vector pointing to this starting point in the data Vector *mfem_vec = new Vector(); mfem_vec->MakeRef(*(this->wrapped_vec.get()), start, data_size); mfem_vec->UseDevice(this->wrapped_vec.get()->UseDevice()); // If this function is called, Ginkgo is creating this // object and should control the memory, so ownership is // set to true (but MFEM doesn't own and won't delete // the data, at it's only a reference to the parent Vector) return VectorWrapper::create( this->get_executor(), data_size, mfem_vec, true); } private: std::unique_ptr> wrapped_vec; }; /** * This class wraps an MFEM Operator for Ginkgo, to make its Mult() * function available to Ginkgo, provided the input and output vectors * are of the VectorWrapper type. * Note that this class does NOT take ownership of the MFEM Operator. * * @ingroup Ginkgo */ class OperatorWrapper : public gko::EnableLinOp, public gko::EnableCreateMethod { public: OperatorWrapper(std::shared_ptr exec, gko::size_type size = 0, const Operator *oper = NULL) : gko::EnableLinOp(exec, gko::dim<2> {size, size}), gko::EnableCreateMethod() { this->wrapped_oper = oper; } protected: void apply_impl(const gko::LinOp *b, gko::LinOp *x) const override; void apply_impl(const gko::LinOp *alpha, const gko::LinOp *b, const gko::LinOp *beta, gko::LinOp *x) const override; private: const Operator *wrapped_oper; }; // Utility function which gets the scalar value of a Ginkgo gko::matrix::Dense // matrix representing the norm of a vector. template double get_norm(const gko::matrix::Dense *norm) { // Put the value on CPU thanks to the master executor auto cpu_norm = clone(norm->get_executor()->get_master(), norm); // Return the scalar value contained at position (0, 0) return cpu_norm->at(0, 0); } // Utility function which computes the norm of a Ginkgo gko::matrix::Dense // vector. template double compute_norm(const gko::matrix::Dense *b) { // Get the executor of the vector auto exec = b->get_executor(); // Initialize a result scalar containing the value 0.0. auto b_norm = gko::initialize>({0.0}, exec); // Use the dense `compute_norm2` function to compute the norm. b->compute_norm2(lend(b_norm)); // Use the other utility function to return the norm contained in `b_norm`` return std::pow(get_norm(lend(b_norm)),2); } /** * Custom logger class which intercepts the residual norm scalar and solution * vector in order to print a table of real vs recurrent (internal to the * solvers) residual norms. * * This has been taken from the custom-logger example of Ginkgo. See the * custom-logger example to understand how to write and modify your own loggers * with Ginkgo. * * @ingroup Ginkgo */ template struct ResidualLogger : gko::log::Logger { // Output the logger's data in a table format void write() const { // Print a header for the table if (compute_real_residual) { mfem::out << "Iteration log with real residual norms:" << std::endl; } else { mfem::out << "Iteration log with residual norms:" << std::endl; } mfem::out << '|' << std::setw(10) << "Iteration" << '|' << std::setw(25) << "Residual Norm" << '|' << std::endl; // Print a separation line. Note that for creating `10` characters // `std::setw()` should be set to `11`. mfem::out << '|' << std::setfill('-') << std::setw(11) << '|' << std::setw(26) << '|' << std::setfill(' ') << std::endl; // Print the data one by one in the form mfem::out << std::scientific; for (std::size_t i = 0; i < iterations.size(); i++) { mfem::out << '|' << std::setw(10) << iterations[i] << '|' << std::setw(25) << residual_norms[i] << '|' << std::endl; } // std::defaultfloat could be used here but some compilers do not support // it properly, e.g. the Intel compiler mfem::out.unsetf(std::ios_base::floatfield); // Print a separation line mfem::out << '|' << std::setfill('-') << std::setw(11) << '|' << std::setw(26) << '|' << std::setfill(' ') << std::endl; } using gko_dense = gko::matrix::Dense; // Customize the logging hook which is called every time an iteration is // completed void on_iteration_complete(const gko::LinOp *, 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 { // If the solver shares the current solution vector and we want to // compute the residual from that if (solution && compute_real_residual) { // Store the matrix's executor auto exec = matrix->get_executor(); // Compute the real residual vector by calling apply on the system // First, compute res = A * x matrix->apply(gko::lend(solution), gko::lend(res)); // Now do res = res - b, depending on which vector/oper type // Check if b is a Ginkgo vector or wrapped MFEM Vector if (dynamic_cast(b)) { const VectorWrapper *b_cast = gko::as(b); // Copy the MFEM Vector stored in b VectorWrapper *res_cast = gko::as(res); res_cast->get_mfem_vec_ref() -= b_cast->get_mfem_vec_const_ref(); } else { // Create a scalar containing the value -1.0 auto neg_one = gko::initialize({-1.0}, exec); res->add_scaled(gko::lend(neg_one), gko::lend(b)); } // Compute the norm of the residual vector and add it to the // `residual_norms` vector residual_norms.push_back(compute_norm(gko::lend(res))); } else { // If the solver shares an implicit or recurrent residual norm, log its value if (implicit_sq_residual_norm) { auto dense_norm = gko::as(implicit_sq_residual_norm); // Add the norm to the `residual_norms` vector residual_norms.push_back(get_norm(dense_norm)); // Otherwise, use the recurrent residual vector } else if (residual_norm) { auto dense_norm = gko::as(residual_norm); // Add the norm to the `residual_norms` vector residual_norms.push_back(get_norm(dense_norm)); // Otherwise, use the recurrent residual vector } else { auto dense_residual = gko::as(residual); // Compute the residual vector's norm auto norm = compute_norm(gko::lend(dense_residual)); // Add the computed norm to the `residual_norms` vector residual_norms.push_back(norm); } } // Add the current iteration number to the `iterations` vector iterations.push_back(iteration); } // Version for solver that doesn't log implicit res norm 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 { on_iteration_complete(op, iteration, residual, solution, residual_norm, nullptr); } // Construct the logger and store the system matrix and b vectors ResidualLogger(std::shared_ptr exec, const gko::LinOp *matrix, const gko_dense *b, bool compute_real_residual=false) : #if MFEM_GINKGO_VERSION < 10500 gko::log::Logger(exec, gko::log::Logger::iteration_complete_mask), #else gko::log::Logger(gko::log::Logger::iteration_complete_mask), #endif matrix {matrix}, b{b}, compute_real_residual{compute_real_residual} { if (compute_real_residual == true) { if (dynamic_cast(b)) { const VectorWrapper *b_cast = gko::as(b); res = std::move(gko_dense::create_with_config_of(b_cast).release()); } else { res = std::move(gko::clone(b).release()); } } } private: // Pointer to the system matrix const gko::LinOp *matrix; // Pointer to the right hand sides const gko_dense *b; // Pointer to the residual workspace vector gko_dense *res; // Vector which stores all the residual norms mutable std::vector residual_norms{}; // Vector which stores all the iteration numbers mutable std::vector iterations{}; // Whether or not to compute the residual at every iteration, // rather than using the recurrent norm const bool compute_real_residual; }; /** * This class wraps a Ginkgo Executor for use in MFEM. * Note that objects in the Ginkgo namespace intended to work * together, e.g. a Ginkgo solver and preconditioner, should use the same * GinkgoExecutor object. In general, most users will want to create * one GinkgoExecutor object for use with all Ginkgo-related objects. * The wrapper can be created to match MFEM's device configuration. */ class GinkgoExecutor { public: // Types of Ginkgo Executors. enum ExecType { /// Reference CPU Executor. REFERENCE = 0, /// OpenMP CPU Executor. OMP = 1, /// CUDA GPU Executor. CUDA = 2, /// HIP GPU Executor. HIP = 3 }; /** * Constructor. * Takes an @p GinkgoExecType argument and creates an Executor. */ GinkgoExecutor(ExecType exec_type); /** * Constructor. * Takes an MFEM @p Device object and creates an Executor * that "matches" (e.g., if MFEM is using the CPU, Ginkgo * will choose the OmpExecutor; if MFEM is using CUDA, * Ginkgo will choose the CudaExecutor). */ GinkgoExecutor(Device &mfem_device); /** * Destructor. */ virtual ~GinkgoExecutor() = default; std::shared_ptr GetExecutor() const { return this->executor; }; private: std::shared_ptr executor; }; /** * This class forms the base class for all of Ginkgo's preconditioners. The * various derived classes only take the additional data that is specific to them. * The entire collection of preconditioners that Ginkgo implements is available * at the Ginkgo documentation and manual pages, * https://ginkgo-project.github.io/ginkgo/doc/develop. * * @ingroup Ginkgo */ class GinkgoPreconditioner : public Solver { public: /** * Constructor. * * The @p exec defines the paradigm where the solution is computed. * Ginkgo currently supports four different executor types: * * + OmpExecutor specifies that the data should be stored and the * associated operations executed on an OpenMP-supporting device (e.g. * host CPU); * + CudaExecutor specifies that the data should be stored and the * operations executed on the NVIDIA GPU accelerator; * + HipExecutor specifies that the data should be stored and the * operations executed on the GPU accelerator using HIP; * + ReferenceExecutor executes a non-optimized reference implementation, * which can be used to debug the library. */ GinkgoPreconditioner(GinkgoExecutor &exec); /** * Destructor. */ virtual ~GinkgoPreconditioner() = default; /** * Generate the preconditioner for the given matrix @p op, * which must be of MFEM SparseMatrix type. * Calling this function is only required when creating a * preconditioner for use with another MFEM solver; to use with * a Ginkgo solver, get the LinOpFactory pointer through @p GetFactory() * and pass to the Ginkgo solver constructor. */ virtual void SetOperator(const Operator &op); /** * Apply the preconditioner to input vector @p x, with out @p y. */ virtual void Mult(const Vector &x, Vector &y) const; /** * Return a pointer to the LinOpFactory that will generate the preconditioner * with the parameters set through the specific constructor. */ const std::shared_ptr GetFactory() const { return this->precond_gen; }; /** * Return a pointer to the generated preconditioner for a specific matrix * (that has previously been set with @p SetOperator). */ const std::shared_ptr GetGeneratedPreconditioner() const { return this->generated_precond; }; /** * Return whether this GinkgoPreconditioner object has an explicitly- * generated preconditioner, built for a specific matrix. */ bool HasGeneratedPreconditioner() const { return this->has_generated_precond; }; protected: /** * The Ginkgo generated solver factory object. */ std::shared_ptr precond_gen; /** * Generated Ginkgo preconditioner for a specific matrix, created through * @p SetOperator(), or a wrapped MFEM preconditioner. * Must exist to use @p Mult(). */ std::shared_ptr generated_precond; /** * The execution paradigm in Ginkgo. The choices are between * `gko::OmpExecutor`, `gko::CudaExecutor` and `gko::ReferenceExecutor` * and more details can be found in Ginkgo's documentation. */ std::shared_ptr executor; /** * Whether or not we have generated a specific preconditioner for * a matrix. */ bool has_generated_precond; }; /** * This class forms the base class for all of Ginkgo's iterative solvers. * It is not intended to be used directly by MFEM applications. The various * derived classes only take the additional data that is specific to them * and solve the given linear system. The entire collection of solvers that * Ginkgo implements is available at the Ginkgo documentation and manual pages, * https://ginkgo-project.github.io/ginkgo/doc/develop. * * @ingroup Ginkgo */ class GinkgoIterativeSolver : public Solver { public: /** * Return a pointer to the LinOpFactory that will generate the solver * with the parameters set through the specific constructor. */ const std::shared_ptr GetFactory() const { return this->solver_gen; }; void SetPrintLevel(int print_lvl) { print_level = print_lvl; } int GetNumIterations() const { return final_iter; } int GetConverged() const { return converged; } double GetFinalNorm() const { return final_norm; } /** * If the Operator is a SparseMatrix, set up a Ginkgo Csr matrix * to use its data directly. If the Operator is not a matrix, * create an OperatorWrapper for it and store. */ virtual void SetOperator(const Operator &op); /** * Solve the linear system Ax=y. Dependent on the information * provided by derived classes one of Ginkgo's linear solvers is chosen. */ virtual void Mult(const Vector &x, Vector &y) const; /** * Return whether this GinkgoIterativeSolver object will use * VectorWrapper types for input and output vectors. * Note that Mult() will automatically create these wrappers if needed. */ bool UsesVectorWrappers() const { return this->needs_wrapped_vecs; }; /** * Destructor. */ virtual ~GinkgoIterativeSolver() = default; protected: /** * Constructor. * * The @p exec defines the paradigm where the solution is computed. * @p use_implicit_res_norm is for internal use by the derived classes * for specific Ginkgo solvers; it indicates whether the solver makes * an implicit residual norm estimate available for convergence checking. * Each derived class automatically sets the correct value when calling this * base class constructor. * */ GinkgoIterativeSolver(GinkgoExecutor &exec, bool use_implicit_res_norm); bool use_implicit_res_norm; int print_level; int max_iter; double rel_tol; double abs_tol; mutable double final_norm; mutable int final_iter; mutable int converged; /** * The Ginkgo solver factory object, to generate specific solvers. */ std::shared_ptr solver_gen; /** * The Ginkgo solver object, generated for a specific operator. */ std::shared_ptr solver; /** * The residual criterion object that controls the reduction of the residual * relative to the initial residual. */ std::shared_ptr::Factory> rel_criterion; /** * The residual criterion object that controls the reduction of the residual * based on an absolute tolerance. */ std::shared_ptr::Factory> abs_criterion; /** * The implicit residual criterion object that controls the reduction of the residual * relative to the initial residual, based on an implicit residual norm value. */ std::shared_ptr::Factory> imp_rel_criterion; /** * The implicit residual criterion object that controls the reduction of the residual * based on an absolute tolerance, based on an implicit residual norm value. */ std::shared_ptr::Factory> imp_abs_criterion; /** * The Ginkgo convergence logger used to check for convergence and other * solver data if needed. */ mutable std::shared_ptr> convergence_logger; /** * The residual logger object used to check for convergence and other solver * data if needed. */ mutable std::shared_ptr> residual_logger; /** * The Ginkgo combined factory object is used to create a combined stopping * criterion to be passed to the solver. */ std::shared_ptr combined_factory; /** * The execution paradigm in Ginkgo. The choices are between * `gko::OmpExecutor`, `gko::CudaExecutor` and `gko::ReferenceExecutor` * and more details can be found in Ginkgo's documentation. */ std::shared_ptr executor; /** * Whether or not we need to use VectorWrapper types with this solver. */ bool needs_wrapped_vecs; /** * Whether or not we need to use VectorWrapper types with the preconditioner * or an inner solver. This value is set upon creation of the * GinkgoIterativeSolver object and should never change. */ bool sub_op_needs_wrapped_vecs; /** Rebuild the Ginkgo stopping criterion factory with the latest values * of rel_tol, abs_tol, and max_iter. */ void update_stop_factory(); private: /** * Initialize the Ginkgo logger object with event masks. Refer to the logging * event masks in Ginkgo's .../include/ginkgo/core/log/logger.hpp. */ void initialize_ginkgo_log(gko::matrix::Dense* b) const; /** * Pointer to either a Ginkgo CSR matrix or an OperatorWrapper wrapping * an MFEM Operator (for matrix-free evaluation). */ std::shared_ptr system_oper; }; /** * This class adds helper functions for updating Ginkgo factories * and solvers, when the full class type is needed. The derived classes * should inherit from this class, rather than from GinkgoIterativeSolver * directly. */ template class EnableGinkgoSolver : public GinkgoIterativeSolver { public: EnableGinkgoSolver(GinkgoExecutor &exec, bool use_implicit_res_norm) : GinkgoIterativeSolver(exec, use_implicit_res_norm) {} void SetRelTol(double rtol) { rel_tol = rtol; this->update_stop_factory(); gko::as(solver_gen)->get_parameters().criteria = { combined_factory }; if (solver) { gko::as(solver)->set_stop_criterion_factory(combined_factory); } } void SetAbsTol(double atol) { abs_tol = atol; this->update_stop_factory(); gko::as(solver_gen)->get_parameters().criteria = { combined_factory }; if (solver) { gko::as(solver)->set_stop_criterion_factory(combined_factory); } } void SetMaxIter(int max_it) { max_iter = max_it; this->update_stop_factory(); gko::as(solver_gen)->get_parameters().criteria = { combined_factory }; if (solver) { gko::as(solver)->set_stop_criterion_factory(combined_factory); } } }; /** * An implementation of the solver interface using the Ginkgo CG solver. * * @ingroup Ginkgo */ class CGSolver : public EnableGinkgoSolver> { public: /** * Constructor. * * @param[in] exec The execution paradigm for the solver. */ CGSolver(GinkgoExecutor &exec); /** * Constructor. * * @param[in] exec The execution paradigm for the solver. * @param[in] preconditioner The preconditioner for the solver. */ CGSolver(GinkgoExecutor &exec, const GinkgoPreconditioner &preconditioner); }; /** * An implementation of the solver interface using the Ginkgo BiCGStab solver. * * @ingroup Ginkgo */ class BICGSTABSolver : public EnableGinkgoSolver> { public: /** * Constructor. * * @param[in] exec The execution paradigm for the solver. */ BICGSTABSolver(GinkgoExecutor &exec); /** * Constructor. * * @param[in] exec The execution paradigm for the solver. * @param[in] preconditioner The preconditioner for the solver. */ BICGSTABSolver(GinkgoExecutor &exec, const GinkgoPreconditioner &preconditioner); }; /** * An implementation of the solver interface using the Ginkgo CGS solver. * * CGS or the conjugate gradient square method is an iterative type Krylov * subspace method which is suitable for general systems. * * @ingroup Ginkgo */ class CGSSolver : public EnableGinkgoSolver> { public: /** * Constructor. * * @param[in] exec The execution paradigm for the solver. */ CGSSolver(GinkgoExecutor &exec); /** * Constructor. * * @param[in] exec The execution paradigm for the solver. * @param[in] preconditioner The preconditioner for the solver. */ CGSSolver(GinkgoExecutor &exec, const GinkgoPreconditioner &preconditioner); }; /** * An implementation of the solver interface using the Ginkgo FCG solver. * * FCG or the flexible conjugate gradient method is an iterative type Krylov * subspace method which is suitable for symmetric positive definite methods. * * Though this method performs very well for symmetric positive definite * matrices, it is in general not suitable for general matrices. * * In contrast to the standard CG based on the Polack-Ribiere formula, the * flexible CG uses the Fletcher-Reeves formula for creating the orthonormal * vectors spanning the Krylov subspace. This increases the computational cost * of every Krylov solver iteration but allows for non-constant preconditioners. * * @ingroup Ginkgo */ class FCGSolver : public EnableGinkgoSolver> { public: /** * Constructor. * * @param[in] exec The execution paradigm for the solver. */ FCGSolver(GinkgoExecutor &exec); /** * Constructor. * * @param[in] exec The execution paradigm for the solver. * @param[in] preconditioner The preconditioner for the solver. */ FCGSolver(GinkgoExecutor &exec, const GinkgoPreconditioner &preconditioner); }; /** * An implementation of the solver interface using the Ginkgo GMRES solver. * * @ingroup Ginkgo */ class GMRESSolver : public EnableGinkgoSolver> { public: /** * Constructor. * * @param[in] exec The execution paradigm for the solver. * @param[in] dim The Krylov dimension of the solver. Value of 0 will * let Ginkgo use its own internal default value. */ GMRESSolver(GinkgoExecutor &exec, int dim = 0); /** * Constructor. * * @param[in] exec The execution paradigm for the solver. * @param[in] preconditioner The preconditioner for the solver. * @param[in] dim The Krylov dimension of the solver. Value of 0 will * let Ginkgo use its own internal default value. */ GMRESSolver(GinkgoExecutor &exec, const GinkgoPreconditioner &preconditioner, int dim = 0); /** * Change the Krylov dimension of the solver. */ void SetKDim(int dim); protected: int m; // Dimension of Krylov subspace }; using gko::solver::cb_gmres::storage_precision; /** * An implementation of the solver interface using the Ginkgo * Compressed Basis GMRES solver. With CB-GMRES, the Krylov basis * is "compressed" by storing in a lower precision. Currently, computations * are always performed in double precision when using this MFEM integration. * The Ginkgo storage precision options are accessed * through Ginkgo::storage_precision::*. The default choice * is Ginkgo::storage_precision::reduce1, i.e., store in float * instead of double. * * @ingroup Ginkgo */ class CBGMRESSolver : public EnableGinkgoSolver> { public: /** * Constructor. * * @param[in] exec The execution paradigm for the solver. * @param[in] dim The Krylov dimension of the solver. Value of 0 will * let Ginkgo use its own internal default value. * @param[in] prec The storage precision used in the CB-GMRES. Options * are: keep (keep double precision), reduce1 (double -> float), * reduce2 (double -> half), integer (double -> int64), * ireduce1 (double -> int32), ireduce2 (double -> int16). * See Ginkgo documentation for more about the CB-GMRES and * these options. */ CBGMRESSolver(GinkgoExecutor &exec, int dim = 0, storage_precision prec = storage_precision::reduce1); /** * Constructor. * * @param[in] exec The execution paradigm for the solver. * @param[in] preconditioner The preconditioner for the solver. * @param[in] dim The Krylov dimension of the solver. Value of 0 will * let Ginkgo use its own internal default value. * @param[in] prec The storage precision used in the CB-GMRES. Options * are: keep (keep double precision), reduce1 (double -> float), * reduce2 (double -> half), integer (double -> int64), * ireduce1 (double -> int32), ireduce2 (double -> int16). * See Ginkgo documentation for more about the CB-GMRES and * these options. */ CBGMRESSolver(GinkgoExecutor &exec, const GinkgoPreconditioner &preconditioner, int dim = 0, storage_precision prec = storage_precision::reduce1); /** * Change the Krylov dimension of the solver. */ void SetKDim(int dim); protected: int m; // Dimension of Krylov subspace }; /** * An implementation of the solver interface using the Ginkgo IR solver. * * Iterative refinement (IR) is an iterative method that uses another coarse * method to approximate the error of the current solution via the current * residual. * * @ingroup Ginkgo */ class IRSolver : public EnableGinkgoSolver> { public: /** * Constructor. * * @param[in] exec The execution paradigm for the solver. */ IRSolver(GinkgoExecutor &exec); /** * Constructor. * * @param[in] exec The execution paradigm for the solver. * @param[in] inner_solver The inner solver for the main solver. */ IRSolver(GinkgoExecutor &exec, const GinkgoIterativeSolver &inner_solver); }; /** * An implementation of the preconditioner interface using the Ginkgo Jacobi * preconditioner. * * @ingroup Ginkgo */ class JacobiPreconditioner : public GinkgoPreconditioner { public: /** * Constructor. * * @param[in] exec The execution paradigm for the preconditioner. * @param[in] storage_opt The storage optimization parameter. * @param[in] accuracy The relative accuracy for the adaptive version. * @param[in] max_block_size Maximum block size. * See the Ginkgo documentation for more information on these parameters. */ JacobiPreconditioner( GinkgoExecutor &exec, const std::string &storage_opt = "none", const double accuracy = 1.e-1, const int max_block_size = 32 ); }; /** * An implementation of the preconditioner interface using the Ginkgo * Incomplete LU preconditioner (ILU(0)). * * @ingroup Ginkgo */ class IluPreconditioner : public GinkgoPreconditioner { public: /** * Constructor. * * @param[in] exec The execution paradigm for the preconditioner. * @param[in] factorization_type The factorization type: "exact" or * "parilu". * @param[in] sweeps The number of sweeps to do in the ParIlu * factorization algorithm. A value of 0 tells Ginkgo to use its * internal default value. This parameter is ignored in the case * of an exact factorization. * @param[in] skip_sort Only set this to true if the input matrix * that will be used to generate this preconditioner is guaranteed * to be sorted by column. * * Note: The use of this preconditioner will sort any input matrix * given to it, potentially changing the order of the stored values. */ IluPreconditioner( GinkgoExecutor &exec, const std::string &factorization_type = "exact", const int sweeps = 0, const bool skip_sort = false ); }; /** * An implementation of the preconditioner interface using the Ginkgo * Incomplete LU-Incomplete Sparse Approximate Inverse preconditioner. * The Ilu-ISAI preconditioner differs from the Ilu preconditioner in * that Incomplete Sparse Approximate Inverses (ISAIs) are formed * to approximate solving the triangular systems defined by L and U. * When the preconditioner is applied, these ISAI matrices are applied * through matrix-vector multiplication. * * @ingroup Ginkgo */ class IluIsaiPreconditioner : public GinkgoPreconditioner { public: /** * Constructor. * * @param[in] exec The execution paradigm for the preconditioner. * @param[in] factorization_type The factorization type: "exact" or * "parilu". * @param[in] sweeps The number of sweeps to do in the ParIlu * factorization algorithm. A value of 0 tells Ginkgo to use its * internal default value. This parameter is ignored in the case * of an exact factorization. * @param[in] sparsity_power Parameter determining the sparsity pattern of * the ISAI approximations. * @param[in] skip_sort Only set this to true if the input matrix * that will be used to generate this preconditioner is guaranteed * to be sorted by column. * See the Ginkgo documentation for more information on these parameters. * * Note: The use of this preconditioner will sort any input matrix * given to it, potentially changing the order of the stored values. */ IluIsaiPreconditioner( GinkgoExecutor &exec, const std::string &factorization_type = "exact", const int sweeps = 0, const int sparsity_power = 1, const bool skip_sort = false ); }; /** * An implementation of the preconditioner interface using the Ginkgo * Incomplete Cholesky preconditioner (IC(0)). * * @ingroup Ginkgo */ class IcPreconditioner : public GinkgoPreconditioner { public: /** * Constructor. * * @param[in] exec The execution paradigm for the preconditioner. * @param[in] factorization_type The factorization type: "exact" or * "paric". * @param[in] sweeps The number of sweeps to do in the ParIc * factorization algorithm. A value of 0 tells Ginkgo to use its * internal default value. This parameter is ignored in the case * of an exact factorization. * @param[in] skip_sort Only set this to true if the input matrix * that will be used to generate this preconditioner is guaranteed * to be sorted by column. * * Note: The use of this preconditioner will sort any input matrix * given to it, potentially changing the order of the stored values. */ IcPreconditioner( GinkgoExecutor &exec, const std::string &factorization_type = "exact", const int sweeps = 0, const bool skip_sort = false ); }; /** * An implementation of the preconditioner interface using the Ginkgo * Incomplete Cholesky-Incomplete Sparse Approximate Inverse preconditioner. * The Ic-ISAI preconditioner differs from the Ic preconditioner in * that Incomplete Sparse Approximate Inverses (ISAIs) are formed * to approximate solving the triangular systems defined by L and L^T. * When the preconditioner is applied, these ISAI matrices are applied * through matrix-vector multiplication. * * @ingroup Ginkgo */ class IcIsaiPreconditioner : public GinkgoPreconditioner { public: /** * Constructor. * * @param[in] exec The execution paradigm for the preconditioner. * @param[in] factorization_type The factorization type: "exact" or * "paric". * @param[in] sweeps The number of sweeps to do in the ParIc * factorization algorithm. A value of 0 tells Ginkgo to use its * internal default value. This parameter is ignored in the case * of an exact factorization. * @param[in] sparsity_power Parameter determining the sparsity pattern of * the ISAI approximations. * @param[in] skip_sort Only set this to true if the input matrix * that will be used to generate this preconditioner is guaranteed * to be sorted by column. * See the Ginkgo documentation for more information on these parameters. * * Note: The use of this preconditioner will sort any input matrix * given to it, potentially changing the order of the stored values. */ IcIsaiPreconditioner( GinkgoExecutor &exec, const std::string &factorization_type = "exact", const int sweeps = 0, const int sparsity_power = 1, const bool skip_sort = false ); }; /** * A wrapper that allows Ginkgo to use MFEM preconditioners. * * @ingroup Ginkgo */ class MFEMPreconditioner : public GinkgoPreconditioner { public: /** * Constructor. * * @param[in] exec The execution paradigm for the preconditioner. * @param[in] mfem_precond The MFEM Preconditioner to wrap. */ MFEMPreconditioner( GinkgoExecutor &exec, const Solver &mfem_precond ); /** * SetOperator is not allowed for this type of preconditioner; * this function overrides the base class in order to give an * error if SetOperator() is called for this class. */ virtual void SetOperator(const Operator &op) { MFEM_ABORT("Ginkgo::MFEMPreconditioner must be constructed " "with the MFEM Operator that it will wrap as an argument;\n" "calling SetOperator() is not allowed."); }; }; } // namespace Ginkgo } // namespace mfem #endif // MFEM_USE_GINKGO #endif // MFEM_GINKGO