From 332f5e60332b2ff1129a617688226c645077283e Mon Sep 17 00:00:00 2001 From: blegouix Date: Tue, 30 Apr 2024 16:23:15 +0200 Subject: [PATCH 01/18] init --- include/ddc/kernels/splines/matrix.hpp | 84 ++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) diff --git a/include/ddc/kernels/splines/matrix.hpp b/include/ddc/kernels/splines/matrix.hpp index fa661d9ff..2903d64eb 100644 --- a/include/ddc/kernels/splines/matrix.hpp +++ b/include/ddc/kernels/splines/matrix.hpp @@ -21,12 +21,34 @@ class Matrix public: explicit Matrix(const int mat_size) : m_n(mat_size) {} + /// @brief Destruct virtual ~Matrix() = default; + /** + * @brief Get an element of the matrix at indexes i, j. + * + * @param i The row index of the desired element. + * @param j The columns index of the desired element. + * + * @return The value of the element of the matrix. + */ virtual double get_element(int i, int j) const = 0; + /** + * @brief Set an element of the matrix at indexes i, j. + * + * @param i The row index of the setted element. + * @param j The column index of the setted element. + * @param aij The value to set in the element of the matrix. + */ virtual void set_element(int i, int j, double aij) = 0; + /** + * @brief Performs a pre-process operation on the Matrix. + * + * Note: this function should be renamed in the future because the pre- + * process operation is not necessarily a factorization. + */ virtual void factorize() { int const info = factorize_method(); @@ -45,6 +67,11 @@ class Matrix } } + /** + * @brief Solve the linear problem Ax=b inplace. + * + * @param[in, out] b A 1D mdpsan storing a right-hand-side of the problem and receiving the corresponding solution. + */ virtual ddc::DSpan1D solve_inplace(ddc::DSpan1D const b) const { assert(int(b.extent(0)) == m_n); @@ -57,6 +84,13 @@ class Matrix return b; } + /** + * @brief Solve the transposed linear problem A^tx=b inplace. + * + * @param[in, out] b A 1D mdpsan storing a right-hand-side of the problem and receiving the corresponding solution. + * + * @return b + */ virtual ddc::DSpan1D solve_transpose_inplace(ddc::DSpan1D const b) const { assert(int(b.extent(0)) == m_n); @@ -69,6 +103,14 @@ class Matrix return b; } + /** + * @brief Solve the multiple-right-hand-side linear problem Ax=b inplace. + * + * @param[in, out] bx A 2D mdpsan storing the multiple right-hand-sides of the problem and receiving the corresponding solution. Important note: the convention is the reverse of the common matrix one, + * row number is second index and column number the first one. It should be changed in the future. + * + * @return bx + */ virtual ddc::DSpan2D solve_multiple_inplace(ddc::DSpan2D const bx) const { assert(int(bx.extent(1)) == m_n); @@ -81,6 +123,14 @@ class Matrix return bx; } + /** + * @brief Solve the transposed multiple-right-hand-side linear problem A^tx=b inplace. + * + * @param[in, out] bx A 2D mdpsan storing the multiple right-hand-sides of the problem and receiving the corresponding solution. Important note: the convention is the reverse of the common matrix one, + * row number is second index and column number the first one. It should be changed in the future. + * + * @return bx + */ virtual ddc::DSpan2D solve_multiple_transpose_inplace(ddc::DSpan2D const bx) const { assert(int(bx.extent(1)) == m_n); @@ -93,6 +143,14 @@ class Matrix return bx; } + /** + * @brief Solve the multiple-right-hand-side linear problem Ax=b inplace. + * + * @param[in, out] bx A 2D Kokkos::View storing the multiple right-hand-sides of the problem and receiving the corresponding solution. Important note: the convention is the reverse of the common matrix one, + * row number is second index and column number the first one. It should be changed in the future. + * + * @return bx + */ template Kokkos::View solve_batch_inplace( Kokkos::View const bx) const @@ -107,11 +165,23 @@ class Matrix return bx; } + /** + * @brief Get the size of the matrix (which is necessarily squared) in one of its dimensions. + * + * @return The size of the matrix in one of its dimensions. + */ int get_size() const { return m_n; } + /** + * @brief Prints a Matrix in a std::ostream. It will segfault is the Matrix is on GPU. + * + * @param out The stream in which the matrix is printed. + * + * @return The stream in which the matrix is printed. + **/ std::ostream& operator<<(std::ostream& os) { int const n = get_size(); @@ -125,8 +195,22 @@ class Matrix } protected: + /** + * @brief A function called by factorize() to actually perform the pre-process operation. + * + * @return The error code of the function. + */ virtual int factorize_method() = 0; + /** + * @brief A function called by factorize() to actually perform the pre-process operation. + * + * @param b A double* to a contiguous array containing the (eventually multiple) right-hand-sides. Memory layout depends on the + * derivated class (and the underlying algorithm). + * @param transpose A character identifying if the normal ('N') or transposed ('T') linear system is solved. + * @param n_equations The number of multiple-right-hand-sides (number of columns of b). + * @return The error code of the function. + */ virtual int solve_inplace_method(double* b, char transpose, int n_equations) const = 0; }; From 832a594f8ae01079447754e7a81fea4a35db6f29 Mon Sep 17 00:00:00 2001 From: blegouix Date: Tue, 30 Apr 2024 16:25:08 +0200 Subject: [PATCH 02/18] clang-format --- include/ddc/kernels/splines/matrix.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/ddc/kernels/splines/matrix.hpp b/include/ddc/kernels/splines/matrix.hpp index 2903d64eb..80b5c5f8b 100644 --- a/include/ddc/kernels/splines/matrix.hpp +++ b/include/ddc/kernels/splines/matrix.hpp @@ -175,7 +175,7 @@ class Matrix return m_n; } - /** + /** * @brief Prints a Matrix in a std::ostream. It will segfault is the Matrix is on GPU. * * @param out The stream in which the matrix is printed. From f669f07d28bd12c48210b0a630b8d86dd0e8cb1b Mon Sep 17 00:00:00 2001 From: blegouix Date: Thu, 2 May 2024 11:40:59 +0200 Subject: [PATCH 03/18] wip --- include/ddc/kernels/splines/matrix.hpp | 7 ++- include/ddc/kernels/splines/matrix_sparse.hpp | 57 ++++++++++++++++++- 2 files changed, 61 insertions(+), 3 deletions(-) diff --git a/include/ddc/kernels/splines/matrix.hpp b/include/ddc/kernels/splines/matrix.hpp index 80b5c5f8b..b6f7ef5cf 100644 --- a/include/ddc/kernels/splines/matrix.hpp +++ b/include/ddc/kernels/splines/matrix.hpp @@ -14,6 +14,11 @@ namespace ddc::detail { +/** + * @brief The parent class for matrices. + * + * It represents a Matrix storage format, methods to fill it and multiple-right-hand-size linear solvers. The matrix is necessarily squared. + */ class Matrix { int m_n; @@ -203,7 +208,7 @@ class Matrix virtual int factorize_method() = 0; /** - * @brief A function called by factorize() to actually perform the pre-process operation. + * @brief A function called by solve_inplace() and similar functions to actually perform the linear solve operation. * * @param b A double* to a contiguous array containing the (eventually multiple) right-hand-sides. Memory layout depends on the * derivated class (and the underlying algorithm). diff --git a/include/ddc/kernels/splines/matrix_sparse.hpp b/include/ddc/kernels/splines/matrix_sparse.hpp index b7f70ebcb..983aa277b 100644 --- a/include/ddc/kernels/splines/matrix_sparse.hpp +++ b/include/ddc/kernels/splines/matrix_sparse.hpp @@ -20,6 +20,8 @@ namespace ddc::detail { /** + * @brief Convert KokkosView to Ginkgo Dense matrix. + * * @param gko_exec[in] A Ginkgo executor that has access to the Kokkos::View memory space * @param view[in] A 2-D Kokkos::View with unit stride in the second dimension * @return A Ginkgo Dense matrix view over the Kokkos::View data @@ -41,6 +43,14 @@ auto to_gko_dense(std::shared_ptr const& gko_exec, KokkosVi view.stride_0()); } +/** + * @brief Return the default value of the hyperparameter cols_per_chunk for a given Kokkos::ExecutionSpace. + * + * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). They have been tuned on the basis of ddc/benchmarks/splines.cpp results. + * + * @tparam ExecSpace The Kokkos::ExecutionSpace type. + * @return The default value for the hyperparameter cols_per_chunk. + */ template int default_cols_per_chunk() noexcept { @@ -67,6 +77,14 @@ int default_cols_per_chunk() noexcept return 1; } +/** + * @brief Return the default value of the hyperparameter preconditionner_max_block_size for a given Kokkos::ExecutionSpace. + * + * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). They have been tuned on the basis of ddc/benchmarks/splines.cpp results. + * + * @tparam ExecSpace The Kokkos::ExecutionSpace type. + * @return The default value for the hyperparameter preconditionner_max_block_size. + */ template unsigned int default_preconditionner_max_block_size() noexcept { @@ -93,7 +111,13 @@ unsigned int default_preconditionner_max_block_size() noexcept return 1u; } -// Matrix class for sparse storage and iterative solve +/** + * @brief A sparse Matrix. + * + * The current storage format is Csr. Ginkgo is used to perform every matrix and linear solver-related operations. + * + * @tparam ExecSpace The Kokkos::ExecutionSpace on which operations related to the matrix are performed. + */ template class Matrix_Sparse : public Matrix { @@ -121,7 +145,13 @@ class Matrix_Sparse : public Matrix unsigned int m_preconditionner_max_block_size; // Maximum size of Jacobi-block preconditionner public: - // Constructor + /** + * @brief Matrix_Sparse constructor. + * + * @param mat_size The size of one of the dimensions of the matrix (which is necessarily squared). + * @param cols_per_chunk An optional hyperparameter used to define the number of right-hand-sides to pass to Ginkgo solver calls. see default_cols_per_chunk. + * @param preconditionner_max_block_size An optional hyperparameter used to define the maximum size of a block used by the block-Jacobi preconditionner. see default_preconditionner_max_block_size. + */ explicit Matrix_Sparse( const int mat_size, std::optional cols_per_chunk = std::nullopt, @@ -149,6 +179,15 @@ class Matrix_Sparse : public Matrix m_matrix_dense->at(i, j) = aij; } + /** + * @brief A function called by factorize() to perform the pre-process operation. + * + * Removes the zeros from the Csr object and instantiate a Ginkgo solver. It also construct a transposed version of the solver. + * + * The stopping criterion is a reduction factor ||x||/||b||<1e-15 with 1000 maximum iterations. + * + * @return The error code of the function. + */ int factorize_method() override { // Remove zeros @@ -185,6 +224,20 @@ class Matrix_Sparse : public Matrix return 0; } + /** + * @brief A function called by solve_inplace() and similar functions to actually perform the linear solve operation. + * + * The solver method is currently Bicgstab on CPU Serial and GPU and Gmres on OMP (because of Ginkgo issue #1563). + * + * Multiple-right-hand-sides are sliced in chunks of size cols_per_chunk which are passed one-after-the-other to Ginkgo. + * + * @param b A double* to a contiguous array containing the (eventually multiple) right-hand-sides. + * Because of Ginkgo limitation, it has to be right layout (row-major). + * @param transpose A character identifying if the normal ('N') or transposed ('T') linear system is solved. + * @param n_equations The number of multiple-right-hand-sides (number of columns in b). + * + * @return The error code of the function. + */ virtual int solve_inplace_method(double* b, char transpose, int n_equations) const override { std::shared_ptr const gko_exec = m_solver->get_executor(); From 5a64daabe3b213cd01e9e9cec52df40fea45bce3 Mon Sep 17 00:00:00 2001 From: blegouix Date: Thu, 2 May 2024 12:06:05 +0200 Subject: [PATCH 04/18] MatrixMaker --- include/ddc/kernels/splines/matrix_maker.hpp | 22 ++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/include/ddc/kernels/splines/matrix_maker.hpp b/include/ddc/kernels/splines/matrix_maker.hpp index 2064af749..f6c60f2d1 100644 --- a/include/ddc/kernels/splines/matrix_maker.hpp +++ b/include/ddc/kernels/splines/matrix_maker.hpp @@ -11,9 +11,31 @@ namespace ddc::detail { +/** + * @brief A static factory for Matrix classes. + * + * At the moment only the sparse matrix is available. + */ class MatrixMaker { public: + /** + * @brief Construct a sparse matrix + * + * @tparam the Kokkos::ExecutionSpace on which matrix-related operation will be performed. + * @param n The size of one of the dimensions of the matrix (which is necessarily squared). + * @param cols_per_chunk A hyperparameter used by the slicer (internal to the solver) to define the size + * of a chunk of right-hand-sides of the linear problem to be computed in parallel (chunks are treated + * by the linear solver one-after-the-other). + * + * This value is optional. If no value is provided then the default value is chosen by the requested solver. + * @param preconditionner_max_block_size A hyperparameter used by the slicer (internal to the solver) to + * define the size of a block used by the Block-Jacobi preconditioner. + * + * This value is optional. If no value is provided then the default value is chosen by the requested solver. + * + * @return The Matrix instance. + */ template static std::unique_ptr make_new_sparse( int const n, From 11486a0ee09108d5d4692cd17727b73dd36a9254 Mon Sep 17 00:00:00 2001 From: blegouix Date: Mon, 6 May 2024 10:58:59 +0200 Subject: [PATCH 05/18] Emily's review --- include/ddc/kernels/splines/matrix_maker.hpp | 3 +-- include/ddc/kernels/splines/matrix_sparse.hpp | 8 ++++---- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/include/ddc/kernels/splines/matrix_maker.hpp b/include/ddc/kernels/splines/matrix_maker.hpp index f6c60f2d1..6f3970c06 100644 --- a/include/ddc/kernels/splines/matrix_maker.hpp +++ b/include/ddc/kernels/splines/matrix_maker.hpp @@ -27,11 +27,10 @@ class MatrixMaker * @param cols_per_chunk A hyperparameter used by the slicer (internal to the solver) to define the size * of a chunk of right-hand-sides of the linear problem to be computed in parallel (chunks are treated * by the linear solver one-after-the-other). - * * This value is optional. If no value is provided then the default value is chosen by the requested solver. + * * @param preconditionner_max_block_size A hyperparameter used by the slicer (internal to the solver) to * define the size of a block used by the Block-Jacobi preconditioner. - * * This value is optional. If no value is provided then the default value is chosen by the requested solver. * * @return The Matrix instance. diff --git a/include/ddc/kernels/splines/matrix_sparse.hpp b/include/ddc/kernels/splines/matrix_sparse.hpp index 983aa277b..6f8815546 100644 --- a/include/ddc/kernels/splines/matrix_sparse.hpp +++ b/include/ddc/kernels/splines/matrix_sparse.hpp @@ -46,7 +46,7 @@ auto to_gko_dense(std::shared_ptr const& gko_exec, KokkosVi /** * @brief Return the default value of the hyperparameter cols_per_chunk for a given Kokkos::ExecutionSpace. * - * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). They have been tuned on the basis of ddc/benchmarks/splines.cpp results. + * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). They have been tuned on the basis of ddc/benchmarks/splines.cpp results on 4xIntel 6230 + Nvidia V100. * * @tparam ExecSpace The Kokkos::ExecutionSpace type. * @return The default value for the hyperparameter cols_per_chunk. @@ -80,7 +80,7 @@ int default_cols_per_chunk() noexcept /** * @brief Return the default value of the hyperparameter preconditionner_max_block_size for a given Kokkos::ExecutionSpace. * - * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). They have been tuned on the basis of ddc/benchmarks/splines.cpp results. + * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). They have been tuned on the basis of ddc/benchmarks/splines.cpp results on 4xIntel 6230 + Nvidia V100. * * @tparam ExecSpace The Kokkos::ExecutionSpace type. * @return The default value for the hyperparameter preconditionner_max_block_size. @@ -114,7 +114,7 @@ unsigned int default_preconditionner_max_block_size() noexcept /** * @brief A sparse Matrix. * - * The current storage format is Csr. Ginkgo is used to perform every matrix and linear solver-related operations. + * The current storage format is CSR. Ginkgo is used to perform every matrix and linear solver-related operations. * * @tparam ExecSpace The Kokkos::ExecutionSpace on which operations related to the matrix are performed. */ @@ -182,7 +182,7 @@ class Matrix_Sparse : public Matrix /** * @brief A function called by factorize() to perform the pre-process operation. * - * Removes the zeros from the Csr object and instantiate a Ginkgo solver. It also construct a transposed version of the solver. + * Removes the zeros from the CSR object and instantiate a Ginkgo solver. It also constructs a transposed version of the solver. * * The stopping criterion is a reduction factor ||x||/||b||<1e-15 with 1000 maximum iterations. * From 2dac78a004a8ac24ec52b3f7a89c7195a1bd5707 Mon Sep 17 00:00:00 2001 From: blegouix Date: Mon, 6 May 2024 11:17:05 +0200 Subject: [PATCH 06/18] minor --- include/ddc/kernels/splines/matrix.hpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/include/ddc/kernels/splines/matrix.hpp b/include/ddc/kernels/splines/matrix.hpp index b6f7ef5cf..44808d4a3 100644 --- a/include/ddc/kernels/splines/matrix.hpp +++ b/include/ddc/kernels/splines/matrix.hpp @@ -111,8 +111,8 @@ class Matrix /** * @brief Solve the multiple-right-hand-side linear problem Ax=b inplace. * - * @param[in, out] bx A 2D mdpsan storing the multiple right-hand-sides of the problem and receiving the corresponding solution. Important note: the convention is the reverse of the common matrix one, - * row number is second index and column number the first one. It should be changed in the future. + * @param[in, out] bx A 2D mdpsan storing the multiple right-hand-sides of the problem and receiving the corresponding solution. + * Important note: the convention is the reverse of the common matrix one, row number is second index and column number the first one. It should be changed in the future. * * @return bx */ @@ -131,8 +131,8 @@ class Matrix /** * @brief Solve the transposed multiple-right-hand-side linear problem A^tx=b inplace. * - * @param[in, out] bx A 2D mdpsan storing the multiple right-hand-sides of the problem and receiving the corresponding solution. Important note: the convention is the reverse of the common matrix one, - * row number is second index and column number the first one. It should be changed in the future. + * @param[in, out] bx A 2D mdpsan storing the multiple right-hand-sides of the problem and receiving the corresponding solution. + * Important note: the convention is the reverse of the common matrix one, row number is second index and column number the first one. It should be changed in the future. * * @return bx */ @@ -151,8 +151,8 @@ class Matrix /** * @brief Solve the multiple-right-hand-side linear problem Ax=b inplace. * - * @param[in, out] bx A 2D Kokkos::View storing the multiple right-hand-sides of the problem and receiving the corresponding solution. Important note: the convention is the reverse of the common matrix one, - * row number is second index and column number the first one. It should be changed in the future. + * @param[in, out] bx A 2D Kokkos::View storing the multiple right-hand-sides of the problem and receiving the corresponding solution. + * Important note: the convention is the reverse of the common matrix one, row number is second index and column number the first one. It should be changed in the future. * * @return bx */ From 65944860045d54b657b60d7e52f36af12d330426 Mon Sep 17 00:00:00 2001 From: blegouix Date: Mon, 6 May 2024 12:40:44 +0200 Subject: [PATCH 07/18] limit number characters per line --- include/ddc/kernels/splines/matrix.hpp | 9 ++++++--- include/ddc/kernels/splines/matrix_sparse.hpp | 12 ++++++++---- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/include/ddc/kernels/splines/matrix.hpp b/include/ddc/kernels/splines/matrix.hpp index 44808d4a3..4a511d4ff 100644 --- a/include/ddc/kernels/splines/matrix.hpp +++ b/include/ddc/kernels/splines/matrix.hpp @@ -112,7 +112,8 @@ class Matrix * @brief Solve the multiple-right-hand-side linear problem Ax=b inplace. * * @param[in, out] bx A 2D mdpsan storing the multiple right-hand-sides of the problem and receiving the corresponding solution. - * Important note: the convention is the reverse of the common matrix one, row number is second index and column number the first one. It should be changed in the future. + * Important note: the convention is the reverse of the common matrix one, row number is second index and column number + * the first one. It should be changed in the future. * * @return bx */ @@ -132,7 +133,8 @@ class Matrix * @brief Solve the transposed multiple-right-hand-side linear problem A^tx=b inplace. * * @param[in, out] bx A 2D mdpsan storing the multiple right-hand-sides of the problem and receiving the corresponding solution. - * Important note: the convention is the reverse of the common matrix one, row number is second index and column number the first one. It should be changed in the future. + * Important note: the convention is the reverse of the common matrix one, row number is second index and column number + * the first one. It should be changed in the future. * * @return bx */ @@ -152,7 +154,8 @@ class Matrix * @brief Solve the multiple-right-hand-side linear problem Ax=b inplace. * * @param[in, out] bx A 2D Kokkos::View storing the multiple right-hand-sides of the problem and receiving the corresponding solution. - * Important note: the convention is the reverse of the common matrix one, row number is second index and column number the first one. It should be changed in the future. + * Important note: the convention is the reverse of the common matrix one, row number is second index and column number + * the first one. It should be changed in the future. * * @return bx */ diff --git a/include/ddc/kernels/splines/matrix_sparse.hpp b/include/ddc/kernels/splines/matrix_sparse.hpp index 6f8815546..8f1d5eeb2 100644 --- a/include/ddc/kernels/splines/matrix_sparse.hpp +++ b/include/ddc/kernels/splines/matrix_sparse.hpp @@ -46,7 +46,8 @@ auto to_gko_dense(std::shared_ptr const& gko_exec, KokkosVi /** * @brief Return the default value of the hyperparameter cols_per_chunk for a given Kokkos::ExecutionSpace. * - * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). They have been tuned on the basis of ddc/benchmarks/splines.cpp results on 4xIntel 6230 + Nvidia V100. + * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). + * They have been tuned on the basis of ddc/benchmarks/splines.cpp results on 4xIntel 6230 + Nvidia V100. * * @tparam ExecSpace The Kokkos::ExecutionSpace type. * @return The default value for the hyperparameter cols_per_chunk. @@ -80,7 +81,8 @@ int default_cols_per_chunk() noexcept /** * @brief Return the default value of the hyperparameter preconditionner_max_block_size for a given Kokkos::ExecutionSpace. * - * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). They have been tuned on the basis of ddc/benchmarks/splines.cpp results on 4xIntel 6230 + Nvidia V100. + * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). + * They have been tuned on the basis of ddc/benchmarks/splines.cpp results on 4xIntel 6230 + Nvidia V100. * * @tparam ExecSpace The Kokkos::ExecutionSpace type. * @return The default value for the hyperparameter preconditionner_max_block_size. @@ -149,8 +151,10 @@ class Matrix_Sparse : public Matrix * @brief Matrix_Sparse constructor. * * @param mat_size The size of one of the dimensions of the matrix (which is necessarily squared). - * @param cols_per_chunk An optional hyperparameter used to define the number of right-hand-sides to pass to Ginkgo solver calls. see default_cols_per_chunk. - * @param preconditionner_max_block_size An optional hyperparameter used to define the maximum size of a block used by the block-Jacobi preconditionner. see default_preconditionner_max_block_size. + * @param cols_per_chunk An optional hyperparameter used to define the number of right-hand-sides to pass to + * Ginkgo solver calls. see default_cols_per_chunk. + * @param preconditionner_max_block_size An optional hyperparameter used to define the maximum size of a block + * used by the block-Jacobi preconditionner. see default_preconditionner_max_block_size. */ explicit Matrix_Sparse( const int mat_size, From a663982208ea197c47565b13055a5cdf929a51ac Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Tue, 14 May 2024 18:21:58 +0200 Subject: [PATCH 08/18] Apply suggestions from code review Co-authored-by: Thomas Padioleau --- include/ddc/kernels/splines/matrix.hpp | 33 +++++++++---------- include/ddc/kernels/splines/matrix_maker.hpp | 10 +++--- include/ddc/kernels/splines/matrix_sparse.hpp | 23 +++++++------ 3 files changed, 31 insertions(+), 35 deletions(-) diff --git a/include/ddc/kernels/splines/matrix.hpp b/include/ddc/kernels/splines/matrix.hpp index 4a511d4ff..13b6181ae 100644 --- a/include/ddc/kernels/splines/matrix.hpp +++ b/include/ddc/kernels/splines/matrix.hpp @@ -17,7 +17,7 @@ namespace ddc::detail { /** * @brief The parent class for matrices. * - * It represents a Matrix storage format, methods to fill it and multiple-right-hand-size linear solvers. The matrix is necessarily squared. + * It represents a square Matrix. Implementations may have different storage formats, filling methods and multiple right-hand sides linear solvers. */ class Matrix { @@ -30,17 +30,17 @@ class Matrix virtual ~Matrix() = default; /** - * @brief Get an element of the matrix at indexes i, j. + * @brief Get an element of the matrix at indexes i, j. It must not be called after `factorize`. * * @param i The row index of the desired element. - * @param j The columns index of the desired element. + * @param j The column index of the desired element. * * @return The value of the element of the matrix. */ virtual double get_element(int i, int j) const = 0; /** - * @brief Set an element of the matrix at indexes i, j. + * @brief Set an element of the matrix at indexes i, j. It must not be called after `factorize`. * * @param i The row index of the setted element. * @param j The column index of the setted element. @@ -75,7 +75,7 @@ class Matrix /** * @brief Solve the linear problem Ax=b inplace. * - * @param[in, out] b A 1D mdpsan storing a right-hand-side of the problem and receiving the corresponding solution. + * @param[in, out] b A 1D mdpsan storing a right-hand side of the problem and receiving the corresponding solution. */ virtual ddc::DSpan1D solve_inplace(ddc::DSpan1D const b) const { @@ -92,7 +92,7 @@ class Matrix /** * @brief Solve the transposed linear problem A^tx=b inplace. * - * @param[in, out] b A 1D mdpsan storing a right-hand-side of the problem and receiving the corresponding solution. + * @param[in, out] b A 1D mdpsan storing a right-hand side of the problem and receiving the corresponding solution. * * @return b */ @@ -109,11 +109,11 @@ class Matrix } /** - * @brief Solve the multiple-right-hand-side linear problem Ax=b inplace. + * @brief Solve the multiple right-hand sides linear problem Ax=b inplace. * - * @param[in, out] bx A 2D mdpsan storing the multiple right-hand-sides of the problem and receiving the corresponding solution. + * @param[in, out] bx A 2D mdpsan storing the multiple right-hand sides of the problem and receiving the corresponding solution. * Important note: the convention is the reverse of the common matrix one, row number is second index and column number - * the first one. It should be changed in the future. + * the first one. This means when solving `A x_i = b_i`, element `(b_i)_j` is stored in `b(j, i)`. * * @return bx */ @@ -130,9 +130,9 @@ class Matrix } /** - * @brief Solve the transposed multiple-right-hand-side linear problem A^tx=b inplace. + * @brief Solve the transposed multiple right-hand sides linear problem A^tx=b inplace. * - * @param[in, out] bx A 2D mdpsan storing the multiple right-hand-sides of the problem and receiving the corresponding solution. + * @param[in, out] bx A 2D mdspan storing the multiple right-hand sides of the problem and receiving the corresponding solution. * Important note: the convention is the reverse of the common matrix one, row number is second index and column number * the first one. It should be changed in the future. * @@ -151,9 +151,9 @@ class Matrix } /** - * @brief Solve the multiple-right-hand-side linear problem Ax=b inplace. + * @brief Solve the multiple-right-hand sides linear problem Ax=b inplace. * - * @param[in, out] bx A 2D Kokkos::View storing the multiple right-hand-sides of the problem and receiving the corresponding solution. + * @param[in, out] bx A 2D Kokkos::View storing the multiple right-hand sides of the problem and receiving the corresponding solution. * Important note: the convention is the reverse of the common matrix one, row number is second index and column number * the first one. It should be changed in the future. * @@ -174,7 +174,7 @@ class Matrix } /** - * @brief Get the size of the matrix (which is necessarily squared) in one of its dimensions. + * @brief Get the size of the square matrix in one of its dimensions. * * @return The size of the matrix in one of its dimensions. */ @@ -184,7 +184,7 @@ class Matrix } /** - * @brief Prints a Matrix in a std::ostream. It will segfault is the Matrix is on GPU. + * @brief Prints a Matrix in a std::ostream. It must not be called after `factorize`. * * @param out The stream in which the matrix is printed. * @@ -213,8 +213,7 @@ class Matrix /** * @brief A function called by solve_inplace() and similar functions to actually perform the linear solve operation. * - * @param b A double* to a contiguous array containing the (eventually multiple) right-hand-sides. Memory layout depends on the - * derivated class (and the underlying algorithm). + * @param b A double* to a contiguous array containing the (eventually multiple) right-hand-sides. The memory layout is right. * @param transpose A character identifying if the normal ('N') or transposed ('T') linear system is solved. * @param n_equations The number of multiple-right-hand-sides (number of columns of b). * @return The error code of the function. diff --git a/include/ddc/kernels/splines/matrix_maker.hpp b/include/ddc/kernels/splines/matrix_maker.hpp index 6f3970c06..189649960 100644 --- a/include/ddc/kernels/splines/matrix_maker.hpp +++ b/include/ddc/kernels/splines/matrix_maker.hpp @@ -13,8 +13,6 @@ namespace ddc::detail { /** * @brief A static factory for Matrix classes. - * - * At the moment only the sparse matrix is available. */ class MatrixMaker { @@ -23,13 +21,13 @@ class MatrixMaker * @brief Construct a sparse matrix * * @tparam the Kokkos::ExecutionSpace on which matrix-related operation will be performed. - * @param n The size of one of the dimensions of the matrix (which is necessarily squared). - * @param cols_per_chunk A hyperparameter used by the slicer (internal to the solver) to define the size - * of a chunk of right-hand-sides of the linear problem to be computed in parallel (chunks are treated + * @param n The size of one of the dimensions of the square matrix. + * @param cols_per_chunk A parameter used by the slicer (internal to the solver) to define the size + * of a chunk of right-hand sides of the linear problem to be computed in parallel (chunks are treated * by the linear solver one-after-the-other). * This value is optional. If no value is provided then the default value is chosen by the requested solver. * - * @param preconditionner_max_block_size A hyperparameter used by the slicer (internal to the solver) to + * @param preconditionner_max_block_size A parameter used by the slicer (internal to the solver) to * define the size of a block used by the Block-Jacobi preconditioner. * This value is optional. If no value is provided then the default value is chosen by the requested solver. * diff --git a/include/ddc/kernels/splines/matrix_sparse.hpp b/include/ddc/kernels/splines/matrix_sparse.hpp index 8f1d5eeb2..fb3f08950 100644 --- a/include/ddc/kernels/splines/matrix_sparse.hpp +++ b/include/ddc/kernels/splines/matrix_sparse.hpp @@ -44,13 +44,13 @@ auto to_gko_dense(std::shared_ptr const& gko_exec, KokkosVi } /** - * @brief Return the default value of the hyperparameter cols_per_chunk for a given Kokkos::ExecutionSpace. + * @brief Return the default value of the parameter cols_per_chunk for a given Kokkos::ExecutionSpace. * * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). * They have been tuned on the basis of ddc/benchmarks/splines.cpp results on 4xIntel 6230 + Nvidia V100. * * @tparam ExecSpace The Kokkos::ExecutionSpace type. - * @return The default value for the hyperparameter cols_per_chunk. + * @return The default value for the parameter cols_per_chunk. */ template int default_cols_per_chunk() noexcept @@ -79,13 +79,13 @@ int default_cols_per_chunk() noexcept } /** - * @brief Return the default value of the hyperparameter preconditionner_max_block_size for a given Kokkos::ExecutionSpace. + * @brief Return the default value of the parameter preconditionner_max_block_size for a given Kokkos::ExecutionSpace. * * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). * They have been tuned on the basis of ddc/benchmarks/splines.cpp results on 4xIntel 6230 + Nvidia V100. * * @tparam ExecSpace The Kokkos::ExecutionSpace type. - * @return The default value for the hyperparameter preconditionner_max_block_size. + * @return The default value for the parameter preconditionner_max_block_size. */ template unsigned int default_preconditionner_max_block_size() noexcept @@ -116,7 +116,7 @@ unsigned int default_preconditionner_max_block_size() noexcept /** * @brief A sparse Matrix. * - * The current storage format is CSR. Ginkgo is used to perform every matrix and linear solver-related operations. + * The storage format is CSR. Ginkgo is used to perform every matrix and linear solver-related operations. * * @tparam ExecSpace The Kokkos::ExecutionSpace on which operations related to the matrix are performed. */ @@ -150,10 +150,10 @@ class Matrix_Sparse : public Matrix /** * @brief Matrix_Sparse constructor. * - * @param mat_size The size of one of the dimensions of the matrix (which is necessarily squared). - * @param cols_per_chunk An optional hyperparameter used to define the number of right-hand-sides to pass to + * @param mat_size The size of one of the dimensions of the square matrix. + * @param cols_per_chunk An optional parameter used to define the number of right-hand sides to pass to * Ginkgo solver calls. see default_cols_per_chunk. - * @param preconditionner_max_block_size An optional hyperparameter used to define the maximum size of a block + * @param preconditionner_max_block_size An optional parameter used to define the maximum size of a block * used by the block-Jacobi preconditionner. see default_preconditionner_max_block_size. */ explicit Matrix_Sparse( @@ -233,12 +233,11 @@ class Matrix_Sparse : public Matrix * * The solver method is currently Bicgstab on CPU Serial and GPU and Gmres on OMP (because of Ginkgo issue #1563). * - * Multiple-right-hand-sides are sliced in chunks of size cols_per_chunk which are passed one-after-the-other to Ginkgo. + * Multiple right-hand sides are sliced in chunks of size cols_per_chunk which are passed one-after-the-other to Ginkgo. * - * @param b A double* to a contiguous array containing the (eventually multiple) right-hand-sides. - * Because of Ginkgo limitation, it has to be right layout (row-major). + * @param b A double* to a contiguous array containing the (eventually multiple) right-hand sides. * @param transpose A character identifying if the normal ('N') or transposed ('T') linear system is solved. - * @param n_equations The number of multiple-right-hand-sides (number of columns in b). + * @param n_equations The number of multiple right-hand sides (number of columns in b). * * @return The error code of the function. */ From dfee1122b6ca68dd8afcf39f32e3ae832276ce69 Mon Sep 17 00:00:00 2001 From: blegouix Date: Wed, 15 May 2024 09:44:24 +0200 Subject: [PATCH 09/18] wip --- include/ddc/kernels/splines/matrix.hpp | 84 +++++--------------------- 1 file changed, 16 insertions(+), 68 deletions(-) diff --git a/include/ddc/kernels/splines/matrix.hpp b/include/ddc/kernels/splines/matrix.hpp index 13b6181ae..7d34b3eb7 100644 --- a/include/ddc/kernels/splines/matrix.hpp +++ b/include/ddc/kernels/splines/matrix.hpp @@ -15,19 +15,26 @@ namespace ddc::detail { /** - * @brief The parent class for matrices. + * @brief The parent class for a linear problem dedicated to compute a spline approximation. * - * It represents a square Matrix. Implementations may have different storage formats, filling methods and multiple right-hand sides linear solvers. + * It represents a square Matrix and provides methods to solve a multiple right-hand sides linear problem. + * Implementations may have different storage formats, filling methods and multiple right-hand sides linear solvers. */ -class Matrix +class SplinesLinearProblem { - int m_n; - public: - explicit Matrix(const int mat_size) : m_n(mat_size) {} + /// @brief The type of a Kokkos::View storing multiple right-hand sides. + using MultiRHS = Kokkos::View; + +private: + std::size_t m_size; + +protected: + explicit SplinesLinearProblem(const std::size_t size) : m_size(size) {} +public: /// @brief Destruct - virtual ~Matrix() = default; + virtual ~SplinesLinearProblem() = default; /** * @brief Get an element of the matrix at indexes i, j. It must not be called after `factorize`. @@ -54,7 +61,7 @@ class Matrix * Note: this function should be renamed in the future because the pre- * process operation is not necessarily a factorization. */ - virtual void factorize() + virtual void finished_filling() { int const info = factorize_method(); @@ -72,42 +79,6 @@ class Matrix } } - /** - * @brief Solve the linear problem Ax=b inplace. - * - * @param[in, out] b A 1D mdpsan storing a right-hand side of the problem and receiving the corresponding solution. - */ - virtual ddc::DSpan1D solve_inplace(ddc::DSpan1D const b) const - { - assert(int(b.extent(0)) == m_n); - int const info = solve_inplace_method(b.data_handle(), 'N', 1); - - if (info < 0) { - std::cerr << -info << "-th argument had an illegal value" << std::endl; - // TODO: Add LOG_FATAL_ERROR - } - return b; - } - - /** - * @brief Solve the transposed linear problem A^tx=b inplace. - * - * @param[in, out] b A 1D mdpsan storing a right-hand side of the problem and receiving the corresponding solution. - * - * @return b - */ - virtual ddc::DSpan1D solve_transpose_inplace(ddc::DSpan1D const b) const - { - assert(int(b.extent(0)) == m_n); - int const info = solve_inplace_method(b.data_handle(), 'T', 1); - - if (info < 0) { - std::cerr << -info << "-th argument had an illegal value" << std::endl; - // TODO: Add LOG_FATAL_ERROR - } - return b; - } - /** * @brief Solve the multiple right-hand sides linear problem Ax=b inplace. * @@ -150,35 +121,12 @@ class Matrix return bx; } - /** - * @brief Solve the multiple-right-hand sides linear problem Ax=b inplace. - * - * @param[in, out] bx A 2D Kokkos::View storing the multiple right-hand sides of the problem and receiving the corresponding solution. - * Important note: the convention is the reverse of the common matrix one, row number is second index and column number - * the first one. It should be changed in the future. - * - * @return bx - */ - template - Kokkos::View solve_batch_inplace( - Kokkos::View const bx) const - { - assert(int(bx.extent(0)) == m_n); - int const info = solve_inplace_method(bx.data(), 'N', bx.extent(1)); - - if (info < 0) { - std::cerr << -info << "-th argument had an illegal value" << std::endl; - // TODO: Add LOG_FATAL_ERROR - } - return bx; - } - /** * @brief Get the size of the square matrix in one of its dimensions. * * @return The size of the matrix in one of its dimensions. */ - int get_size() const + int size() const { return m_n; } From 37e8f40aa299f43467239613c04bfa1780cbd714 Mon Sep 17 00:00:00 2001 From: blegouix Date: Wed, 15 May 2024 12:51:10 +0200 Subject: [PATCH 10/18] init --- include/ddc/kernels/splines/matrix.hpp | 91 +++---------------- include/ddc/kernels/splines/matrix_maker.hpp | 12 +-- include/ddc/kernels/splines/matrix_sparse.hpp | 76 ++++++++-------- .../ddc/kernels/splines/spline_builder.hpp | 10 +- tests/splines/matrix.cpp | 41 ++++++--- 5 files changed, 86 insertions(+), 144 deletions(-) diff --git a/include/ddc/kernels/splines/matrix.hpp b/include/ddc/kernels/splines/matrix.hpp index 7d34b3eb7..7bf108885 100644 --- a/include/ddc/kernels/splines/matrix.hpp +++ b/include/ddc/kernels/splines/matrix.hpp @@ -15,16 +15,17 @@ namespace ddc::detail { /** - * @brief The parent class for a linear problem dedicated to compute a spline approximation. + * @brief The parent class for a linear problem dedicated to the computation of a spline approximation. * * It represents a square Matrix and provides methods to solve a multiple right-hand sides linear problem. * Implementations may have different storage formats, filling methods and multiple right-hand sides linear solvers. */ +template class SplinesLinearProblem { public: /// @brief The type of a Kokkos::View storing multiple right-hand sides. - using MultiRHS = Kokkos::View; + using MultiRHS = Kokkos::View; private: std::size_t m_size; @@ -56,70 +57,18 @@ class SplinesLinearProblem virtual void set_element(int i, int j, double aij) = 0; /** - * @brief Performs a pre-process operation on the Matrix. - * - * Note: this function should be renamed in the future because the pre- - * process operation is not necessarily a factorization. + * @brief Performs a pre-process operation on the solver. */ - virtual void finished_filling() - { - int const info = factorize_method(); - - if (info < 0) { - std::cerr << -info << "-th argument had an illegal value" << std::endl; - // TODO: Add LOG_FATAL_ERROR - } else if (info > 0) { - std::cerr << "U(" << info << "," << info << ") is exactly zero."; - std::cerr << " The factorization has been completed, but the factor"; - std::cerr << " U is exactly singular, and division by zero will occur " - "if " - "it is used to"; - std::cerr << " solve a system of equations."; - // TODO: Add LOG_FATAL_ERROR - } - } + virtual void setup_solver() = 0; /** * @brief Solve the multiple right-hand sides linear problem Ax=b inplace. * - * @param[in, out] bx A 2D mdpsan storing the multiple right-hand sides of the problem and receiving the corresponding solution. - * Important note: the convention is the reverse of the common matrix one, row number is second index and column number - * the first one. This means when solving `A x_i = b_i`, element `(b_i)_j` is stored in `b(j, i)`. - * - * @return bx - */ - virtual ddc::DSpan2D solve_multiple_inplace(ddc::DSpan2D const bx) const - { - assert(int(bx.extent(1)) == m_n); - int const info = solve_inplace_method(bx.data_handle(), 'N', bx.extent(0)); - - if (info < 0) { - std::cerr << -info << "-th argument had an illegal value" << std::endl; - // TODO: Add LOG_FATAL_ERROR - } - return bx; - } - - /** - * @brief Solve the transposed multiple right-hand sides linear problem A^tx=b inplace. - * - * @param[in, out] bx A 2D mdspan storing the multiple right-hand sides of the problem and receiving the corresponding solution. + * @param[in, out] multi_rhs A 2D Kokkos::View storing the multiple right-hand sides of the problem and receiving the corresponding solution. * Important note: the convention is the reverse of the common matrix one, row number is second index and column number - * the first one. It should be changed in the future. - * - * @return bx + * the first one. This means when solving `A x_i = b_i`, element `b_{ij}` is stored in `b(j, i)`. */ - virtual ddc::DSpan2D solve_multiple_transpose_inplace(ddc::DSpan2D const bx) const - { - assert(int(bx.extent(1)) == m_n); - int const info = solve_inplace_method(bx.data_handle(), 'T', bx.extent(0)); - - if (info < 0) { - std::cerr << -info << "-th argument had an illegal value" << std::endl; - // TODO: Add LOG_FATAL_ERROR - } - return bx; - } + virtual void solve(MultiRHS b, bool transpose = false) const = 0; /** * @brief Get the size of the square matrix in one of its dimensions. @@ -128,11 +77,11 @@ class SplinesLinearProblem */ int size() const { - return m_n; + return m_size; } /** - * @brief Prints a Matrix in a std::ostream. It must not be called after `factorize`. + * @brief Prints a Matrix in a std::ostream. It must not be called after `setup_solver`. * * @param out The stream in which the matrix is printed. * @@ -140,7 +89,7 @@ class SplinesLinearProblem **/ std::ostream& operator<<(std::ostream& os) { - int const n = get_size(); + int const n = size(); for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { os << std::fixed << std::setprecision(3) << std::setw(10) << get_element(i, j); @@ -149,24 +98,6 @@ class SplinesLinearProblem } return os; } - -protected: - /** - * @brief A function called by factorize() to actually perform the pre-process operation. - * - * @return The error code of the function. - */ - virtual int factorize_method() = 0; - - /** - * @brief A function called by solve_inplace() and similar functions to actually perform the linear solve operation. - * - * @param b A double* to a contiguous array containing the (eventually multiple) right-hand-sides. The memory layout is right. - * @param transpose A character identifying if the normal ('N') or transposed ('T') linear system is solved. - * @param n_equations The number of multiple-right-hand-sides (number of columns of b). - * @return The error code of the function. - */ - virtual int solve_inplace_method(double* b, char transpose, int n_equations) const = 0; }; } // namespace ddc::detail diff --git a/include/ddc/kernels/splines/matrix_maker.hpp b/include/ddc/kernels/splines/matrix_maker.hpp index 189649960..f9dbf7373 100644 --- a/include/ddc/kernels/splines/matrix_maker.hpp +++ b/include/ddc/kernels/splines/matrix_maker.hpp @@ -12,9 +12,9 @@ namespace ddc::detail { /** - * @brief A static factory for Matrix classes. + * @brief A static factory for SplinesLinearProblem classes. */ -class MatrixMaker +class SplinesLinearProblemMaker { public: /** @@ -31,16 +31,16 @@ class MatrixMaker * define the size of a block used by the Block-Jacobi preconditioner. * This value is optional. If no value is provided then the default value is chosen by the requested solver. * - * @return The Matrix instance. + * @return The SplinesLinearProblem instance. */ template - static std::unique_ptr make_new_sparse( + static std::unique_ptr> make_new_sparse( int const n, std::optional cols_per_chunk = std::nullopt, std::optional preconditionner_max_block_size = std::nullopt) { - return std::make_unique< - Matrix_Sparse>(n, cols_per_chunk, preconditionner_max_block_size); + return std::make_unique>(n, cols_per_chunk, preconditionner_max_block_size); } }; diff --git a/include/ddc/kernels/splines/matrix_sparse.hpp b/include/ddc/kernels/splines/matrix_sparse.hpp index fb3f08950..82b49069e 100644 --- a/include/ddc/kernels/splines/matrix_sparse.hpp +++ b/include/ddc/kernels/splines/matrix_sparse.hpp @@ -46,7 +46,7 @@ auto to_gko_dense(std::shared_ptr const& gko_exec, KokkosVi /** * @brief Return the default value of the parameter cols_per_chunk for a given Kokkos::ExecutionSpace. * - * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). + * The values are hardware-specific (but they can be overriden in the constructor of SplinesLinearProblemSparse). * They have been tuned on the basis of ddc/benchmarks/splines.cpp results on 4xIntel 6230 + Nvidia V100. * * @tparam ExecSpace The Kokkos::ExecutionSpace type. @@ -81,7 +81,7 @@ int default_cols_per_chunk() noexcept /** * @brief Return the default value of the parameter preconditionner_max_block_size for a given Kokkos::ExecutionSpace. * - * The values are hardware-specific (but they can be overriden in the constructor of MatrixSparse). + * The values are hardware-specific (but they can be overriden in the constructor of SplinesLinearProblemSparse). * They have been tuned on the basis of ddc/benchmarks/splines.cpp results on 4xIntel 6230 + Nvidia V100. * * @tparam ExecSpace The Kokkos::ExecutionSpace type. @@ -114,15 +114,21 @@ unsigned int default_preconditionner_max_block_size() noexcept } /** - * @brief A sparse Matrix. + * @brief A sparse linear problem dedicated to the computation of a spline approximation. * * The storage format is CSR. Ginkgo is used to perform every matrix and linear solver-related operations. * * @tparam ExecSpace The Kokkos::ExecutionSpace on which operations related to the matrix are performed. */ template -class Matrix_Sparse : public Matrix +class SplinesLinearProblemSparse : public SplinesLinearProblem { +public: + using SplinesLinearProblem::MultiRHS; + using SplinesLinearProblem::size; + using SplinesLinearProblem::operator<<; + +private: using matrix_sparse_type = gko::matrix::Csr; #ifdef KOKKOS_ENABLE_OPENMP using solver_type = std::conditional_t< @@ -148,7 +154,7 @@ class Matrix_Sparse : public Matrix public: /** - * @brief Matrix_Sparse constructor. + * @brief SplinesLinearProblemSparse constructor. * * @param mat_size The size of one of the dimensions of the square matrix. * @param cols_per_chunk An optional parameter used to define the number of right-hand sides to pass to @@ -156,11 +162,11 @@ class Matrix_Sparse : public Matrix * @param preconditionner_max_block_size An optional parameter used to define the maximum size of a block * used by the block-Jacobi preconditionner. see default_preconditionner_max_block_size. */ - explicit Matrix_Sparse( + explicit SplinesLinearProblemSparse( const int mat_size, std::optional cols_per_chunk = std::nullopt, std::optional preconditionner_max_block_size = std::nullopt) - : Matrix(mat_size) + : SplinesLinearProblem(mat_size) , m_cols_per_chunk(cols_per_chunk.value_or(default_cols_per_chunk())) , m_preconditionner_max_block_size(preconditionner_max_block_size.value_or( default_preconditionner_max_block_size())) @@ -174,8 +180,9 @@ class Matrix_Sparse : public Matrix virtual double get_element([[maybe_unused]] int i, [[maybe_unused]] int j) const override { - throw std::runtime_error("MatrixSparse::get_element() is not implemented because no API is " - "provided by Ginkgo"); + throw std::runtime_error( + "SplinesLinearProblemSparse::get_element() is not implemented because no API is " + "provided by Ginkgo"); } virtual void set_element(int i, int j, double aij) override @@ -192,10 +199,10 @@ class Matrix_Sparse : public Matrix * * @return The error code of the function. */ - int factorize_method() override + void setup_solver() override { // Remove zeros - gko::matrix_data matrix_data(gko::dim<2>(get_size(), get_size())); + gko::matrix_data matrix_data(gko::dim<2>(size(), size())); m_matrix_dense->write(matrix_data); m_matrix_dense.reset(); matrix_data.remove_zeros(); @@ -224,8 +231,6 @@ class Matrix_Sparse : public Matrix m_solver = solver_factory->generate(m_matrix_sparse); m_solver_tr = m_solver->transpose(); gko_exec->synchronize(); - - return 0; } /** @@ -238,20 +243,19 @@ class Matrix_Sparse : public Matrix * @param b A double* to a contiguous array containing the (eventually multiple) right-hand sides. * @param transpose A character identifying if the normal ('N') or transposed ('T') linear system is solved. * @param n_equations The number of multiple right-hand sides (number of columns in b). - * - * @return The error code of the function. */ - virtual int solve_inplace_method(double* b, char transpose, int n_equations) const override + void solve(typename SplinesLinearProblem::MultiRHS b, bool const transpose) + const override { + std::size_t const n_equations = b.extent(1); + assert(int(b.extent(0)) == size()); + std::shared_ptr const gko_exec = m_solver->get_executor(); std::shared_ptr const convergence_logger = gko::log::Convergence::create(); - int const main_chunk_size = std::min(m_cols_per_chunk, n_equations); + int const main_chunk_size = std::min(m_cols_per_chunk, (int)n_equations); - Kokkos::View const - b_view(b, get_size(), n_equations); - Kokkos::View const - x_view("", get_size(), main_chunk_size); + Kokkos::View const x("", size(), main_chunk_size); int const iend = (n_equations + main_chunk_size - 1) / main_chunk_size; for (int i = 0; i < iend; ++i) { @@ -259,37 +263,31 @@ class Matrix_Sparse : public Matrix int const subview_end = (i + 1 == iend) ? n_equations : (subview_begin + main_chunk_size); - auto const b_subview = Kokkos:: - subview(b_view, Kokkos::ALL, Kokkos::pair(subview_begin, subview_end)); - auto const x_subview = Kokkos:: - subview(x_view, Kokkos::ALL, Kokkos::pair(0, subview_end - subview_begin)); + auto const b_chunk + = Kokkos::subview(b, Kokkos::ALL, Kokkos::pair(subview_begin, subview_end)); + auto const x_chunk + = Kokkos::subview(x, Kokkos::ALL, Kokkos::pair(0, subview_end - subview_begin)); - Kokkos::deep_copy(x_subview, b_subview); + Kokkos::deep_copy(x_chunk, b_chunk); - if (transpose == 'N') { + if (!transpose) { m_solver->add_logger(convergence_logger); - m_solver - ->apply(to_gko_dense(gko_exec, b_subview), - to_gko_dense(gko_exec, x_subview)); + m_solver->apply(to_gko_dense(gko_exec, b_chunk), to_gko_dense(gko_exec, x_chunk)); m_solver->remove_logger(convergence_logger); - } else if (transpose == 'T') { + } else { m_solver_tr->add_logger(convergence_logger); m_solver_tr - ->apply(to_gko_dense(gko_exec, b_subview), - to_gko_dense(gko_exec, x_subview)); + ->apply(to_gko_dense(gko_exec, b_chunk), to_gko_dense(gko_exec, x_chunk)); m_solver_tr->remove_logger(convergence_logger); - } else { - throw std::domain_error("transpose option not recognized"); } if (!convergence_logger->has_converged()) { - throw std::runtime_error("Ginkgo did not converged in ddc::detail::Matrix_Sparse"); + throw std::runtime_error( + "Ginkgo did not converged in ddc::detail::SplinesLinearProblemSparse"); } - Kokkos::deep_copy(b_subview, x_subview); + Kokkos::deep_copy(b_chunk, x_chunk); } - - return 1; } }; diff --git a/include/ddc/kernels/splines/spline_builder.hpp b/include/ddc/kernels/splines/spline_builder.hpp index 173fcc455..84434ab21 100644 --- a/include/ddc/kernels/splines/spline_builder.hpp +++ b/include/ddc/kernels/splines/spline_builder.hpp @@ -136,7 +136,7 @@ class SplineBuilder double m_dx; // average cell size for normalization of derivatives // interpolator specific - std::unique_ptr matrix; + std::unique_ptr> matrix; /// Calculate offset so that the matrix is diagonally dominant int compute_offset(interpolation_domain_type const& interpolation_domain); @@ -264,7 +264,7 @@ class SplineBuilder * * @return A reference to the interpolation matrix. */ - const ddc::detail::Matrix& get_interpolation_matrix() const noexcept + const ddc::detail::SplinesLinearProblem& get_interpolation_matrix() const noexcept { return *matrix; } @@ -487,14 +487,14 @@ void SplineBuilder< return; */ - matrix = ddc::detail::MatrixMaker::make_new_sparse( + matrix = ddc::detail::SplinesLinearProblemMaker::make_new_sparse( ddc::discrete_space().nbasis(), cols_per_chunk, preconditionner_max_block_size); build_matrix_system(); - matrix->factorize(); + matrix->setup_solver(); } template < @@ -724,7 +724,7 @@ operator()( ddc::discrete_space().nbasis(), batch_domain().size()); // Compute spline coef - matrix->solve_batch_inplace(bcoef_section); + matrix->solve(bcoef_section); // Transpose back spline_tr in spline ddc::parallel_for_each( exec_space(), diff --git a/tests/splines/matrix.cpp b/tests/splines/matrix.cpp index 440ff6dc8..39c3e20f4 100644 --- a/tests/splines/matrix.cpp +++ b/tests/splines/matrix.cpp @@ -19,7 +19,8 @@ namespace { -void fill_identity(ddc::DSpan2D mat) +void fill_identity( + ddc::detail::SplinesLinearProblem::MultiRHS mat) { assert(mat.extent(0) == mat.extent(1)); for (std::size_t i(0); i < mat.extent(0); ++i) { @@ -30,7 +31,7 @@ void fill_identity(ddc::DSpan2D mat) } /* -void copy_matrix(ddc::DSpan2D copy, std::unique_ptr& mat) +void copy_matrix(ddc::detail::SplinesLinearProblem::MultiRHS copy, std::unique_ptr>& mat) { assert(mat->get_size() == int(copy.extent(0))); assert(mat->get_size() == int(copy.extent(1))); @@ -43,7 +44,9 @@ void copy_matrix(ddc::DSpan2D copy, std::unique_ptr& mat) } */ -void check_inverse(ddc::DSpan2D matrix, ddc::DSpan2D inv) +void check_inverse( + ddc::detail::SplinesLinearProblem::MultiRHS matrix, + ddc::detail::SplinesLinearProblem::MultiRHS inv) { double TOL = 1e-10; std::size_t N = matrix.extent(0); @@ -52,14 +55,16 @@ void check_inverse(ddc::DSpan2D matrix, ddc::DSpan2D inv) for (std::size_t j(0); j < N; ++j) { double id_val = 0.0; for (std::size_t k(0); k < N; ++k) { - id_val += matrix(i, k) * inv(j, k); + id_val += matrix(i, k) * inv(k, j); } EXPECT_NEAR(id_val, static_cast(i == j), TOL); } } } -void check_inverse_transpose(ddc::DSpan2D matrix, ddc::DSpan2D inv) +void check_inverse_transpose( + ddc::detail::SplinesLinearProblem::MultiRHS matrix, + ddc::detail::SplinesLinearProblem::MultiRHS inv) { double TOL = 1e-10; std::size_t N = matrix.extent(0); @@ -68,7 +73,7 @@ void check_inverse_transpose(ddc::DSpan2D matrix, ddc::DSpan2D inv) for (std::size_t j(0); j < N; ++j) { double id_val = 0.0; for (std::size_t k(0); k < N; ++k) { - id_val += matrix(i, k) * inv(k, j); + id_val += matrix(i, k) * inv(j, k); } EXPECT_NEAR(id_val, static_cast(i == j), TOL); } @@ -85,11 +90,13 @@ TEST_P(MatrixSizesFixture, Sparse) auto const [N, k] = GetParam(); - std::unique_ptr matrix - = ddc::detail::MatrixMaker::make_new_sparse(N); + std::unique_ptr> matrix + = ddc::detail::SplinesLinearProblemMaker::make_new_sparse< + Kokkos::DefaultExecutionSpace>(N); std::vector val_ptr(N * N); - ddc::DSpan2D val(val_ptr.data(), N, N); + ddc::detail::SplinesLinearProblem::MultiRHS + val(val_ptr.data(), N, N); for (std::size_t i(0); i < N; ++i) { for (std::size_t j(0); j < N; ++j) { if (i == j) { @@ -105,23 +112,29 @@ TEST_P(MatrixSizesFixture, Sparse) } // copy_matrix(val, matrix); // copy_matrix is not available for sparse matrix because of a limitation of Ginkgo API (get_element is not implemented). The workaround is to fill val directly in the loop - matrix->factorize(); + matrix->setup_solver(); Kokkos::DualView inv_ptr("inv_ptr", N * N); - ddc::DSpan2D inv(inv_ptr.h_view.data(), N, N); + ddc::detail::SplinesLinearProblem::MultiRHS + inv(inv_ptr.h_view.data(), N, N); fill_identity(inv); inv_ptr.modify_host(); inv_ptr.sync_device(); - matrix->solve_multiple_inplace(ddc::DSpan2D(inv_ptr.d_view.data(), N, N)); + matrix->solve(ddc::detail::SplinesLinearProblem< + Kokkos::DefaultExecutionSpace>::MultiRHS(inv_ptr.d_view.data(), N, N)); inv_ptr.modify_device(); inv_ptr.sync_host(); Kokkos::DualView inv_tr_ptr("inv_tr_ptr", N * N); - ddc::DSpan2D inv_tr(inv_tr_ptr.h_view.data(), N, N); + ddc::detail::SplinesLinearProblem::MultiRHS + inv_tr(inv_tr_ptr.h_view.data(), N, N); fill_identity(inv_tr); inv_tr_ptr.modify_host(); inv_tr_ptr.sync_device(); - matrix->solve_multiple_transpose_inplace(ddc::DSpan2D(inv_tr_ptr.d_view.data(), N, N)); + matrix + ->solve(ddc::detail::SplinesLinearProblem:: + MultiRHS(inv_tr_ptr.d_view.data(), N, N), + true); inv_tr_ptr.modify_device(); inv_tr_ptr.sync_host(); From aba480de3d2396f9f225d01ab7819a341e99c475 Mon Sep 17 00:00:00 2001 From: blegouix Date: Wed, 15 May 2024 13:04:25 +0200 Subject: [PATCH 11/18] doc --- include/ddc/kernels/splines/matrix.hpp | 13 ++++++------- include/ddc/kernels/splines/matrix_sparse.hpp | 11 ++++------- 2 files changed, 10 insertions(+), 14 deletions(-) diff --git a/include/ddc/kernels/splines/matrix.hpp b/include/ddc/kernels/splines/matrix.hpp index 7bf108885..5dcc93199 100644 --- a/include/ddc/kernels/splines/matrix.hpp +++ b/include/ddc/kernels/splines/matrix.hpp @@ -17,7 +17,7 @@ namespace ddc::detail { /** * @brief The parent class for a linear problem dedicated to the computation of a spline approximation. * - * It represents a square Matrix and provides methods to solve a multiple right-hand sides linear problem. + * It stores a square matrix and provides methods to solve a multiple right-hand sides linear problem. * Implementations may have different storage formats, filling methods and multiple right-hand sides linear solvers. */ template @@ -38,7 +38,7 @@ class SplinesLinearProblem virtual ~SplinesLinearProblem() = default; /** - * @brief Get an element of the matrix at indexes i, j. It must not be called after `factorize`. + * @brief Get an element of the matrix at indexes i, j. It must not be called after `setup_solver`. * * @param i The row index of the desired element. * @param j The column index of the desired element. @@ -48,7 +48,7 @@ class SplinesLinearProblem virtual double get_element(int i, int j) const = 0; /** - * @brief Set an element of the matrix at indexes i, j. It must not be called after `factorize`. + * @brief Set an element of the matrix at indexes i, j. It must not be called after `setup_solver`. * * @param i The row index of the setted element. * @param j The column index of the setted element. @@ -57,16 +57,15 @@ class SplinesLinearProblem virtual void set_element(int i, int j, double aij) = 0; /** - * @brief Performs a pre-process operation on the solver. + * @brief Perform a pre-process operation on the solver. Must be called after filling the matrix. */ virtual void setup_solver() = 0; /** - * @brief Solve the multiple right-hand sides linear problem Ax=b inplace. + * @brief Solve the multiple right-hand sides linear problem Ax=b or its transposed version A^tx=b inplace. * * @param[in, out] multi_rhs A 2D Kokkos::View storing the multiple right-hand sides of the problem and receiving the corresponding solution. - * Important note: the convention is the reverse of the common matrix one, row number is second index and column number - * the first one. This means when solving `A x_i = b_i`, element `b_{ij}` is stored in `b(j, i)`. + * @param transpose Choose between the direct or transposed version of the linear problem. */ virtual void solve(MultiRHS b, bool transpose = false) const = 0; diff --git a/include/ddc/kernels/splines/matrix_sparse.hpp b/include/ddc/kernels/splines/matrix_sparse.hpp index 82b49069e..0d8931b9d 100644 --- a/include/ddc/kernels/splines/matrix_sparse.hpp +++ b/include/ddc/kernels/splines/matrix_sparse.hpp @@ -191,13 +191,11 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem } /** - * @brief A function called by factorize() to perform the pre-process operation. + * @brief Perform a pre-process operation on the solver. Must be called after filling the matrix. * * Removes the zeros from the CSR object and instantiate a Ginkgo solver. It also constructs a transposed version of the solver. * * The stopping criterion is a reduction factor ||x||/||b||<1e-15 with 1000 maximum iterations. - * - * @return The error code of the function. */ void setup_solver() override { @@ -234,15 +232,14 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem } /** - * @brief A function called by solve_inplace() and similar functions to actually perform the linear solve operation. + * @brief Solve the multiple right-hand sides linear problem Ax=b or its transposed version A^tx=b inplace. * * The solver method is currently Bicgstab on CPU Serial and GPU and Gmres on OMP (because of Ginkgo issue #1563). * * Multiple right-hand sides are sliced in chunks of size cols_per_chunk which are passed one-after-the-other to Ginkgo. * - * @param b A double* to a contiguous array containing the (eventually multiple) right-hand sides. - * @param transpose A character identifying if the normal ('N') or transposed ('T') linear system is solved. - * @param n_equations The number of multiple right-hand sides (number of columns in b). + * @param[in, out] multi_rhs A 2D Kokkos::View storing the multiple right-hand sides of the problem and receiving the corresponding solution. + * @param transpose Choose between the direct or transposed version of the linear problem. */ void solve(typename SplinesLinearProblem::MultiRHS b, bool const transpose) const override From 126b0bc9679940d1c59ddee2b989e9a6eee645a9 Mon Sep 17 00:00:00 2001 From: blegouix Date: Wed, 15 May 2024 13:07:19 +0200 Subject: [PATCH 12/18] minor --- include/ddc/kernels/splines/matrix_sparse.hpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/include/ddc/kernels/splines/matrix_sparse.hpp b/include/ddc/kernels/splines/matrix_sparse.hpp index 0d8931b9d..400c4da4f 100644 --- a/include/ddc/kernels/splines/matrix_sparse.hpp +++ b/include/ddc/kernels/splines/matrix_sparse.hpp @@ -244,21 +244,20 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem void solve(typename SplinesLinearProblem::MultiRHS b, bool const transpose) const override { - std::size_t const n_equations = b.extent(1); assert(int(b.extent(0)) == size()); std::shared_ptr const gko_exec = m_solver->get_executor(); std::shared_ptr const convergence_logger = gko::log::Convergence::create(); - int const main_chunk_size = std::min(m_cols_per_chunk, (int)n_equations); + int const main_chunk_size = std::min(m_cols_per_chunk, (int)b.extent(1)); Kokkos::View const x("", size(), main_chunk_size); - int const iend = (n_equations + main_chunk_size - 1) / main_chunk_size; + int const iend = (b.extent(1) + main_chunk_size - 1) / main_chunk_size; for (int i = 0; i < iend; ++i) { int const subview_begin = i * main_chunk_size; int const subview_end - = (i + 1 == iend) ? n_equations : (subview_begin + main_chunk_size); + = (i + 1 == iend) ? b.extent(1) : (subview_begin + main_chunk_size); auto const b_chunk = Kokkos::subview(b, Kokkos::ALL, Kokkos::pair(subview_begin, subview_end)); From b66d43fc6d874efbdf8a01845b5392917ac124cb Mon Sep 17 00:00:00 2001 From: blegouix Date: Wed, 15 May 2024 13:33:04 +0200 Subject: [PATCH 13/18] std::size_t --- benchmarks/splines.cpp | 6 ++-- include/ddc/kernels/splines/matrix.hpp | 16 ++++----- include/ddc/kernels/splines/matrix_maker.hpp | 2 +- include/ddc/kernels/splines/matrix_sparse.hpp | 33 ++++++++++--------- .../ddc/kernels/splines/spline_builder.hpp | 6 ++-- .../ddc/kernels/splines/spline_builder_2d.hpp | 2 +- 6 files changed, 34 insertions(+), 31 deletions(-) diff --git a/benchmarks/splines.cpp b/benchmarks/splines.cpp index 18b61a10b..08bb55f83 100644 --- a/benchmarks/splines.cpp +++ b/benchmarks/splines.cpp @@ -185,15 +185,15 @@ static void characteristics_advection(benchmark::State& state) #ifdef KOKKOS_ENABLE_CUDA std::string chip = "gpu"; -int cols_per_chunk_ref = 65535; +std::size_t cols_per_chunk_ref = 65535; unsigned int preconditionner_max_block_size_ref = 1u; #elif defined(KOKKOS_ENABLE_OPENMP) std::string chip = "cpu"; -int cols_per_chunk_ref = 8192; +std::size_t cols_per_chunk_ref = 8192; unsigned int preconditionner_max_block_size_ref = 32u; #elif defined(KOKKOS_ENABLE_SERIAL) std::string chip = "cpu"; -int cols_per_chunk_ref = 8192; +std::size_t cols_per_chunk_ref = 8192; unsigned int preconditionner_max_block_size_ref = 32u; #endif diff --git a/include/ddc/kernels/splines/matrix.hpp b/include/ddc/kernels/splines/matrix.hpp index 5dcc93199..76ccc32a5 100644 --- a/include/ddc/kernels/splines/matrix.hpp +++ b/include/ddc/kernels/splines/matrix.hpp @@ -15,9 +15,9 @@ namespace ddc::detail { /** - * @brief The parent class for a linear problem dedicated to the computation of a spline approximation. + * @brief The parent class for linear problems dedicated to the computation of spline approximations. * - * It stores a square matrix and provides methods to solve a multiple right-hand sides linear problem. + * Store a square matrix and provide method to solve a multiple right-hand sides linear problem. * Implementations may have different storage formats, filling methods and multiple right-hand sides linear solvers. */ template @@ -45,7 +45,7 @@ class SplinesLinearProblem * * @return The value of the element of the matrix. */ - virtual double get_element(int i, int j) const = 0; + virtual double get_element(std::size_t i, std::size_t j) const = 0; /** * @brief Set an element of the matrix at indexes i, j. It must not be called after `setup_solver`. @@ -54,7 +54,7 @@ class SplinesLinearProblem * @param j The column index of the setted element. * @param aij The value to set in the element of the matrix. */ - virtual void set_element(int i, int j, double aij) = 0; + virtual void set_element(std::size_t i, std::size_t j, double aij) = 0; /** * @brief Perform a pre-process operation on the solver. Must be called after filling the matrix. @@ -74,7 +74,7 @@ class SplinesLinearProblem * * @return The size of the matrix in one of its dimensions. */ - int size() const + std::size_t size() const { return m_size; } @@ -88,9 +88,9 @@ class SplinesLinearProblem **/ std::ostream& operator<<(std::ostream& os) { - int const n = size(); - for (int i = 0; i < n; ++i) { - for (int j = 0; j < n; ++j) { + std::size_t const n = size(); + for (std::size_t i = 0; i < n; ++i) { + for (std::size_t j = 0; j < n; ++j) { os << std::fixed << std::setprecision(3) << std::setw(10) << get_element(i, j); } os << std::endl; diff --git a/include/ddc/kernels/splines/matrix_maker.hpp b/include/ddc/kernels/splines/matrix_maker.hpp index f9dbf7373..06e76baad 100644 --- a/include/ddc/kernels/splines/matrix_maker.hpp +++ b/include/ddc/kernels/splines/matrix_maker.hpp @@ -36,7 +36,7 @@ class SplinesLinearProblemMaker template static std::unique_ptr> make_new_sparse( int const n, - std::optional cols_per_chunk = std::nullopt, + std::optional cols_per_chunk = std::nullopt, std::optional preconditionner_max_block_size = std::nullopt) { return std::make_unique const& gko_exec, KokkosVi * @return The default value for the parameter cols_per_chunk. */ template -int default_cols_per_chunk() noexcept +std::size_t default_cols_per_chunk() noexcept { #ifdef KOKKOS_ENABLE_SERIAL if (std::is_same_v) { @@ -129,7 +129,7 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem using SplinesLinearProblem::operator<<; private: - using matrix_sparse_type = gko::matrix::Csr; + using matrix_sparse_type = gko::matrix::Csr; #ifdef KOKKOS_ENABLE_OPENMP using solver_type = std::conditional_t< std::is_same_v, @@ -148,7 +148,7 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem std::shared_ptr m_solver; std::shared_ptr m_solver_tr; - int m_cols_per_chunk; // Maximum number of columns of B to be passed to a Ginkgo solver + std::size_t m_cols_per_chunk; // Maximum number of columns of B to be passed to a Ginkgo solver unsigned int m_preconditionner_max_block_size; // Maximum size of Jacobi-block preconditionner @@ -163,8 +163,8 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem * used by the block-Jacobi preconditionner. see default_preconditionner_max_block_size. */ explicit SplinesLinearProblemSparse( - const int mat_size, - std::optional cols_per_chunk = std::nullopt, + const std::size_t mat_size, + std::optional cols_per_chunk = std::nullopt, std::optional preconditionner_max_block_size = std::nullopt) : SplinesLinearProblem(mat_size) , m_cols_per_chunk(cols_per_chunk.value_or(default_cols_per_chunk())) @@ -178,14 +178,15 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem m_matrix_sparse = matrix_sparse_type::create(gko_exec, gko::dim<2>(mat_size, mat_size)); } - virtual double get_element([[maybe_unused]] int i, [[maybe_unused]] int j) const override + virtual double get_element([[maybe_unused]] std::size_t i, [[maybe_unused]] std::size_t j) + const override { throw std::runtime_error( "SplinesLinearProblemSparse::get_element() is not implemented because no API is " "provided by Ginkgo"); } - virtual void set_element(int i, int j, double aij) override + virtual void set_element(std::size_t i, std::size_t j, double aij) override { m_matrix_dense->at(i, j) = aij; } @@ -244,25 +245,27 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem void solve(typename SplinesLinearProblem::MultiRHS b, bool const transpose) const override { - assert(int(b.extent(0)) == size()); + assert(b.extent(0) == size()); std::shared_ptr const gko_exec = m_solver->get_executor(); std::shared_ptr const convergence_logger = gko::log::Convergence::create(); - int const main_chunk_size = std::min(m_cols_per_chunk, (int)b.extent(1)); + std::size_t const main_chunk_size = std::min(m_cols_per_chunk, b.extent(1)); Kokkos::View const x("", size(), main_chunk_size); - int const iend = (b.extent(1) + main_chunk_size - 1) / main_chunk_size; - for (int i = 0; i < iend; ++i) { - int const subview_begin = i * main_chunk_size; - int const subview_end + std::size_t const iend = (b.extent(1) + main_chunk_size - 1) / main_chunk_size; + for (std::size_t i = 0; i < iend; ++i) { + std::size_t const subview_begin = i * main_chunk_size; + std::size_t const subview_end = (i + 1 == iend) ? b.extent(1) : (subview_begin + main_chunk_size); auto const b_chunk = Kokkos::subview(b, Kokkos::ALL, Kokkos::pair(subview_begin, subview_end)); - auto const x_chunk - = Kokkos::subview(x, Kokkos::ALL, Kokkos::pair(0, subview_end - subview_begin)); + auto const x_chunk = Kokkos:: + subview(x, + Kokkos::ALL, + Kokkos::pair(std::size_t(0), subview_end - subview_begin)); Kokkos::deep_copy(x_chunk, b_chunk); diff --git a/include/ddc/kernels/splines/spline_builder.hpp b/include/ddc/kernels/splines/spline_builder.hpp index 84434ab21..47106b57a 100644 --- a/include/ddc/kernels/splines/spline_builder.hpp +++ b/include/ddc/kernels/splines/spline_builder.hpp @@ -144,7 +144,7 @@ class SplineBuilder public: explicit SplineBuilder( batched_interpolation_domain_type const& batched_interpolation_domain, - std::optional cols_per_chunk = std::nullopt, + std::optional cols_per_chunk = std::nullopt, std::optional preconditionner_max_block_size = std::nullopt) : m_batched_interpolation_domain(batched_interpolation_domain) , m_offset(compute_offset(interpolation_domain())) @@ -311,7 +311,7 @@ class SplineBuilder void allocate_matrix( int lower_block_size, int upper_block_size, - std::optional cols_per_chunk = std::nullopt, + std::optional cols_per_chunk = std::nullopt, std::optional preconditionner_max_block_size = std::nullopt); void build_matrix_system(); @@ -476,7 +476,7 @@ void SplineBuilder< allocate_matrix( [[maybe_unused]] int lower_block_size, [[maybe_unused]] int upper_block_size, - std::optional cols_per_chunk, + std::optional cols_per_chunk, std::optional preconditionner_max_block_size) { // Special case: linear spline diff --git a/include/ddc/kernels/splines/spline_builder_2d.hpp b/include/ddc/kernels/splines/spline_builder_2d.hpp index eebaf1f98..1dcf486d2 100644 --- a/include/ddc/kernels/splines/spline_builder_2d.hpp +++ b/include/ddc/kernels/splines/spline_builder_2d.hpp @@ -155,7 +155,7 @@ class SplineBuilder2D */ explicit SplineBuilder2D( batched_interpolation_domain_type const& batched_interpolation_domain, - std::optional cols_per_chunk = std::nullopt, + std::optional cols_per_chunk = std::nullopt, std::optional preconditionner_max_block_size = std::nullopt) : m_spline_builder1( batched_interpolation_domain, From 70eb7ee701a33a8eeb81384c741ea0268264b140 Mon Sep 17 00:00:00 2001 From: blegouix Date: Wed, 15 May 2024 15:15:17 +0200 Subject: [PATCH 14/18] rename files --- include/ddc/kernels/splines.hpp | 6 +++--- include/ddc/kernels/splines/spline_builder.hpp | 1 + .../splines/{matrix.hpp => splines_linear_solver.hpp} | 0 .../{matrix_maker.hpp => splines_linear_solver_maker.hpp} | 2 +- .../{matrix_sparse.hpp => splines_linear_solver_sparse.hpp} | 2 +- 5 files changed, 6 insertions(+), 5 deletions(-) rename include/ddc/kernels/splines/{matrix.hpp => splines_linear_solver.hpp} (100%) rename include/ddc/kernels/splines/{matrix_maker.hpp => splines_linear_solver_maker.hpp} (97%) rename include/ddc/kernels/splines/{matrix_sparse.hpp => splines_linear_solver_sparse.hpp} (99%) diff --git a/include/ddc/kernels/splines.hpp b/include/ddc/kernels/splines.hpp index fc67f83fd..ec1651f46 100644 --- a/include/ddc/kernels/splines.hpp +++ b/include/ddc/kernels/splines.hpp @@ -11,9 +11,6 @@ #include "splines/greville_interpolation_points.hpp" #include "splines/knots_as_interpolation_points.hpp" #include "splines/math_tools.hpp" -#include "splines/matrix.hpp" -#include "splines/matrix_maker.hpp" -#include "splines/matrix_sparse.hpp" #include "splines/null_extrapolation_rule.hpp" #include "splines/periodic_extrapolation_rule.hpp" #include "splines/spline_boundary_conditions.hpp" @@ -21,4 +18,7 @@ #include "splines/spline_builder_2d.hpp" #include "splines/spline_evaluator.hpp" #include "splines/spline_evaluator_2d.hpp" +#include "splines/splines_linear_solver.hpp" +#include "splines/splines_linear_solver_maker.hpp" +#include "splines/splines_linear_solver_sparse.hpp" #include "splines/view.hpp" diff --git a/include/ddc/kernels/splines/spline_builder.hpp b/include/ddc/kernels/splines/spline_builder.hpp index 47106b57a..e445785be 100644 --- a/include/ddc/kernels/splines/spline_builder.hpp +++ b/include/ddc/kernels/splines/spline_builder.hpp @@ -8,6 +8,7 @@ #include "ddc/kokkos_allocator.hpp" #include "deriv.hpp" +#include "splines_linear_solver_maker.hpp" namespace ddc { enum class SplineSolver { diff --git a/include/ddc/kernels/splines/matrix.hpp b/include/ddc/kernels/splines/splines_linear_solver.hpp similarity index 100% rename from include/ddc/kernels/splines/matrix.hpp rename to include/ddc/kernels/splines/splines_linear_solver.hpp diff --git a/include/ddc/kernels/splines/matrix_maker.hpp b/include/ddc/kernels/splines/splines_linear_solver_maker.hpp similarity index 97% rename from include/ddc/kernels/splines/matrix_maker.hpp rename to include/ddc/kernels/splines/splines_linear_solver_maker.hpp index 06e76baad..56b284f66 100644 --- a/include/ddc/kernels/splines/matrix_maker.hpp +++ b/include/ddc/kernels/splines/splines_linear_solver_maker.hpp @@ -7,7 +7,7 @@ #include #include -#include "matrix_sparse.hpp" +#include "splines_linear_solver_sparse.hpp" namespace ddc::detail { diff --git a/include/ddc/kernels/splines/matrix_sparse.hpp b/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp similarity index 99% rename from include/ddc/kernels/splines/matrix_sparse.hpp rename to include/ddc/kernels/splines/splines_linear_solver_sparse.hpp index f160879fe..3636eefb4 100644 --- a/include/ddc/kernels/splines/matrix_sparse.hpp +++ b/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp @@ -15,7 +15,7 @@ #include #include "ginkgo_executors.hpp" -#include "matrix.hpp" +#include "splines_linear_solver.hpp" namespace ddc::detail { From 63d7c723cbe748d66be50bb037b72e3b2e8d19d9 Mon Sep 17 00:00:00 2001 From: blegouix Date: Thu, 16 May 2024 09:10:21 +0200 Subject: [PATCH 15/18] attempt to fix clang --- include/ddc/kernels/splines/splines_linear_solver_sparse.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp b/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp index 3636eefb4..00492a0b7 100644 --- a/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp +++ b/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp @@ -124,7 +124,7 @@ template class SplinesLinearProblemSparse : public SplinesLinearProblem { public: - using SplinesLinearProblem::MultiRHS; + using typename SplinesLinearProblem::MultiRHS; using SplinesLinearProblem::size; using SplinesLinearProblem::operator<<; @@ -242,8 +242,7 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem * @param[in, out] multi_rhs A 2D Kokkos::View storing the multiple right-hand sides of the problem and receiving the corresponding solution. * @param transpose Choose between the direct or transposed version of the linear problem. */ - void solve(typename SplinesLinearProblem::MultiRHS b, bool const transpose) - const override + void solve(MultiRHS b, bool const transpose) const override { assert(b.extent(0) == size()); From 6033f089a786f793080abb3b8566a0781fd34ff3 Mon Sep 17 00:00:00 2001 From: blegouix Date: Fri, 17 May 2024 12:53:21 +0200 Subject: [PATCH 16/18] restore SplinesLinearSolverSparse::get_element --- .../splines/splines_linear_solver_sparse.hpp | 4 +--- tests/splines/matrix.cpp | 16 ++++++---------- 2 files changed, 7 insertions(+), 13 deletions(-) diff --git a/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp b/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp index 00492a0b7..bf33e94a4 100644 --- a/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp +++ b/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp @@ -181,9 +181,7 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem virtual double get_element([[maybe_unused]] std::size_t i, [[maybe_unused]] std::size_t j) const override { - throw std::runtime_error( - "SplinesLinearProblemSparse::get_element() is not implemented because no API is " - "provided by Ginkgo"); + return m_matrix_dense->at(i, j); } virtual void set_element(std::size_t i, std::size_t j, double aij) override diff --git a/tests/splines/matrix.cpp b/tests/splines/matrix.cpp index 39c3e20f4..4352e4b9b 100644 --- a/tests/splines/matrix.cpp +++ b/tests/splines/matrix.cpp @@ -30,11 +30,12 @@ void fill_identity( } } -/* -void copy_matrix(ddc::detail::SplinesLinearProblem::MultiRHS copy, std::unique_ptr>& mat) +void copy_matrix( + ddc::detail::SplinesLinearProblem::MultiRHS copy, + std::unique_ptr>& mat) { - assert(mat->get_size() == int(copy.extent(0))); - assert(mat->get_size() == int(copy.extent(1))); + assert(mat->size() == int(copy.extent(0))); + assert(mat->size() == int(copy.extent(1))); for (std::size_t i(0); i < copy.extent(0); ++i) { for (std::size_t j(0); j < copy.extent(1); ++j) { @@ -42,7 +43,6 @@ void copy_matrix(ddc::detail::SplinesLinearProblem::MultiRHS matrix, @@ -101,16 +101,12 @@ TEST_P(MatrixSizesFixture, Sparse) for (std::size_t j(0); j < N; ++j) { if (i == j) { matrix->set_element(i, j, 3. / 4); - val(i, j) = 3. / 4; } else if (std::abs((std::ptrdiff_t)(j - i)) <= (std::ptrdiff_t)k) { matrix->set_element(i, j, -(1. / 4) / k); - val(i, j) = -(1. / 4) / k; - } else { - val(i, j) = 0.; } } } - // copy_matrix(val, matrix); // copy_matrix is not available for sparse matrix because of a limitation of Ginkgo API (get_element is not implemented). The workaround is to fill val directly in the loop + copy_matrix(val, matrix); matrix->setup_solver(); From 6a560906d026faa04f8b9f571c3717a18d72fdea Mon Sep 17 00:00:00 2001 From: Thomas Padioleau Date: Fri, 17 May 2024 16:18:01 +0200 Subject: [PATCH 17/18] Update include/ddc/kernels/splines/splines_linear_solver_sparse.hpp --- include/ddc/kernels/splines/splines_linear_solver_sparse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp b/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp index bf33e94a4..7ae29c74b 100644 --- a/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp +++ b/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp @@ -178,7 +178,7 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem m_matrix_sparse = matrix_sparse_type::create(gko_exec, gko::dim<2>(mat_size, mat_size)); } - virtual double get_element([[maybe_unused]] std::size_t i, [[maybe_unused]] std::size_t j) + virtual double get_element(std::size_t i, std::size_t j) const override { return m_matrix_dense->at(i, j); From 4dd4564fba85284c2a3198512ae3e6137a167728 Mon Sep 17 00:00:00 2001 From: Thomas Padioleau Date: Fri, 17 May 2024 16:19:53 +0200 Subject: [PATCH 18/18] Update include/ddc/kernels/splines/splines_linear_solver_sparse.hpp --- include/ddc/kernels/splines/splines_linear_solver_sparse.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp b/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp index 7ae29c74b..3d38ccb20 100644 --- a/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp +++ b/include/ddc/kernels/splines/splines_linear_solver_sparse.hpp @@ -178,8 +178,7 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem m_matrix_sparse = matrix_sparse_type::create(gko_exec, gko::dim<2>(mat_size, mat_size)); } - virtual double get_element(std::size_t i, std::size_t j) - const override + virtual double get_element(std::size_t i, std::size_t j) const override { return m_matrix_dense->at(i, j); }