GetFEM  5.4.2
getfem_accumulated_distro.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2018 Andriy Andreykiv
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file getfem_accumulated_distro.h
33 @author Andriy Andreykiv <andriy.andreykiv@gmail.com>
34 @date November 21, 2018.
35 @brief Distribution of assembly results (matrices/vectors) for parallel
36  assembly.
37 
38 This is the kernel of getfem.
39 */
40 #pragma once
41 
42 #include <gmm/gmm_def.h>
43 
44 #include "getfem_omp.h"
45 
46 namespace getfem
47 {
48 
49 namespace detail {
50 
51  //T is a single vector
52  template<class T>
53  void add_spec(const T &a, T &b,
54  gmm::abstract_null_type, gmm::abstract_vector){
55  gmm::add(a, b);
56  }
57 
58  //T is a single matrix
59  template<class T>
60  void add_spec(const T &a, T &b,
61  gmm::abstract_null_type, gmm::abstract_matrix){
62  gmm::add(a, b);
63  }
64 
65  //T is a vector of either matrices or vectors
66  template<class T>
67  void add_list(const T &a, T &b){
68  GMM_ASSERT2(a.size() == b.size(), "size mismatch");
69  auto ita = begin(a);
70  auto itb = begin(b);
71  auto ita_end = end(a);
72  for (;ita != ita_end; ++ita, ++itb) gmm::add(*ita, *itb);
73  }
74 
75  //T is a vector of vectors
76  template<class T>
77  void add_spec(const T &a, T &b,
78  gmm::abstract_vector, gmm::abstract_vector){
79  add_list(a, b);
80  }
81 
82  //T is a vector of matrices
83  template<class T>
84  void add_spec(const T &a, T &b,
85  gmm::abstract_matrix, gmm::abstract_vector){
86  add_list(a, b);
87  }
88 
89  //T is a single vector
90  template<class T>
91  void equal_resize_spec(T &a, const T &b,
92  gmm::abstract_null_type, gmm::abstract_vector){
93  gmm::resize(a, gmm::vect_size(b));
94  }
95 
96  //T is a single matrix
97  template<class T>
98  void equal_resize_spec(T &a, const T &b,
99  gmm::abstract_null_type, gmm::abstract_matrix){
100  gmm::resize(a, gmm::mat_nrows(b), gmm::mat_ncols(b));
101  }
102 
103  //T is a vector of either matrices or vectors
104  template<class T>
105  void equal_resize_list(T &a, const T &b){
106  GMM_ASSERT2(a.empty(), "the first list should be still empty");
107  a.resize(b.size());
108  auto ita = begin(a);
109  auto itb = begin(b);
110  auto ita_end = end(a);
111  using Component = typename T::value_type;
112  using AlgoC = typename gmm::linalg_traits<Component>::linalg_type;
113  for (;ita != ita_end; ++ita, ++itb){
114  equal_resize_spec(*ita, *itb, gmm::abstract_null_type{}, AlgoC{});
115  }
116  }
117 
118  //T is a vector of vectors
119  template<class T>
120  void equal_resize_spec(T &a, const T &b,
121  gmm::abstract_vector, gmm::abstract_vector){
122  equal_resize_list(a, b);
123  }
124 
125  //T is a vector of matrices
126  template<class T>
127  void equal_resize_spec(T &a, const T &b,
128  gmm::abstract_matrix, gmm::abstract_vector){
129  equal_resize_list(a, b);
130  }
131 
132 } // namespace detail
133 
134  /**Generic addition for gmm types as well as vectors of gmm types*/
135  template<class T>
136  void gen_add(const T &a, T &b){
137  using AlgoT = typename gmm::linalg_traits<T>::linalg_type;
138  using Component = typename T::value_type;
139  using AlgoC = typename gmm::linalg_traits<Component>::linalg_type;
140  detail::add_spec(a, b, AlgoC{}, AlgoT{});
141  }
142 
143  /**Resize 'a' to the same size as 'b'.
144  a and b can be matrices, vectors, or lists of matrices or vectors
145  */
146  template<class T>
147  void equal_resize(T &a, const T &b){
148  using AlgoT = typename gmm::linalg_traits<T>::linalg_type;
149  using Component = typename T::value_type;
150  using AlgoC = typename gmm::linalg_traits<Component>::linalg_type;
151  detail::equal_resize_spec(a, b, AlgoC{}, AlgoT{});
152  }
153 
154  /**Takes a matrix or vector, or vector of matrices or vectors and
155  creates an empty copy on each thread. When the
156  thread computations are done (in the destructor), accumulates
157  the assembled copies into the original*/
158  template <class T>
160  {
161  T& original;
162  omp_distribute<T> distributed;
163 
164  public:
165 
166  explicit accumulated_distro(T& l)
167  : original{l}{
168  if (distributed.num_threads() == 1) return;
169  //intentionally skipping thread 0, as accumulated_distro will
170  //use original for it
171  for(size_type t = 1; t != distributed.num_threads(); ++t){
172  equal_resize(distributed(t), original);
173  }
174  }
175 
176  T& get(){
177  if (distributed.num_threads() == 1 ||
178  distributed.this_thread() == 0) return original;
179  else return distributed;
180  }
181 
182  operator T&(){ // implicit conversion
183  return get();
184  }
185 
186  T& operator = (const T &x){
187  return distributed = x;
188  }
189 
191  if (distributed.num_threads() == 1) return;
192 
193  if (me_is_multithreaded_now()) {
194  // GMM_ASSERT1 not convenient here
195  cerr << "Accumulation distribution should not run in parallel";
196  exit(1);
197  }
198 
199  using namespace std;
200  auto to_add = vector<T*>{};
201  to_add.push_back(&original);
202  for (size_type t = 1; t != distributed.num_threads(); ++t){
203  to_add.push_back(&distributed(t));
204  }
205 
206  //Accumulation in parallel.
207  //Adding, for instance, elements 1 to 0, 2 to 3, 5 to 4 and 7 to 6
208  //on separate 4 threads in case of parallelization of the assembly
209  //on 8 threads.
210  while (to_add.size() > 1){
212  auto i = distributed.this_thread() * 2;
213  if (i + 1 < to_add.size()){
214  auto &target = *to_add[i];
215  auto &source = *to_add[i + 1];
216  gen_add(source, target);
217  }
218  )
219  //erase every second item , as it was already added
220  for (auto it = begin(to_add), ite = end(to_add);
221  it != end(to_add) && next(it) != end(to_add);
222  it = to_add.erase(next(it)));
223  }
224  }
225  };
226 
227 } // namespace getfem
getfem::me_is_multithreaded_now
bool me_is_multithreaded_now()
is the program running in the parallel section
Definition: getfem_omp.cc:106
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm_def.h
Basic definitions and tools of GMM.
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::gen_add
void gen_add(const T &a, T &b)
Generic addition for gmm types as well as vectors of gmm types.
Definition: getfem_accumulated_distro.h:136
getfem_omp.h
Tools for multithreaded, OpenMP and Boost based parallelization.
getfem::omp_distribute< T >
GETFEM_OMP_PARALLEL
#define GETFEM_OMP_PARALLEL(body)
Organizes a proper parallel omp section:
Definition: getfem_omp.h:483
getfem::equal_resize
void equal_resize(T &a, const T &b)
Resize 'a' to the same size as 'b'.
Definition: getfem_accumulated_distro.h:147
getfem::accumulated_distro
Takes a matrix or vector, or vector of matrices or vectors and creates an empty copy on each thread.
Definition: getfem_accumulated_distro.h:159