Matrix classes doc

include/ddc/kernels/splines/matrix.hpp

@@ -14,19 +14,46 @@
 
 namespace ddc::detail {
 
+/**
+ * @brief The parent class for matrices.
+ *
+ * It represents a square Matrix. Implementations may have different storage formats, filling methods and multiple right-hand sides linear solvers.
+ */
 class Matrix
 {
     int m_n;
 
 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. It must not be called after `factorize`.
+     *
+     * @param i The row 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. 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.
+     * @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 +72,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 +89,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 +108,15 @@ class Matrix
         return b;
     }
 
+    /**
+     * @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);
@@ -81,6 +129,15 @@ class Matrix
         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.
+     * 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 +150,15 @@ 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 <class... Args>
     Kokkos::View<double**, Args...> solve_batch_inplace(
             Kokkos::View<double**, Args...> const bx) const
@@ -107,11 +173,23 @@ class Matrix
         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
     {
         return m_n;
     }
 
+    /**
+     * @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.
+     *
+     * @return The stream in which the matrix is printed.
+     **/
     std::ostream& operator<<(std::ostream& os)
     {
         int const n = get_size();
@@ -125,8 +203,21 @@ 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 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;
 };
 

include/ddc/kernels/splines/matrix_maker.hpp

@@ -11,9 +11,28 @@
 
 namespace ddc::detail {
 
+/**
+ * @brief A static factory for Matrix classes.
+ */
 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 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 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.
+     *
+     * @return The Matrix instance.
+     */
     template <typename ExecSpace>
     static std::unique_ptr<Matrix> make_new_sparse(
             int const n,

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,15 @@ auto to_gko_dense(std::shared_ptr<const gko::Executor> const& gko_exec, KokkosVi
                    view.stride_0());
 }
 
+/**
+ * @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 parameter cols_per_chunk.
+ */
 template <class ExecSpace>
 int default_cols_per_chunk() noexcept
 {
@@ -67,6 +78,15 @@ int default_cols_per_chunk() noexcept
     return 1;
 }
 
+/**
+ * @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 parameter preconditionner_max_block_size.
+ */
 template <class ExecSpace>
 unsigned int default_preconditionner_max_block_size() noexcept
 {
@@ -93,7 +113,13 @@ unsigned int default_preconditionner_max_block_size() noexcept
     return 1u;
 }
 
-// Matrix class for sparse storage and iterative solve
+/**
+ * @brief A sparse Matrix.
+ *
+ * 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 ExecSpace>
 class Matrix_Sparse : public Matrix
 {
@@ -121,7 +147,15 @@ 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 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 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(
             const int mat_size,
             std::optional<int> cols_per_chunk = std::nullopt,
@@ -149,6 +183,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 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.
+     */
     int factorize_method() override
     {
         // Remove zeros
@@ -185,6 +228,19 @@ 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. 
+     * @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();