Add distributed MPI support for solving classes (#248)

Author: TLCFEM

.clang-format

@@ -159,7 +159,7 @@ ReferenceAlignment: Pointer
 ReflowComments: true
 RemoveBracesLLVM: false
 RemoveSemicolon: false
-RequiresClausePosition: OwnLine
+RequiresClausePosition: SingleLine
 RequiresExpressionIndentation: OuterScope
 SeparateDefinitionBlocks: Leave
 ShortNamespaceLines: 1

CMakeLists.txt

@@ -141,6 +141,7 @@ if (USE_MPI)
     elseif (MPI_C_INCLUDE_DIRS_LOWER MATCHES "mpich")
         set(MKL_MPI mpich)
     endif ()
+    add_compile_definitions(SUANPAN_DISTRIBUTED)
 endif ()
 
 if (USE_MKL)

Domain/DomainState.cpp

@@ -34,8 +34,6 @@
 #include "ILU.hpp"
 #include "Jacobi.hpp"
 
-template<typename T, typename U> concept ArmaContainer = std::is_floating_point_v<U> && (std::is_convertible_v<T, Mat<U>> || std::is_convertible_v<T, SpMat<U>>) ;
-
 template<sp_d T> class MetaMat;
 
 template<sp_d T> class op_add {
@@ -199,9 +197,9 @@ template<sp_d T> class MetaMat {
 
     virtual void operator*=(T) = 0;
 
-    template<ArmaContainer<T> C> int solve(Mat<T>& X, C&& B) { return IterativeSolver::NONE == this->setting.iterative_solver ? this->direct_solve(X, std::forward<C>(B)) : this->iterative_solve(X, std::forward<C>(B)); }
+    template<typename C> requires is_arma_mat<T, C> int solve(Mat<T>& X, C&& B) { return IterativeSolver::NONE == this->setting.iterative_solver ? this->direct_solve(X, std::forward<C>(B)) : this->iterative_solve(X, std::forward<C>(B)); }
 
-    template<ArmaContainer<T> C> Mat<T> solve(C&& B) {
+    template<typename C> requires is_arma_mat<T, C> Mat<T> solve(C&& B) {
         Mat<T> X;
 
         if(SUANPAN_SUCCESS != this->solve(X, std::forward<C>(B))) throw std::runtime_error("fail to solve the system");

Domain/MetaMat/MetaMat.hpp

@@ -257,25 +257,32 @@ Element::Element(const unsigned T, const unsigned NN, const unsigned ND, uvec&&
 Element::Element(const unsigned T, const unsigned NN, const unsigned ND, uvec&& NT, uvec&& MT, const bool F, const MaterialType MTP, std::vector<DOF>&& DI)
     : DataElement{std::move(NT), std::move(MT), uvec{}, F}
     , ElementBase(T)
+    , Distributed(T)
     , num_node(NN)
     , num_dof(ND)
     , material_type(MTP)
     , section_type(SectionType::D0)
-    , dof_identifier(std::move(DI)) { suanpan_assert([&] { if(!dof_identifier.empty() && num_dof != dof_identifier.size()) throw invalid_argument("size of dof identifier must meet number of dofs"); }); }
+    , dof_identifier(std::move(DI)) {
+    suanpan_assert([&] { if(!dof_identifier.empty() && num_dof != dof_identifier.size()) throw invalid_argument("size of dof identifier must meet number of dofs"); });
+}
 
 Element::Element(const unsigned T, const unsigned NN, const unsigned ND, uvec&& NT, uvec&& ST, const bool F, const SectionType STP, std::vector<DOF>&& DI)
     : DataElement{std::move(NT), uvec{}, std::move(ST), F}
     , ElementBase(T)
+    , Distributed(T)
     , num_node(NN)
     , num_dof(ND)
     , material_type(MaterialType::D0)
     , section_type(STP)
-    , dof_identifier(std::move(DI)) { suanpan_assert([&] { if(!dof_identifier.empty() && num_dof != dof_identifier.size()) throw invalid_argument("size of dof identifier must meet number of dofs"); }); }
+    , dof_identifier(std::move(DI)) {
+    suanpan_assert([&] { if(!dof_identifier.empty() && num_dof != dof_identifier.size()) throw invalid_argument("size of dof identifier must meet number of dofs"); });
+}
 
 // for contact elements that use node groups
 Element::Element(const unsigned T, const unsigned ND, uvec&& GT)
     : DataElement{std::move(GT), {}, {}, false}
     , ElementBase(T)
+    , Distributed(T)
     , num_node(static_cast<unsigned>(-1))
     , num_dof(ND)
     , use_group(true)
@@ -286,6 +293,7 @@ Element::Element(const unsigned T, const unsigned ND, uvec&& GT)
 Element::Element(const unsigned T, const unsigned ND, const unsigned ET, const unsigned NT)
     : DataElement{{NT}, {}, {}, false}
     , ElementBase(T)
+    , Distributed(ET)
     , num_node(static_cast<unsigned>(-1))
     , num_dof(ND)
     , use_other(ET)

Element/Element.cpp

@@ -29,6 +29,7 @@
 #define ELEMENT_H
 
 #include <Element/ElementBase.h>
+#include <MPI/Distributed.h>
 
 enum class MaterialType : unsigned;
 enum class SectionType : unsigned;
@@ -114,7 +115,7 @@ struct DataElement {
     const double characteristic_length = 1.;
 };
 
-class Element : protected DataElement, public ElementBase {
+class Element : protected DataElement, public ElementBase, public Distributed {
     const unsigned num_node;                      // number of nodes
     const unsigned num_dof;                       // number of DoFs
     const unsigned num_size = num_dof * num_node; // number of size

Element/Element.h

@@ -17,3 +17,6 @@ add_dependencies(tester.pardiso solver.pardiso)
 add_dependencies(suanPan solver.pardiso)
 
 install(TARGETS solver.pardiso DESTINATION bin)
+
+add_executable(distributed_obj distributed_obj.cpp)
+target_link_libraries(distributed_obj MKL::MKL_SCALAPACK MPI::MPI_C MPI::MPI_CXX MPI::MPI_Fortran)

MPI/CMakeLists.txt

@@ -0,0 +1,63 @@
+/*******************************************************************************
+ * Copyright (C) 2017-2025 Theodore Chang
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ ******************************************************************************/
+
+#ifndef DISTRIBUTED_H
+#define DISTRIBUTED_H
+
+#include <suanPan.h>
+
+class Distributed {
+    static constexpr int root_rank{0};
+
+    static auto assign_process(const int obj_tag) { return obj_tag % comm_size; }
+
+public:
+    const int tag, process_rank;
+
+    const bool is_local;
+
+    explicit Distributed(const int obj_tag)
+        : tag(obj_tag)
+        , process_rank(assign_process(tag))
+        , is_local(comm_rank == process_rank) {}
+
+#ifdef SUANPAN_DISTRIBUTED
+    /**
+     * @brief Performs a gather operation on a distributed matrix object.
+     *
+     * This function initiates a non-blocking gather operation. If the calling process is the root process,
+     * it receives data from the specified process. If the calling process is not the root process, it sends
+     * its data to the root process.
+     *
+     * @tparam DT The data type of the matrix elements.
+     * @param object The matrix object to be gathered.
+     * @return An optional non-blocking request handle for the gather operation.
+     */
+    template<mpl_data_t DT> std::optional<mpl::irequest> gather(const Mat<DT>& object) {
+        if(root_rank == comm_rank) {
+            if(!is_local) return comm_world.irecv(const_cast<DT*>(object.memptr()), mpl::contiguous_layout<DT>{object.n_elem}, process_rank, mpl::tag_t{tag});
+        }
+        else if(is_local) return comm_world.isend(object.memptr(), mpl::contiguous_layout<DT>{object.n_elem}, root_rank, mpl::tag_t{tag});
+
+        return {};
+    }
+#else
+    template<mpl_data_t DT> static auto gather(const Mat<DT>&) {}
+#endif
+};
+
+#endif // DISTRIBUTED_H

MPI/Distributed.h

@@ -0,0 +1,39 @@
+#include "Distributed.h"
+
+#include <mpl/mpl.hpp>
+
+class Object : public Distributed {
+    mat data{-999, -999};
+
+public:
+    explicit Object(int tag)
+        : Distributed(tag) {}
+
+    auto generate() {
+        if(is_local) data.fill(fill::randn);
+    }
+
+    auto gather_data() { return gather(data); }
+
+    auto print() { printf("Object %d on process %d data: %+.6f %+.6f.\n", tag, mpl::environment::comm_world().rank(), data(0), data(1)); }
+};
+
+int main() {
+    arma_rng::set_seed_random();
+
+    std::vector<Object> objects;
+
+    for(auto i = 0; i < 100; ++i) objects.emplace_back(i);
+
+    for(auto& i : objects) i.generate();
+
+    mpl::irequest_pool requests;
+    for(auto& i : objects)
+        if(auto req = i.gather_data(); req.has_value()) requests.push(std::move(req).value());
+    requests.waitall();
+
+    if(0 == mpl::environment::comm_world().rank())
+        for(auto& i : objects) i.print();
+
+    return 0;
+}

MPI/distributed_obj.cpp

@@ -16,13 +16,16 @@
  ******************************************************************************/
 
 #include "Integrator.h"
+
 #include <Domain/DomainBase.h>
 #include <Domain/Factory.hpp>
 
 Integrator::Integrator(const unsigned T)
     : Tag(T) {}
 
-void Integrator::set_domain(const weak_ptr<DomainBase>& D) { if(database.lock() != D.lock()) database = D; }
+void Integrator::set_domain(const weak_ptr<DomainBase>& D) {
+    if(database.lock() != D.lock()) database = D;
+}
 
 shared_ptr<DomainBase> Integrator::get_domain() const { return database.lock(); }
 
@@ -260,13 +263,69 @@ mat Integrator::solve(sp_mat&& B) {
     return X;
 }
 
-int Integrator::solve(mat& X, const mat& B) { return database.lock()->get_factory()->get_stiffness()->solve(X, B); }
+int Integrator::solve(mat& X, const mat& B) {
+    int info{0};
+    // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+    if(0 == comm_rank) info = database.lock()->get_factory()->get_stiffness()->solve(X, B);
+
+    if(SUANPAN_SUCCESS == bcast_from_root(info)) {
+        // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+        // ReSharper disable once CppDFAUnreachableCode
+        if(0 != comm_rank) X.set_size(B.n_rows, B.n_cols);
+        bcast_from_root(X);
+    }
+
+    return info;
+}
+
+int Integrator::solve(mat& X, const sp_mat& B) {
+    int info{0};
+    // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+    if(0 == comm_rank) info = database.lock()->get_factory()->get_stiffness()->solve(X, B);
 
-int Integrator::solve(mat& X, const sp_mat& B) { return database.lock()->get_factory()->get_stiffness()->solve(X, B); }
+    if(SUANPAN_SUCCESS == bcast_from_root(info)) {
+        // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+        // ReSharper disable once CppDFAUnreachableCode
+        if(0 != comm_rank) X.set_size(B.n_rows, B.n_cols);
+        bcast_from_root(X);
+    }
+
+    return info;
+}
 
-int Integrator::solve(mat& X, mat&& B) { return database.lock()->get_factory()->get_stiffness()->solve(X, std::move(B)); }
+int Integrator::solve(mat& X, mat&& B) {
+    const auto n_rows = B.n_rows, n_cols = B.n_cols;
+
+    int info{0};
+    // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+    if(0 == comm_rank) info = database.lock()->get_factory()->get_stiffness()->solve(X, std::move(B));
+
+    if(SUANPAN_SUCCESS == bcast_from_root(info)) {
+        // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+        // ReSharper disable once CppDFAUnreachableCode
+        if(0 != comm_rank) X.set_size(n_rows, n_cols);
+        bcast_from_root(X);
+    }
+
+    return info;
+}
+
+int Integrator::solve(mat& X, sp_mat&& B) {
+    const auto n_rows = B.n_rows, n_cols = B.n_cols;
+
+    int info{0};
+    // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+    if(0 == comm_rank) info = database.lock()->get_factory()->get_stiffness()->solve(X, std::move(B));
+
+    if(SUANPAN_SUCCESS == bcast_from_root(info)) {
+        // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+        // ReSharper disable once CppDFAUnreachableCode
+        if(0 != comm_rank) X.set_size(n_rows, n_cols);
+        bcast_from_root(X);
+    }
 
-int Integrator::solve(mat& X, sp_mat&& B) { return database.lock()->get_factory()->get_stiffness()->solve(X, std::move(B)); }
+    return info;
+}
 
 /**
  * Avoid machine error accumulation.
@@ -363,13 +422,69 @@ void ExplicitIntegrator::update_from_ninja() {
     W->update_trial_acceleration_by(W->get_ninja());
 }
 
-int ExplicitIntegrator::solve(mat& X, const mat& B) { return get_domain()->get_factory()->get_mass()->solve(X, B); }
+int ExplicitIntegrator::solve(mat& X, const mat& B) {
+    int info{0};
+    // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+    if(0 == comm_rank) info = get_domain()->get_factory()->get_mass()->solve(X, B);
+
+    if(SUANPAN_SUCCESS == bcast_from_root(info)) {
+        // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+        // ReSharper disable once CppDFAUnreachableCode
+        if(0 != comm_rank) X.set_size(B.n_rows, B.n_cols);
+        bcast_from_root(X);
+    }
+
+    return info;
+}
+
+int ExplicitIntegrator::solve(mat& X, const sp_mat& B) {
+    int info{0};
+    // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+    if(0 == comm_rank) info = get_domain()->get_factory()->get_mass()->solve(X, B);
+
+    if(SUANPAN_SUCCESS == bcast_from_root(info)) {
+        // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+        // ReSharper disable once CppDFAUnreachableCode
+        if(0 != comm_rank) X.set_size(B.n_rows, B.n_cols);
+        bcast_from_root(X);
+    }
+
+    return info;
+}
 
-int ExplicitIntegrator::solve(mat& X, const sp_mat& B) { return get_domain()->get_factory()->get_mass()->solve(X, B); }
+int ExplicitIntegrator::solve(mat& X, mat&& B) {
+    const auto n_rows = B.n_rows, n_cols = B.n_cols;
 
-int ExplicitIntegrator::solve(mat& X, mat&& B) { return get_domain()->get_factory()->get_mass()->solve(X, std::move(B)); }
+    int info{0};
+    // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+    if(0 == comm_rank) info = get_domain()->get_factory()->get_mass()->solve(X, std::move(B));
 
-int ExplicitIntegrator::solve(mat& X, sp_mat&& B) { return get_domain()->get_factory()->get_mass()->solve(X, std::move(B)); }
+    if(SUANPAN_SUCCESS == bcast_from_root(info)) {
+        // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+        // ReSharper disable once CppDFAUnreachableCode
+        if(0 != comm_rank) X.set_size(n_rows, n_cols);
+        bcast_from_root(X);
+    }
+
+    return info;
+}
+
+int ExplicitIntegrator::solve(mat& X, sp_mat&& B) {
+    const auto n_rows = B.n_rows, n_cols = B.n_cols;
+
+    int info{0};
+    // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+    if(0 == comm_rank) info = get_domain()->get_factory()->get_mass()->solve(X, std::move(B));
+
+    if(SUANPAN_SUCCESS == bcast_from_root(info)) {
+        // ReSharper disable once CppIfCanBeReplacedByConstexprIf
+        // ReSharper disable once CppDFAUnreachableCode
+        if(0 != comm_rank) X.set_size(n_rows, n_cols);
+        bcast_from_root(X);
+    }
+
+    return info;
+}
 
 vec ExplicitIntegrator::from_incre_velocity(const vec&, const uvec&) { throw invalid_argument("support velocity cannot be used with explicit integrator"); }