From: Pratik Nayak Date: Wed, 30 Jan 2019 11:50:06 +0000 (+0100) Subject: Add the ginkgo source and include files. X-Git-Tag: v9.1.0-rc1~387^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=57381dfaaeb088cda88e1628a7336c56dfc6ba18;p=dealii.git Add the ginkgo source and include files. --- diff --git a/include/deal.II/lac/ginkgo_solver.h b/include/deal.II/lac/ginkgo_solver.h new file mode 100644 index 0000000000..3ac50c9128 --- /dev/null +++ b/include/deal.II/lac/ginkgo_solver.h @@ -0,0 +1,243 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_ginkgo_solver_h +# define dealii_ginkgo_solver_h + + +# include + +# ifdef DEAL_II_WITH_GINKGO + +# include +# include +# include +# include +# include + +# include + +# include + +DEAL_II_NAMESPACE_OPEN + +namespace GinkgoWrappers +{ + /** + * This class forms the base class for all of Ginkgo's iterative solvers. + * 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 documentation + * and manual pages. + * + * @ingroup GinkgoWrappers + */ + template + class SolverBase + { + public: + /** + * Constructor. + * + * The @p executor defines the paradigm where the solution is computed. + * Ginkgo currently supports three 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); + * ``` + * auto omp = gko::create(); + * ``` + * + CudaExecutor specifies that the data should be stored and the + * operations executed on the NVIDIA GPU accelerator; + * ``` + * if(gko::CudaExecutor::get_num_devices() > 0 ) { + * auto cuda = gko::create(); + * } + * ``` + * + ReferenceExecutor executes a non-optimized reference implementation, + * which can be used to debug the library. + * ``` + * auto ref = gko::create(); + * ``` + * + * The following code snippet demonstrates the using of the OpenMP executor + * to create a solver which would use the OpenMP paradigm to the solve the + * system on the CPU. + * + * ``` + * auto omp = gko::create(); + * using cg = gko::solver::Cg<>; + * auto solver_gen = + * cg::build() + * .with_criteria( + * gko::stop::Iteration::build().with_max_iters(20u).on(omp), + * gko::stop::ResidualNormReduction<>::build() + * .with_reduction_factor(1e-6) + * .on(omp)) + * .on(omp); + * auto solver = solver_gen->generate(system_matrix); + * + * solver->apply(lend(rhs), lend(solution)); + * ``` + * + * + * The @p solver_control object is the same as for other + * deal.II iterative solvers. + */ + SolverBase(SolverControl & solver_control, + std::shared_ptr executor); + + /** + * Destructor. + */ + virtual ~SolverBase() = default; + + /** + * Initialize the matrix and copy over its data to Ginkgo's data structures. + */ + void + initialize(const SparseMatrix &matrix); + + /** + * Solve the linear system Ax=b. Dependent on the information + * provided by derived classes one of Ginkgo's linear solvers is + * chosen. + */ + void + apply(Vector &solution, const Vector &rhs); + + /** + * Solve the linear system Ax=b. Dependent on the information + * provided by derived classes one of Ginkgo's linear solvers is + * chosen. + */ + void + solve(const SparseMatrix &matrix, + Vector & solution, + const Vector & rhs); + + /** + * Access to the object that controls convergence. + */ + SolverControl & + control() const; + + + protected: + /** + * Reference to the object that controls convergence of the iterative + * solvers. + */ + SolverControl &solver_control; + + /** + * The Ginkgo generated solver factory object. + */ + std::shared_ptr solver_gen; + + /** + * The residual criterion object that controls the reduction of the residual + * based on the tolerance set in the solver_control member. + */ + std::shared_ptr::Factory> + residual_criterion; + + /** + * The Ginkgo convergence logger used to check for convergence and other + * solver data if needed. + */ + std::shared_ptr> convergence_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; + + private: + /** + * Initialize the Ginkgo logger object with event masks. Refer to + * Ginkgo's + * logging event masks. + */ + void + initialize_ginkgo_log(); + + /** + * Ginkgo matrix data structure. First template parameter is for storing the + * array of the non-zeros of the matrix. The second is for the row pointers + * and the column indices. + * + * @todo Templatize based on Matrix type. + */ + std::shared_ptr> system_matrix; + }; + + + /** + * An implementation of the solver interface using the Ginkgo CG solver. + * + * @ingroup GinkgoWrappers + */ + template + class SolverCG : public SolverBase + { + public: + /** + * A standardized data struct to pipe additional data to the solver. + */ + struct AdditionalData + {}; + + /** + * Constructor. + * + * @p solver_control The solver control object is then used to set the + * parameters and setup the CG solver from the CG factory which solves the + * linear system. + * + * @p executor The execution paradigm for the CG solver. + */ + SolverCG(SolverControl & solver_control, + std::shared_ptr executor, + const AdditionalData & data = AdditionalData()); + + protected: + /** + * Store a copy of the settings for this particular solver. + */ + const AdditionalData additional_data; + }; + + +} // namespace GinkgoWrappers + +DEAL_II_NAMESPACE_CLOSE + +# endif // DEAL_II_WITH_GINKGO + +#endif +/*---------------------------- ginkgo_solver.h ---------------------------*/ diff --git a/source/lac/ginkgo_solver.cc b/source/lac/ginkgo_solver.cc new file mode 100644 index 0000000000..b1e08892ee --- /dev/null +++ b/source/lac/ginkgo_solver.cc @@ -0,0 +1,313 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#include +#include + +#include + +#ifdef DEAL_II_WITH_GINKGO + +# include + +# include + + +DEAL_II_NAMESPACE_OPEN + +namespace GinkgoWrappers +{ + template + SolverBase::SolverBase( + SolverControl & solver_control, + std::shared_ptr executor) + : solver_control(solver_control) + , executor(executor) + { + using ResidualCriterionFactory = gko::stop::ResidualNormReduction<>; + residual_criterion = ResidualCriterionFactory::build() + .with_reduction_factor(solver_control.tolerance()) + .on(executor); + + combined_factory = + gko::stop::Combined::build() + .with_criteria(residual_criterion, + gko::stop::Iteration::build() + .with_max_iters(solver_control.max_steps()) + .on(executor)) + .on(executor); + } + + template + void + SolverBase::initialize_ginkgo_log() + { + // Add the logger object. See the different masks available in Ginkgo's + // documentation + convergence_logger = gko::log::Convergence<>::create( + executor, gko::log::Logger::criterion_check_completed_mask); + } + + template + void + SolverBase::apply(Vector & solution, + const Vector &rhs) + { + // some shortcuts. + using val_array = gko::Array; + using vec = gko::matrix::Dense; + + Assert(system_matrix, ExcNotInitialized()); + Assert(executor, ExcNotInitialized()); + Assert(rhs.size() == solution.size(), + ExcDimensionMismatch(rhs.size(), solution.size())); + + // Generate the solver from the solver using the system matrix. + auto solver = solver_gen->generate(system_matrix); + + // Create the rhs vector in Ginkgo's format. + std::vector f(rhs.size()); + std::copy(rhs.begin(), rhs.begin() + rhs.size(), f.begin()); + auto b = + vec::create(executor, + gko::dim<2>(rhs.size(), 1), + val_array::view(executor->get_master(), rhs.size(), f.data()), + 1); + + // Create the solution vector in Ginkgo's format. + std::vector u(solution.size()); + std::copy(solution.begin(), solution.begin() + solution.size(), u.begin()); + auto x = vec::create(executor, + gko::dim<2>(solution.size(), 1), + val_array::view(executor->get_master(), + solution.size(), + u.data()), + 1); + + // Create the logger object to log some data from the solvers to confirm + // convergence. + initialize_ginkgo_log(); + + Assert(convergence_logger, ExcNotInitialized()); + // Add the convergence logger object to the combined factory to retrieve the + // solver and other data + combined_factory->add_logger(convergence_logger); + + // Finally, apply the solver to b and get the solution in x. + solver->apply(gko::lend(b), gko::lend(x)); + + // The convergence_logger object contains the residual vector after the + // solver has returned. use this vector to compute the residual norm of the + // solution. Get the residual norm from the logger. As the convergence + // logger returns a `linop`, it is necessary to convert it to a Dense + // matrix. Additionally, if the logger is logging on the gpu, it is + // necessary to copy the data to the host and hence the + // `residual_norm_d_master` + auto residual_norm = convergence_logger->get_residual_norm(); + auto residual_norm_d = + gko::as>(residual_norm); + auto residual_norm_d_master = + gko::matrix::Dense::create(executor->get_master(), + gko::dim<2>{1, 1}); + residual_norm_d_master->copy_from(residual_norm_d); + + // Get the number of iterations taken to converge to the solution. + auto num_iteration = convergence_logger->get_num_iterations(); + + // Ginkgo works with a relative residual norm through its + // ResidualNormReduction criterion. Therefore, to get the normalized + // residual, we divide by the norm of the rhs. + auto b_norm = gko::matrix::Dense::create(executor->get_master(), + gko::dim<2>{1, 1}); + if (executor != executor->get_master()) + { + auto b_master = vec::create(executor->get_master(), + gko::dim<2>(rhs.size(), 1), + val_array::view(executor->get_master(), + rhs.size(), + f.data()), + 1); + b_master->compute_norm2(b_norm.get()); + } + else + { + b->compute_norm2(b_norm.get()); + } + + Assert(b_norm.get()->at(0, 0) != 0.0, ExcDivideByZero()); + // Pass the number of iterations and residual norm to the solver_control + // object. As both `residual_norm_d_master` and `b_norm` are seen as Dense + // matrices, we use the `at` function to get the first value here. In case + // of multiple right hand sides, this will need to be modified. + const SolverControl::State state = + solver_control.check(num_iteration, + residual_norm_d_master->at(0, 0) / b_norm->at(0, 0)); + + // in case of failure: throw exception + if (state != SolverControl::success) + AssertThrow(false, + SolverControl::NoConvergence(solver_control.last_step(), + solver_control.last_value())); + + // Check if the solution is on a CUDA device, if so, copy it over to the + // host. + if (executor != executor->get_master()) + { + auto x_master = vec::create(executor->get_master(), + gko::dim<2>(solution.size(), 1), + val_array::view(executor, + solution.size(), + x->get_values()), + 1); + x.reset(x_master.release()); + } + // Finally copy over the solution vector to deal.II's solution vector. + std::copy(x->get_values(), + x->get_values() + solution.size(), + solution.begin()); + } + + + template + SolverControl & + SolverBase::control() const + { + return solver_control; + } + + + template + void + SolverBase::initialize( + const SparseMatrix &matrix) + { + // Needs to be a square matrix + Assert(matrix.m() == matrix.n(), ExcNotQuadratic()); + + using size_type = dealii::types::global_dof_index; + const size_type N = matrix.m(); + + using mtx = gko::matrix::Csr; + std::shared_ptr system_matrix_compute; + system_matrix_compute = mtx::create(executor->get_master(), + gko::dim<2>(N), + matrix.n_nonzero_elements()); + ValueType *mat_values = system_matrix_compute->get_values(); + IndexType *mat_row_ptrs = system_matrix_compute->get_row_ptrs(); + IndexType *mat_col_idxs = system_matrix_compute->get_col_idxs(); + + // Copy over the data from the matrix to the data structures Ginkgo needs. + // + // Final note: if the matrix has entries in the sparsity pattern that are + // actually occupied by entries that have a zero numerical value, then we + // keep them anyway. people are supposed to provide accurate sparsity + // patterns. + + // first fill row lengths array + mat_row_ptrs[0] = 0; + for (size_type row = 1; row <= N; ++row) + mat_row_ptrs[row] = + mat_row_ptrs[row - 1] + matrix.get_row_length(row - 1); + + // Copy over matrix elements. note that for sparse matrices, + // iterators are sorted so that they traverse each row from start to end + // before moving on to the next row. however, this isn't true for block + // matrices, so we have to do a bit of book keeping + { + // Have an array that for each row points to the first entry not yet + // written to + std::vector row_pointers(N + 1); + std::copy(system_matrix_compute->get_row_ptrs(), + system_matrix_compute->get_row_ptrs() + N + 1, + row_pointers.begin()); + + // Loop over the elements of the matrix row by row, as suggested in the + // documentation of the sparse matrix iterator class + for (size_type row = 0; row < N; ++row) + { + for (typename SparseMatrix::const_iterator p = + matrix.begin(row); + p != matrix.end(row); + ++p) + { + // Write entry into the first free one for this row + mat_col_idxs[row_pointers[row]] = p->column(); + mat_values[row_pointers[row]] = p->value(); + + // Then move pointer ahead + ++row_pointers[row]; + } + } + + // At the end, we should have written all rows completely + for (size_type i = 0; i < N - 1; ++i) + Assert(row_pointers[i] == mat_row_ptrs[i + 1], ExcInternalError()); + } + system_matrix = + mtx::create(executor, gko::dim<2>(N), matrix.n_nonzero_elements()); + system_matrix->copy_from(system_matrix_compute.get()); + } + + + template + void + SolverBase::solve(const SparseMatrix &matrix, + Vector & solution, + const Vector &rhs) + { + initialize(matrix); + apply(solution, rhs); + } + + + /* ---------------------- SolverCG ------------------------ */ + + template + SolverCG::SolverCG( + SolverControl & solver_control, + std::shared_ptr executor, + const AdditionalData & data) + : SolverBase(solver_control, executor) + , additional_data(data) + { + using cg = gko::solver::Cg; + this->solver_gen = + cg::build().with_criteria(this->combined_factory).on(executor); + } + + // Explicit instantiations in GinkgoWrappers +# define GKO_DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(_macro) \ + template _macro(float, int32_t); \ + template _macro(double, int32_t); \ + template _macro(float, int64_t); \ + template _macro(double, int64_t); + +# define DECLARE_SOLVER_BASE(ValueType, IndexType) \ + class SolverBase + GKO_DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(DECLARE_SOLVER_BASE); +# undef DECLARE_SOLVER_BASE + +# define DECLARE_SOLVER_CG(ValueType, IndexType) \ + class SolverCG + GKO_DEALII_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(DECLARE_SOLVER_CG); +# undef DECLARE_SOLVER_CG + +} // namespace GinkgoWrappers + + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_GINKGO