distributed_vector.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 #include "libmesh/libmesh_common.h"
21 
22 // C++ includes
23 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
24 #include <cmath> // for std::abs
25 #include <limits> // std::numeric_limits<T>::min()
26 
27 // Local Includes
29 #include "libmesh/dense_vector.h"
31 #include "libmesh/parallel.h"
32 #include "libmesh/parallel_sync.h"
33 #include "libmesh/tensor_tools.h"
34 
35 namespace libMesh
36 {
37 
38 
39 
40 //--------------------------------------------------------------------------
41 // DistributedVector methods
42 template <typename T>
44 {
45  // This function must be run on all processors at once
46  parallel_object_only();
47 
48  libmesh_assert (this->initialized());
49  libmesh_assert_equal_to (_values.size(), _local_size);
50  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
51 
52  T local_sum = 0.;
53 
54  for (numeric_index_type i=0; i<local_size(); i++)
55  local_sum += _values[i];
56 
57  this->comm().sum(local_sum);
58 
59  return local_sum;
60 }
61 
62 
63 
64 template <typename T>
66 {
67  // This function must be run on all processors at once
68  parallel_object_only();
69 
70  libmesh_assert (this->initialized());
71  libmesh_assert_equal_to (_values.size(), _local_size);
72  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
73 
74  double local_l1 = 0.;
75 
76  for (numeric_index_type i=0; i<local_size(); i++)
77  local_l1 += std::abs(_values[i]);
78 
79  this->comm().sum(local_l1);
80 
81  return local_l1;
82 }
83 
84 
85 
86 template <typename T>
88 {
89  // This function must be run on all processors at once
90  parallel_object_only();
91 
92  libmesh_assert (this->initialized());
93  libmesh_assert_equal_to (_values.size(), _local_size);
94  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
95 
96  double local_l2 = 0.;
97 
98  for (numeric_index_type i=0; i<local_size(); i++)
99  local_l2 += TensorTools::norm_sq(_values[i]);
100 
101  this->comm().sum(local_l2);
102 
103  return std::sqrt(local_l2);
104 }
105 
106 
107 
108 template <typename T>
110 {
111  // This function must be run on all processors at once
112  parallel_object_only();
113 
114  libmesh_assert (this->initialized());
115  libmesh_assert_equal_to (_values.size(), _local_size);
116  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
117 
118  Real local_linfty = 0.;
119 
120  for (numeric_index_type i=0; i<local_size(); i++)
121  local_linfty = std::max(local_linfty,
122  static_cast<Real>(std::abs(_values[i]))
123  ); // Note we static_cast so that both
124  // types are the same, as required
125  // by std::max
126 
127  this->comm().max(local_linfty);
128 
129  return local_linfty;
130 }
131 
132 
133 
134 template <typename T>
136 {
137  libmesh_assert (this->closed());
138  libmesh_assert (this->initialized());
139  libmesh_assert_equal_to (_values.size(), _local_size);
140  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
141 
142  add(1., v);
143 
144  return *this;
145 }
146 
147 
148 
149 template <typename T>
151 {
152  libmesh_assert (this->closed());
153  libmesh_assert (this->initialized());
154  libmesh_assert_equal_to (_values.size(), _local_size);
155  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
156 
157  add(-1., v);
158 
159  return *this;
160 }
161 
162 
163 
164 template <typename T>
166 {
167  libmesh_assert_equal_to(size(), v.size());
168 
169  const DistributedVector<T> & v_vec = cast_ref<const DistributedVector<T> &>(v);
170 
171  std::size_t val_size = _values.size();
172 
173  for (std::size_t i=0; i<val_size; i++)
174  _values[i] = _values[i] / v_vec._values[i];
175 
176  return *this;
177 }
178 
179 
180 
181 
182 template <typename T>
184 {
185  for (numeric_index_type i=0; i<local_size(); i++)
186  {
187  // Don't divide by zero
188  libmesh_assert_not_equal_to (_values[i], T(0));
189 
190  _values[i] = 1. / _values[i];
191  }
192 }
193 
194 
195 
196 
197 template <typename T>
199 {
200  for (numeric_index_type i=0; i<local_size(); i++)
201  {
202  // Replace values by complex conjugate
203  _values[i] = libmesh_conj(_values[i]);
204  }
205 }
206 
207 
208 
209 
210 
211 template <typename T>
213 {
214  libmesh_assert (this->initialized());
215  libmesh_assert_equal_to (_values.size(), _local_size);
216  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
217 
218  for (numeric_index_type i=0; i<local_size(); i++)
219  _values[i] += v;
220 }
221 
222 
223 
224 template <typename T>
226 {
227  libmesh_assert (this->initialized());
228  libmesh_assert_equal_to (_values.size(), _local_size);
229  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
230 
231  add (1., v);
232 }
233 
234 
235 
236 template <typename T>
237 void DistributedVector<T>::add (const T a, const NumericVector<T> & v)
238 {
239  libmesh_assert (this->initialized());
240  libmesh_assert_equal_to (_values.size(), _local_size);
241  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
242 
243  add(a, v);
244 }
245 
246 
247 
248 template <typename T>
249 void DistributedVector<T>::scale (const T factor)
250 {
251  libmesh_assert (this->initialized());
252  libmesh_assert_equal_to (_values.size(), _local_size);
253  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
254 
255  for (std::size_t i=0; i<local_size(); i++)
256  _values[i] *= factor;
257 }
258 
259 template <typename T>
261 {
262  libmesh_assert (this->initialized());
263  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
264 
265  for (numeric_index_type i=0; i<local_size(); i++)
266  this->set(i,std::abs(_values[i]));
267 }
268 
269 
270 
271 
272 
273 template <typename T>
275 {
276  // This function must be run on all processors at once
277  parallel_object_only();
278 
279  // Make sure the NumericVector passed in is really a DistributedVector
280  const DistributedVector<T> * v = cast_ptr<const DistributedVector<T> *>(&V);
281 
282  // Make sure that the two vectors are distributed in the same way.
283  libmesh_assert_equal_to ( this->first_local_index(), v->first_local_index() );
284  libmesh_assert_equal_to ( this->last_local_index(), v->last_local_index() );
285 
286  // The result of dotting together the local parts of the vector.
287  T local_dot = 0;
288 
289  for (std::size_t i=0; i<this->local_size(); i++)
290  local_dot += this->_values[i] * v->_values[i];
291 
292  // The local dot products are now summed via MPI
293  this->comm().sum(local_dot);
294 
295  return local_dot;
296 }
297 
298 
299 
300 template <typename T>
303 {
304  libmesh_assert (this->initialized());
305  libmesh_assert_equal_to (_values.size(), _local_size);
306  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
307 
308  for (std::size_t i=0; i<local_size(); i++)
309  _values[i] = s;
310 
311  return *this;
312 }
313 
314 
315 
316 template <typename T>
319 {
320  // Make sure the NumericVector passed in is really a DistributedVector
321  const DistributedVector<T> * v = cast_ptr<const DistributedVector<T> *>(&v_in);
322 
323  *this = *v;
324 
325  return *this;
326 }
327 
328 
329 
330 template <typename T>
333 {
335  this->_is_closed = v._is_closed;
336 
337  _global_size = v._global_size;
338  _local_size = v._local_size;
339  _first_local_index = v._first_local_index;
340  _last_local_index = v._last_local_index;
341 
342  if (v.local_size() == this->local_size())
343  _values = v._values;
344 
345  else
346  libmesh_error_msg("v.local_size() = " << v.local_size() << " must be equal to this->local_size() = " << this->local_size());
347 
348  return *this;
349 }
350 
351 
352 
353 template <typename T>
355 DistributedVector<T>::operator = (const std::vector<T> & v)
356 {
357  libmesh_assert (this->initialized());
358  libmesh_assert_equal_to (_values.size(), _local_size);
359  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
360 
361  if (v.size() == local_size())
362  _values = v;
363 
364  else if (v.size() == size())
365  for (std::size_t i=first_local_index(); i<last_local_index(); i++)
366  _values[i-first_local_index()] = v[i];
367 
368  else
369  libmesh_error_msg("Incompatible sizes in DistributedVector::operator=");
370 
371  return *this;
372 }
373 
374 
375 
376 template <typename T>
378 
379 {
380  libmesh_assert (this->initialized());
381  libmesh_assert_equal_to (_values.size(), _local_size);
382  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
383 
384  DistributedVector<T> * v_local = cast_ptr<DistributedVector<T> *>(&v_local_in);
385 
386  v_local->_first_local_index = 0;
387 
388  v_local->_global_size =
389  v_local->_local_size =
390  v_local->_last_local_index = size();
391 
392  v_local->_is_initialized =
393  v_local->_is_closed = true;
394 
395  // Call localize on the vector's values. This will help
396  // prevent code duplication
397  localize (v_local->_values);
398 
399 #ifndef LIBMESH_HAVE_MPI
400 
401  libmesh_assert_equal_to (local_size(), size());
402 
403 #endif
404 }
405 
406 
407 
408 template <typename T>
410  const std::vector<numeric_index_type> &) const
411 {
412  libmesh_assert (this->initialized());
413  libmesh_assert_equal_to (_values.size(), _local_size);
414  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
415 
416  // TODO: We don't yet support the send list; this is inefficient:
417  localize (v_local_in);
418 }
419 
420 
421 
422 template <typename T>
423 void DistributedVector<T>::localize (std::vector<T> & v_local,
424  const std::vector<numeric_index_type> & indices) const
425 {
426  // Resize v_local so there is enough room to hold all the local values.
427  v_local.resize(indices.size());
428 
429  // We need to know who has the values we want, so get everyone's _local_size
430  std::vector<numeric_index_type> local_sizes;
431  this->comm().allgather (_local_size, local_sizes);
432 
433  // Make a vector of partial sums of local sizes
434  std::vector<numeric_index_type> local_size_sums(this->n_processors());
435  local_size_sums[0] = local_sizes[0];
436  for (numeric_index_type i=1; i<local_sizes.size(); i++)
437  local_size_sums[i] = local_size_sums[i-1] + local_sizes[i];
438 
439  // We now fill in 'requested_ids' based on the indices. Also keep
440  // track of the local index (in the indices vector) for each of
441  // these, since we need that when unpacking.
442  std::map<processor_id_type, std::vector<numeric_index_type>>
443  requested_ids, local_requested_ids;
444 
445  // We'll use this typedef a couple of times below.
446  typedef typename std::vector<numeric_index_type>::iterator iter_t;
447 
448  // For each index in indices, determine which processor it is on.
449  // This is an O(N*log(p)) algorithm that uses std::upper_bound().
450  // Note: upper_bound() returns an iterator to the first entry which is
451  // greater than the given value.
452  for (numeric_index_type i=0; i<indices.size(); i++)
453  {
454  iter_t ub = std::upper_bound(local_size_sums.begin(),
455  local_size_sums.end(),
456  indices[i]);
457 
458  processor_id_type on_proc = cast_int<processor_id_type>
459  (std::distance(local_size_sums.begin(), ub));
460 
461  requested_ids[on_proc].push_back(indices[i]);
462  local_requested_ids[on_proc].push_back(i);
463  }
464 
465  auto gather_functor =
466  [this]
467  (processor_id_type, const std::vector<dof_id_type> & ids,
468  std::vector<T> & values)
469  {
470  // The first send/receive we did was for indices, the second one will be
471  // for corresponding floating point values, so create storage for that now...
472  const std::size_t ids_size = ids.size();
473  values.resize(ids_size);
474 
475  for (std::size_t i=0; i != ids_size; i++)
476  {
477  // The index of the requested value
478  const numeric_index_type requested_index = ids[i];
479 
480  // Transform into local numbering, and get requested value.
481  values[i] = _values[requested_index - _first_local_index];
482  }
483  };
484 
485  auto action_functor =
486  [& v_local, & local_requested_ids]
487  (processor_id_type pid,
488  const std::vector<dof_id_type> &,
489  const std::vector<T> & values)
490  {
491  // Now write the received values to the appropriate place(s) in v_local
492  for (std::size_t i=0, vsize = values.size(); i != vsize; i++)
493  {
494  libmesh_assert(local_requested_ids.count(pid));
495  libmesh_assert_less(i, local_requested_ids[pid].size());
496 
497  // Get the index in v_local where this value needs to be inserted.
498  const numeric_index_type local_requested_index =
499  local_requested_ids[pid][i];
500 
501  // Actually set the value in v_local
502  v_local[local_requested_index] = values[i];
503  }
504  };
505 
506  const T * ex = nullptr;
508  (this->comm(), requested_ids, gather_functor, action_functor, ex);
509 }
510 
511 
512 
513 template <typename T>
515  const numeric_index_type last_local_idx,
516  const std::vector<numeric_index_type> & send_list)
517 {
518  // Only good for serial vectors
519  libmesh_assert_equal_to (this->size(), this->local_size());
520  libmesh_assert_greater (last_local_idx, first_local_idx);
521  libmesh_assert_less_equal (send_list.size(), this->size());
522  libmesh_assert_less (last_local_idx, this->size());
523 
524  const numeric_index_type my_size = this->size();
525  const numeric_index_type my_local_size = (last_local_idx - first_local_idx + 1);
526 
527  // Don't bother for serial cases
528  if ((first_local_idx == 0) &&
529  (my_local_size == my_size))
530  return;
531 
532 
533  // Build a parallel vector, initialize it with the local
534  // parts of (*this)
535  DistributedVector<T> parallel_vec(this->comm());
536 
537  parallel_vec.init (my_size, my_local_size, true, PARALLEL);
538 
539  // Copy part of *this into the parallel_vec
540  for (numeric_index_type i=first_local_idx; i<=last_local_idx; i++)
541  parallel_vec._values[i-first_local_idx] = _values[i];
542 
543  // localize like normal
544  parallel_vec.localize (*this, send_list);
545 }
546 
547 
548 
549 template <typename T>
550 void DistributedVector<T>::localize (std::vector<T> & v_local) const
551 {
552  // This function must be run on all processors at once
553  parallel_object_only();
554 
555  libmesh_assert (this->initialized());
556  libmesh_assert_equal_to (_values.size(), _local_size);
557  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
558 
559  v_local = this->_values;
560 
561  this->comm().allgather (v_local);
562 
563 #ifndef LIBMESH_HAVE_MPI
564  libmesh_assert_equal_to (local_size(), size());
565 #endif
566 }
567 
568 
569 
570 template <typename T>
571 void DistributedVector<T>::localize_to_one (std::vector<T> & v_local,
572  const processor_id_type pid) const
573 {
574  // This function must be run on all processors at once
575  parallel_object_only();
576 
577  libmesh_assert (this->initialized());
578  libmesh_assert_equal_to (_values.size(), _local_size);
579  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
580 
581  v_local = this->_values;
582 
583  this->comm().gather (pid, v_local);
584 
585 #ifndef LIBMESH_HAVE_MPI
586  libmesh_assert_equal_to (local_size(), size());
587 #endif
588 }
589 
590 
591 
592 template <typename T>
594  const NumericVector<T> &)
595 //void DistributedVector<T>::pointwise_mult (const NumericVector<T> & vec1,
596 // const NumericVector<T> & vec2)
597 {
598  libmesh_not_implemented();
599 }
600 
601 
602 
603 //--------------------------------------------------------------
604 // Explicit instantiations
605 template class DistributedVector<Number>;
606 
607 } // namespace libMesh
bool closed()
Definition: libmesh.C:265
double abs(double a)
T libmesh_conj(T a)
virtual T dot(const NumericVector< T > &V) const override
numeric_index_type _last_local_index
virtual void scale(const T factor) override
virtual numeric_index_type last_local_index() const override
virtual numeric_index_type size() const =0
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) override
DistributedVector & operator=(const DistributedVector &)
uint8_t processor_id_type
Definition: id_types.h:99
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
virtual NumericVector< T > & operator/=(const NumericVector< T > &v) override
virtual void add(const numeric_index_type i, const T value) override
virtual Real linfty_norm() const override
numeric_index_type _global_size
long double max(long double a, double b)
virtual numeric_index_type local_size() const override
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) override
virtual void reciprocal() override
virtual Real l1_norm() const override
dof_id_type numeric_index_type
Definition: id_types.h:92
void pull_parallel_vector_data(const Communicator &comm, const MapToVectors &queries, RequestContainer &reqs, GatherFunctor &gather_data, ActionFunctor &act_on_data, const datum *example)
virtual numeric_index_type size() const override
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
virtual void conjugate() override
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const override
numeric_index_type _local_size
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual numeric_index_type first_local_index() const override
virtual void localize(std::vector< T > &v_local) const override
virtual Real l2_norm() const override
virtual void abs() override
bool initialized()
Definition: libmesh.C:258
virtual void init(const numeric_index_type N, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC) override
virtual T sum() const override
numeric_index_type _first_local_index