OpenMP implementation #210

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open

rvanvenetie wants to merge 65 commits into master from cpp/spacetime-omp

src/CMakeLists.txt

-Original file line number
+Diff line change
@@ Expand Up / @@ -22,6 +22,13 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS ON) @@
     #     set(CMAKE_CXX_COMPILER_LAUNCHER "${CCACHE_PROGRAM}")
     # endif()
+    # Set OpenMP
+    find_package(OpenMP)
+    include(ProcessorCount)
+    ProcessorCount(MAX_NUMBER_THREADS)
+    message("Setting maximum number of threads to ${MAX_NUMBER_THREADS}")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DMAX_NUMBER_THREADS=${MAX_NUMBER_THREADS}")
     ### Download GoogleTest
     include(FetchContent)
     FetchContent_Declare(
@@ Expand Down @@

src/applications/adaptive.cpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -114,10 +114,12 @@ int main(int argc, char* argv[]) { @@
       std::string problem, domain;
       size_t initial_refines = 0;
       size_t max_dofs = 0;
+      size_t num_threads = 1;
       bool calculate_condition_numbers = false;
       bool print_centers = false;
       bool print_sampling = false;
-      bool print_time_apply = false;
+      bool print_time_apply = true;
+      bool print_bilforms = false;
       std::vector<double> print_time_slices;
       boost::program_options::options_description problem_optdesc(
           "Problem options");
@@ Expand All / @@ -133,7 +135,9 @@ int main(int argc, char* argv[]) { @@
           "print_sampling", po::value<bool>(&print_sampling))(
           "print_time_slices",
           po::value<std::vector<double>>(&print_time_slices)->multitoken())(
-          "print_time_apply", po::value<bool>(&print_time_apply));
+          "print_time_apply", po::value<bool>(&print_time_apply))(
+          "print_bilforms", po::value<bool>(&print_bilforms))(
+          "num_threads", po::value<size_t>(&num_threads));
       std::sort(print_time_slices.begin(), print_time_slices.end());
@@ Expand Down Expand Up / @@ -165,10 +169,20 @@ int main(int argc, char* argv[]) { @@
       po::store(po::command_line_parser(argc, argv).options(cmdline_options).run(),
                 vm);
       po::notify(vm);
+      assert(num_threads > 0);
+      if (num_threads > 1 && adapt_opts.use_cache) {
+        std::cout << "Multithreading is only enabled for no-cache." << std::endl;
+        return 1;
+      }
+      assert(num_threads <= omp_get_max_threads());
+      assert(num_threads <= MAX_NUMBER_THREADS);
+      omp_set_num_threads(num_threads);
       std::cout << "Problem options:" << std::endl;
       std::cout << "\tProblem: " << problem << std::endl;
       std::cout << "\tDomain: " << domain
                 << "; initial-refines: " << initial_refines << std::endl;
+      std::cout << "\tNumber-threads: " << num_threads << std::endl;
       std::cout << std::endl;
       std::cout << adapt_opts << std::endl;
@@ Expand Down Expand Up / @@ -215,9 +229,13 @@ int main(int argc, char* argv[]) { @@
         // A slight overestimate.
         ndof_Xd = vec_Xd->Bfs().size();
+        size_t ndof_Xd_time = vec_Xd->Project_0()->Bfs().size();
+        size_t ndof_Xd_space = vec_Xd->Project_1()->Bfs().size();
         size_t ndof_Xdd = heat_eq.vec_Xdd()->Bfs().size();
         size_t ndof_Ydd = heat_eq.vec_Ydd()->Bfs().size();
         std::cout << "iter: " << ++iter << "\n\tXDelta-size: " << ndof_Xd
+                  << "\n\tXDelta-space-size: " << ndof_Xd_space
+                  << "\n\tXDelta-time-size: " << ndof_Xd_time
                   << "\n\tXDelta-Gradedness: "
                   << vec_Xd->Gradedness(&max_gradedness)
                   << "\n\tXDeltaDelta-size: " << ndof_Xdd
@@ Expand Down Expand Up / @@ -317,18 +335,25 @@ int main(int argc, char* argv[]) { @@
         if (print_time_apply) {
           auto heat_d_dd = heat_eq.heat_d_dd();
-          std::cout << "\n\tA-time-per-apply: " << heat_d_dd->A()->TimePerApply()
-                    << "\n\tB-time-per-apply: " << heat_d_dd->B()->TimePerApply()
-                    << "\n\tBT-time-per-apply: " << heat_d_dd->BT()->TimePerApply()
-                    << "\n\tG-time-per-apply: " << heat_d_dd->G()->TimePerApply()
-                    << "\n\tP_Y-time-per-apply: "
-                    << heat_d_dd->P_Y()->TimePerApply()
-                    << "\n\tP_X-time-per-apply: "
-                    << heat_d_dd->P_X()->TimePerApply()
-                    << "\n\tS-time-per-apply: " << heat_d_dd->S()->TimePerApply()
-                    << "\n\ttotal-time-apply: " << heat_d_dd->TotalTimeApply()
-                    << "\n\ttotal-time-construct: "
-                    << heat_d_dd->TotalTimeConstruct() << std::flush;
+          std::cout
+              << "\n\tA-time-per-apply: " << heat_d_dd->A()->TimePerApply()
+              << "\n\tB-time-per-apply: " << heat_d_dd->B()->TimePerApply()
+              << "\n\tB-A-time-per-apply: " << heat_d_dd->B()->A()->TimePerApply()
+              << "\n\tB-B-time-per-apply: " << heat_d_dd->B()->B()->TimePerApply()
+              << "\n\tBT-time-per-apply: " << heat_d_dd->BT()->TimePerApply()
+              << "\n\tG-time-per-apply: " << heat_d_dd->G()->TimePerApply()
+              << "\n\tP_Y-time-per-apply: " << heat_d_dd->P_Y()->TimePerApply()
+              << "\n\tP_X-time-per-apply: " << heat_d_dd->P_X()->TimePerApply()
+              << "\n\tS-time-per-apply: " << heat_d_dd->S()->TimePerApply()
+              << "\n\ttotal-time-apply: " << heat_d_dd->TotalTimeApply()
+              << "\n\ttotal-time-construct: " << heat_d_dd->TotalTimeConstruct()
+              << std::flush;
+        }
+        if (print_bilforms) {
+          auto heat_d_dd = heat_eq.heat_d_dd();
+          std::cout << "\n\tB-A-bilforms: " << heat_d_dd->B()->A()->Information()
+                    << "\n\tP_Y-bilforms: " << heat_d_dd->P_Y()->Information()
+                    << std::flush;
         }
         if (print_centers) {
@@ Expand Down @@

src/applications/uniform.cpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -53,10 +53,11 @@ int main(int argc, char* argv[]) { @@
       size_t initial_refines = 0;
       size_t max_level = 0;
       size_t max_dofs = 0;
+      size_t num_threads = 1;
       std::string refine;
       bool calculate_condition_PY = false;
       bool calculate_condition_PX = false;
-      bool print_time_apply = false;
+      bool print_time_apply = true;
       bool print_centers = false;
       double solve_rtol = 1e-5;
       boost::program_options::options_description problem_optdesc(
@@ Expand All / @@ -70,6 +71,7 @@ int main(int argc, char* argv[]) { @@
               ->default_value(std::numeric_limits<std::size_t>::max()))(
           "max_dofs", po::value<size_t>(&max_dofs)->default_value(
                           std::numeric_limits<std::size_t>::max()))(
+          "num_threads", po::value<size_t>(&num_threads))(
           "refine", po::value<std::string>(&refine)->default_value("sparse"))(
           "print_centers", po::value<bool>(&print_centers))(
           "print_time_apply", po::value<bool>(&print_time_apply))(
@@ Expand Down Expand Up / @@ -118,6 +120,10 @@ int main(int argc, char* argv[]) { @@
       std::cout << adapt_opts << "\tsolve-rtol: " << solve_rtol << std::endl
                 << std::endl;
+      assert(num_threads > 0 && num_threads <= omp_get_max_threads() &&
+             num_threads <= MAX_NUMBER_THREADS);
+      omp_set_num_threads(num_threads);
       auto T = InitialTriangulation(domain, initial_refines);
       auto B = Time::Bases();
       auto vec_Xd = std::make_shared<
@@ Expand Down Expand Up / @@ -229,18 +235,19 @@ int main(int argc, char* argv[]) { @@
         if (print_time_apply) {
           auto heat_d_dd = heat_eq.heat_d_dd();
-          std::cout << "\n\tA-time-per-apply: " << heat_d_dd->A()->TimePerApply()
-                    << "\n\tB-time-per-apply: " << heat_d_dd->B()->TimePerApply()
-                    << "\n\tBT-time-per-apply: " << heat_d_dd->BT()->TimePerApply()
-                    << "\n\tG-time-per-apply: " << heat_d_dd->G()->TimePerApply()
-                    << "\n\tP_Y-time-per-apply: "
-                    << heat_d_dd->P_Y()->TimePerApply()
-                    << "\n\tP_X-time-per-apply: "
-                    << heat_d_dd->P_X()->TimePerApply()
-                    << "\n\tS-time-per-apply: " << heat_d_dd->S()->TimePerApply()
-                    << "\n\ttotal-time-apply: " << heat_d_dd->TotalTimeApply()
-                    << "\n\ttotal-time-construct: "
-                    << heat_d_dd->TotalTimeConstruct() << std::flush;
+          std::cout
+              << "\n\tA-time-per-apply: " << heat_d_dd->A()->TimePerApply()
+              << "\n\tB-time-per-apply: " << heat_d_dd->B()->TimePerApply()
+              << "\n\tB-A-time-per-apply: " << heat_d_dd->B()->A()->TimePerApply()
+              << "\n\tB-B-time-per-apply: " << heat_d_dd->B()->B()->TimePerApply()
+              << "\n\tBT-time-per-apply: " << heat_d_dd->BT()->TimePerApply()
+              << "\n\tG-time-per-apply: " << heat_d_dd->G()->TimePerApply()
+              << "\n\tP_Y-time-per-apply: " << heat_d_dd->P_Y()->TimePerApply()
+              << "\n\tP_X-time-per-apply: " << heat_d_dd->P_X()->TimePerApply()
+              << "\n\tS-time-per-apply: " << heat_d_dd->S()->TimePerApply()
+              << "\n\ttotal-time-apply: " << heat_d_dd->TotalTimeApply()
+              << "\n\ttotal-time-construct: " << heat_d_dd->TotalTimeConstruct()
+              << std::flush;
         }
         if (print_centers) {
@@ Expand Down @@

src/datastructures/boost.hpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -3,6 +3,7 @@ @@
     #include <boost/container/options.hpp>
     #include <boost/container/small_vector.hpp>
     #include <boost/container/static_vector.hpp>
+    #include <vector>
     template <typename I, size_t N>
     using SmallVector = boost::container::small_vector<I, N>;
@@ Expand Down @@

src/datastructures/multi_tree_view.ipp

-Original file line number
+Diff line change
@@ Expand Up / @@ -87,13 +87,14 @@ std::vector<I*> MultiNodeViewInterface<I, T...>::Union( @@
         // Now do the union magic in all dimensions.
         static_for<dim>([&queue, &my_node, &other_node, &call_filter](auto i) {
           // Get a list of all children of the other_node in axis `i`.
-          static std::vector<I_other*> filtered_children;
+          static thread_local std::vector<I_other*> filtered_children;
           filtered_children.clear();
           for (const auto& other_child_i : other_node->children(i))
             if (call_filter(other_child_i))
               filtered_children.emplace_back(other_child_i);
-          static std::vector<std::tuple_element_t<i, TupleNodes>> other_children_i;
+          static thread_local std::vector<std::tuple_element_t<i, TupleNodes>>
+              other_children_i;
           other_children_i.clear();
           for (const auto& other_child_i : filtered_children)
             other_children_i.emplace_back(std::get<i>(other_child_i->nodes()));
@@ Expand Down @@

src/datastructures/tree.hpp

            
                      Original file line number
                      Diff line number
                      Diff line change
                  
    @@ -1,5 +1,8 @@
  
    #pragma once

    #include <omp.h>

    #include <algorithm>

    #include <array>

    #include <memory>

    #include <queue>

    #include <utility>

    @@ -17,6 +20,12 @@ using T_func_noop = decltype(func_noop);
  
    using T_func_true = decltype(func_true);

    using T_func_false = decltype(func_false);

    // Global variable holding the current thread number.

    #ifndef MAX_NUMBER_THREADS

    #define MAX_NUMBER_THREADS 1

    #endif

    static thread_local int thread_number = omp_get_thread_num();

    template <typename I>

    struct NodeTrait;  // This should define N_children and N_parents.

    @@ -38,8 +47,8 @@ class Node {
  
      }

      int level() const { return level_; }

      bool marked() const { return marked_; }

      void set_marked(bool value) { marked_ = value; }

      bool marked() const { return marked_[thread_number]; }

      void set_marked(bool value) { marked_[thread_number] = value; }

      bool is_leaf() const { return children_.size() == 0; }

      inline bool is_metaroot() const { return (level_ == -1); }

      const auto &parents() const { return parents_; }

    @@ -48,20 +57,19 @@ class Node {
  
      // General data field for universal storage.

      template <typename T>

      T *data() {

        assert(data_ != nullptr);

        return static_cast<T *>(data_);

        assert(data_[thread_number] != nullptr);

        return static_cast<T *>(data_[thread_number]);

      }

      template <typename T>

      void set_data(T *value) {

        assert(data_ == nullptr);

        data_ = static_cast<void *>(value);

        assert(data_[thread_number] == nullptr);

        data_[thread_number] = static_cast<void *>(value);

      }

      void reset_data() {

        assert(data_ != nullptr);

        data_ = nullptr;

        assert(data_[thread_number] != nullptr);

        data_[thread_number] = nullptr;

      }

      bool has_data() { return data_ != nullptr; }

      bool has_data() { return data_[thread_number] != nullptr; }

      template <typename Func = T_func_noop>

      std::vector<I *> Bfs(bool include_metaroot = false,

    @@ -93,9 +101,9 @@ class Node {
  
      }

     protected:

      bool marked_ = false;

      int level_;

      void *data_ = nullptr;

      std::array<unsigned short, MAX_NUMBER_THREADS> marked_{0};

      std::array<void *, MAX_NUMBER_THREADS> data_{nullptr};

      // Store children/parents as raw pointers.

      SmallVector<I *, NodeTrait<I>::N_children> children_;

src/space/CMakeLists.txt

-Original file line number
+Diff line change
@@ -1,4 +1,5 @@
     add_library(space STATIC triangulation.cpp initial_triangulation.cpp basis.cpp operators.cpp triangulation_view.cpp linear_form.cpp integration.cpp)
+    target_link_libraries(space PUBLIC OpenMP::OpenMP_CXX)
     # Executables
     add_executable(space_adaptive adaptive.cpp)
@@ Expand Down @@

src/space/operators.cpp

            
                      Original file line number
                      Diff line number
                      Diff line change
                  
    @@ -213,13 +213,13 @@ void CGInverse<ForwardOp>::ApplySingleScale(Eigen::VectorXd &vec_SS) const {
  
    // Define the class variables.

    template <typename ForwardOp>

    std::vector<std::vector<std::pair<uint, double>>>

    thread_local std::vector<std::vector<std::pair<uint, double>>>

        MultigridPreconditioner<ForwardOp>::row_mat;

    template <typename ForwardOp>

    std::vector<std::vector<Element2D *>>

    thread_local std::vector<std::vector<Element2D *>>

        MultigridPreconditioner<ForwardOp>::patches;

    template <typename ForwardOp>

    std::vector<std::vector<uint>>

    thread_local std::vector<std::vector<uint>>

        MultigridPreconditioner<ForwardOp>::vertices_relaxation;

    template <typename ForwardOp>

    @@ -346,12 +346,12 @@ void MultigridPreconditioner<ForwardOp>::ApplySingleScale(
  
      // Shortcut.

      const uint V = triang_.V;

      // Reuse a static variable for storing the corrections.

      static std::vector<double> e;

      // Reuse a static variable for storing the row of a matrix.

      static thread_local std::vector<double> e;

      e.reserve(V * 3);

      // Reuse a static variable for storing the residual in the downard cycle.

      static std::vector<double> r_down;

      static thread_local std::vector<double> r_down;

      r_down.reserve(V * 3);

      // Initialize the multigrid matrix (row_mat).

src/space/operators.hpp

            
                      Original file line number
                      Diff line number
                      Diff line change
                  
    @@ -221,9 +221,9 @@ class MultigridPreconditioner : public BackwardOperator {
  
      DirectInverse<ForwardOp> initial_triang_solver_;

      // (Static) variables reused for calculation of the multigrid matrix.

      static std::vector<std::vector<std::pair<uint, double>>> row_mat;

      static std::vector<std::vector<Element2D *>> patches;

      static std::vector<std::vector<uint>> vertices_relaxation;

      static thread_local std::vector<std::vector<std::pair<uint, double>>> row_mat;

      static thread_local std::vector<std::vector<Element2D *>> patches;

      static thread_local std::vector<std::vector<uint>> vertices_relaxation;

    };

    template <template <typename> class InverseOp>

src/spacetime/bilinear_form.hpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -57,6 +57,8 @@ class BilinearForm @@
       auto sigma() { return sigma_; }
       auto theta() { return theta_; }
+      std::string Information() final;
      protected:
       // References to in/out vectors.
       DblVecIn *vec_in_;
@@ Expand All / @@ -72,6 +74,7 @@ class BilinearForm @@
       // Debug information.
       using BilinearFormBase<DblVecIn, DblVecOut>::time_construct_;
       using BilinearFormBase<DblVecIn, DblVecOut>::time_apply_;
+      using BilinearFormBase<DblVecIn, DblVecOut>::time_apply_split_;
       using BilinearFormBase<DblVecIn, DblVecOut>::num_apply_;
       // Define frozen templates, useful for storing the bil forms.
@@ Expand All / @@ -87,6 +90,14 @@ class BilinearForm @@
       std::vector<Time::BilinearForm<OperatorTime, FI<0>, FO<0>>> bil_time_low_;
       std::vector<Time::BilinearForm<OperatorTime, FI<0>, FO<0>>> bil_time_upp_;
       std::vector<space::BilinearForm<OperatorSpace, FO<1>, FO<1>>> bil_space_upp_;
+      // Store ordering for spatial parallism.
+      std::vector<FI<0> *> sigma_proj_0_;
+      std::vector<FO<1> *> theta_proj_1_;
+      std::vector<FO<0> *> vec_out_proj_0_;
+      std::vector<FO<1> *> vec_out_proj_1_;
+      std::vector<size_t> ordering_sigma_;
+      std::vector<size_t> ordering_vec_out_;
     };
     // Helper functions.
@@ Expand Down Expand Up / @@ -146,6 +157,7 @@ class BlockDiagonalBilinearForm @@
       Eigen::VectorXd Apply(const Eigen::VectorXd &v) final;
       DblVecIn *vec_in() const final { return vec_in_; }
       DblVecOut *vec_out() const final { return vec_out_; }
+      std::string Information() final;
      protected:
       bool use_cache_;
@@ Expand All / @@ -157,13 +169,17 @@ class BlockDiagonalBilinearForm @@
       // Debug information.
       using BilinearFormBase<DblVecIn, DblVecOut>::time_construct_;
       using BilinearFormBase<DblVecIn, DblVecOut>::time_apply_;
+      using BilinearFormBase<DblVecIn, DblVecOut>::time_apply_split_;
       using BilinearFormBase<DblVecIn, DblVecOut>::num_apply_;
       // The (cached) bilinear forms.
       template <size_t i>
       using FI = datastructures::FrozenDoubleNode<
           datastructures::DoubleNodeVector<BasisTimeIn, BasisSpace>, i>;
       std::vector<space::BilinearForm<OperatorSpace, FI<1>, FI<1>>> space_bilforms_;
+      std::vector<FI<0> *> vec_out_proj_0_;
+      std::vector<size_t> ordering_;
     };
     template <typename OpSpace, typename BTimeIn, typename BTimeOut>
@@ Expand Down @@

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

OpenMP implementation #210

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!

OpenMP implementation #210

Are you sure you want to change the base?

Uh oh!

OpenMP implementation #210

Uh oh!

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!