distributed_vector.h
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 
23 
24 #ifndef LIBMESH_DISTRIBUTED_VECTOR_H
25 #define LIBMESH_DISTRIBUTED_VECTOR_H
26 
27 // Local includes
28 #include "libmesh/numeric_vector.h"
29 #include "libmesh/parallel.h"
30 
31 // C++ includes
32 #include <vector>
33 #include <algorithm>
34 #include <limits>
35 
36 namespace libMesh
37 {
38 
51 template <typename T>
52 class DistributedVector final : public NumericVector<T>
53 {
54 public:
55 
59  explicit
61  const ParallelType = AUTOMATIC);
62 
66  explicit
68  const numeric_index_type n,
69  const ParallelType ptype = AUTOMATIC);
70 
76  const numeric_index_type n,
77  const numeric_index_type n_local,
78  const ParallelType ptype = AUTOMATIC);
79 
86  const numeric_index_type N,
87  const numeric_index_type n_local,
88  const std::vector<numeric_index_type> & ghost,
89  const ParallelType ptype = AUTOMATIC);
90 
99 
104  DistributedVector (DistributedVector &&) = default;
105  DistributedVector (const DistributedVector &) = default;
107  virtual ~DistributedVector () = default;
108 
109  virtual void close () override;
110 
111  virtual void clear () override;
112 
113  virtual void zero () override;
114 
115  virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
116 
117  virtual std::unique_ptr<NumericVector<T>> clone () const override;
118 
119  virtual void init (const numeric_index_type N,
120  const numeric_index_type n_local,
121  const bool fast=false,
122  const ParallelType ptype=AUTOMATIC) override;
123 
124  virtual void init (const numeric_index_type N,
125  const bool fast=false,
126  const ParallelType ptype=AUTOMATIC) override;
127 
128  virtual void init (const numeric_index_type N,
129  const numeric_index_type n_local,
130  const std::vector<numeric_index_type> & ghost,
131  const bool fast = false,
132  const ParallelType = AUTOMATIC) override;
133 
134  virtual void init (const NumericVector<T> & other,
135  const bool fast = false) override;
136 
137  virtual NumericVector<T> & operator= (const T s) override;
138 
139  virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
140 
141  virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
142 
143  virtual Real min () const override;
144 
145  virtual Real max () const override;
146 
147  virtual T sum() const override;
148 
149  virtual Real l1_norm () const override;
150 
151  virtual Real l2_norm () const override;
152 
153  virtual Real linfty_norm () const override;
154 
155  virtual numeric_index_type size () const override;
156 
157  virtual numeric_index_type local_size() const override;
158 
159  virtual numeric_index_type first_local_index() const override;
160 
161  virtual numeric_index_type last_local_index() const override;
162 
163  virtual T operator() (const numeric_index_type i) const override;
164 
165  virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
166 
167  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
168 
169  virtual NumericVector<T> & operator /= (const NumericVector<T> & v) override;
170 
171  virtual void reciprocal() override;
172 
173  virtual void conjugate() override;
174 
175  virtual void set (const numeric_index_type i, const T value) override;
176 
177  virtual void add (const numeric_index_type i, const T value) override;
178 
179  virtual void add (const T s) override;
180 
181  virtual void add (const NumericVector<T> & V) override;
182 
183  virtual void add (const T a, const NumericVector<T> & v) override;
184 
190 
191  virtual void add_vector (const NumericVector<T> &,
192  const SparseMatrix<T> &) override
193  { libmesh_not_implemented(); }
194 
195  virtual void add_vector_transpose (const NumericVector<T> &,
196  const SparseMatrix<T> &) override
197  { libmesh_not_implemented(); }
198 
199  virtual void scale (const T factor) override;
200 
201  virtual void abs() override;
202 
203  virtual T dot(const NumericVector<T> & V) const override;
204 
205  virtual void localize (std::vector<T> & v_local) const override;
206 
207  virtual void localize (NumericVector<T> & v_local) const override;
208 
209  virtual void localize (NumericVector<T> & v_local,
210  const std::vector<numeric_index_type> & send_list) const override;
211 
212  virtual void localize (std::vector<T> & v_local,
213  const std::vector<numeric_index_type> & indices) const override;
214 
215  virtual void localize (const numeric_index_type first_local_idx,
216  const numeric_index_type last_local_idx,
217  const std::vector<numeric_index_type> & send_list) override;
218 
219  virtual void localize_to_one (std::vector<T> & v_local,
220  const processor_id_type proc_id=0) const override;
221 
222  virtual void pointwise_mult (const NumericVector<T> & vec1,
223  const NumericVector<T> & vec2) override;
224 
225  virtual void swap (NumericVector<T> & v) override;
226 
227 private:
228 
232  std::vector<T> _values;
233 
238 
243 
248 
253 };
254 
255 
256 //--------------------------------------------------------------------------
257 // DistributedVector inline methods
258 template <typename T>
259 inline
261  const ParallelType ptype) :
262  NumericVector<T>(comm_in, ptype),
263  _global_size (0),
264  _local_size (0),
265  _first_local_index(0),
266  _last_local_index (0)
267 {
268  this->_type = ptype;
269 }
270 
271 
272 
273 template <typename T>
274 inline
276  const numeric_index_type n,
277  const ParallelType ptype)
278  : NumericVector<T>(comm_in, ptype)
279 {
280  this->init(n, n, false, ptype);
281 }
282 
283 
284 
285 template <typename T>
286 inline
288  const numeric_index_type n,
289  const numeric_index_type n_local,
290  const ParallelType ptype)
291  : NumericVector<T>(comm_in, ptype)
292 {
293  this->init(n, n_local, false, ptype);
294 }
295 
296 
297 
298 template <typename T>
299 inline
301  const numeric_index_type n,
302  const numeric_index_type n_local,
303  const std::vector<numeric_index_type> & ghost,
304  const ParallelType ptype)
305  : NumericVector<T>(comm_in, ptype)
306 {
307  this->init(n, n_local, ghost, false, ptype);
308 }
309 
310 
311 
312 template <typename T>
313 inline
315  const numeric_index_type n_local,
316  const bool fast,
317  const ParallelType ptype)
318 {
319  // This function must be run on all processors at once
320  parallel_object_only();
321 
322  libmesh_assert_less_equal (n_local, n);
323 
324  if (ptype == AUTOMATIC)
325  {
326  if (n == n_local)
327  this->_type = SERIAL;
328  else
329  this->_type = PARALLEL;
330  }
331  else
332  this->_type = ptype;
333 
334  libmesh_assert ((this->_type==SERIAL && n==n_local) ||
335  this->_type==PARALLEL);
336 
337  // Clear the data structures if already initialized
338  if (this->initialized())
339  this->clear();
340 
341  // Initialize data structures
342  _values.resize(n_local);
343  _local_size = n_local;
344  _global_size = n;
345 
346  _first_local_index = 0;
347 
348 #ifdef LIBMESH_HAVE_MPI
349 
350  std::vector<numeric_index_type> local_sizes (this->n_processors(), 0);
351 
352  local_sizes[this->processor_id()] = n_local;
353 
354  this->comm().sum(local_sizes);
355 
356  // _first_local_index is the sum of _local_size
357  // for all processor ids less than ours
358  for (processor_id_type p=0; p!=this->processor_id(); p++)
359  _first_local_index += local_sizes[p];
360 
361 
362 # ifdef DEBUG
363  // Make sure all the local sizes sum up to the global
364  // size, otherwise there is big trouble!
365  numeric_index_type dbg_sum=0;
366 
367  for (processor_id_type p=0; p!=this->n_processors(); p++)
368  dbg_sum += local_sizes[p];
369 
370  libmesh_assert_equal_to (dbg_sum, n);
371 
372 # endif
373 
374 #else
375 
376  // No other options without MPI!
377  if (n != n_local)
378  libmesh_error_msg("ERROR: MPI is required for n != n_local!");
379 
380 #endif
381 
382  _last_local_index = _first_local_index + n_local;
383 
384  // Set the initialized flag
385  this->_is_initialized = true;
386 
387  // Zero the components unless directed otherwise
388  if (!fast)
389  this->zero();
390 }
391 
392 
393 template <typename T>
394 inline
396  const numeric_index_type n_local,
397  const std::vector<numeric_index_type> & /*ghost*/,
398  const bool fast,
399  const ParallelType ptype)
400 {
401  // TODO: we shouldn't ignore the ghost sparsity pattern
402  this->init(n, n_local, fast, ptype);
403 }
404 
405 
406 
407 /* Default implementation for solver packages for which ghosted
408  vectors are not yet implemented. */
409 template <class T>
411  const bool fast)
412 {
413  this->init(other.size(),other.local_size(),fast,other.type());
414 }
415 
416 
417 
418 template <typename T>
419 inline
421  const bool fast,
422  const ParallelType ptype)
423 {
424  this->init(n,n,fast,ptype);
425 }
426 
427 
428 
429 template <typename T>
430 inline
432 {
433  libmesh_assert (this->initialized());
434 
435  this->_is_closed = true;
436 }
437 
438 
439 
440 template <typename T>
441 inline
443 {
444  _values.clear();
445 
446  _global_size =
447  _local_size =
448  _first_local_index =
449  _last_local_index = 0;
450 
451 
452  this->_is_closed = this->_is_initialized = false;
453 }
454 
455 
456 
457 template <typename T>
458 inline
460 {
461  libmesh_assert (this->initialized());
462  libmesh_assert_equal_to (_values.size(), _local_size);
463  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
464 
465  std::fill (_values.begin(),
466  _values.end(),
467  0.);
468 }
469 
470 
471 
472 template <typename T>
473 inline
474 std::unique_ptr<NumericVector<T>> DistributedVector<T>::zero_clone () const
475 {
476  NumericVector<T> * cloned_vector = new DistributedVector<T>(this->comm());
477  cloned_vector->init(*this);
478  return std::unique_ptr<NumericVector<T>>(cloned_vector);
479 }
480 
481 
482 
483 template <typename T>
484 inline
485 std::unique_ptr<NumericVector<T>> DistributedVector<T>::clone () const
486 {
487  NumericVector<T> * cloned_vector = new DistributedVector<T>(this->comm());
488  cloned_vector->init(*this, true);
489  *cloned_vector = *this;
490  return std::unique_ptr<NumericVector<T>>(cloned_vector);
491 }
492 
493 
494 
495 template <typename T>
496 inline
498 {
499  libmesh_assert (this->initialized());
500  libmesh_assert_equal_to (_values.size(), _local_size);
501  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
502 
503  return _global_size;
504 }
505 
506 
507 
508 template <typename T>
509 inline
511 {
512  libmesh_assert (this->initialized());
513  libmesh_assert_equal_to (_values.size(), _local_size);
514  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
515 
516  return _local_size;
517 }
518 
519 
520 
521 template <typename T>
522 inline
524 {
525  libmesh_assert (this->initialized());
526  libmesh_assert_equal_to (_values.size(), _local_size);
527  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
528 
529  return _first_local_index;
530 }
531 
532 
533 
534 template <typename T>
535 inline
537 {
538  libmesh_assert (this->initialized());
539  libmesh_assert_equal_to (_values.size(), _local_size);
540  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
541 
542  return _last_local_index;
543 }
544 
545 
546 
547 template <typename T>
548 inline
550 {
551  libmesh_assert (this->initialized());
552  libmesh_assert_equal_to (_values.size(), _local_size);
553  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
554  libmesh_assert ( ((i >= first_local_index()) &&
555  (i < last_local_index())) );
556 
557  return _values[i - _first_local_index];
558 }
559 
560 
561 
562 template <typename T>
563 inline
565 {
566  libmesh_assert (this->initialized());
567  libmesh_assert_equal_to (_values.size(), _local_size);
568  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
569  libmesh_assert_less (i, size());
570  libmesh_assert_less (i-first_local_index(), local_size());
571 
572  _values[i - _first_local_index] = value;
573 }
574 
575 
576 
577 template <typename T>
578 inline
580 {
581  libmesh_assert (this->initialized());
582  libmesh_assert_equal_to (_values.size(), _local_size);
583  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
584  libmesh_assert_less (i, size());
585  libmesh_assert_less (i-first_local_index(), local_size());
586 
587  _values[i - _first_local_index] += value;
588 }
589 
590 
591 
592 template <typename T>
593 inline
595 {
596  // This function must be run on all processors at once
597  parallel_object_only();
598 
599  libmesh_assert (this->initialized());
600  libmesh_assert_equal_to (_values.size(), _local_size);
601  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
602 
603  Real local_min = _values.size() ?
605  for (numeric_index_type i = 1; i < _values.size(); ++i)
606  local_min = std::min(libmesh_real(_values[i]), local_min);
607 
608  this->comm().min(local_min);
609 
610  return local_min;
611 }
612 
613 
614 
615 template <typename T>
616 inline
618 {
619  // This function must be run on all processors at once
620  parallel_object_only();
621 
622  libmesh_assert (this->initialized());
623  libmesh_assert_equal_to (_values.size(), _local_size);
624  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
625 
626  Real local_max = _values.size() ?
628  for (numeric_index_type i = 1; i < _values.size(); ++i)
629  local_max = std::max(libmesh_real(_values[i]), local_max);
630 
631  this->comm().max(local_max);
632 
633  return local_max;
634 }
635 
636 
637 template <typename T>
638 inline
640 {
641  DistributedVector<T> & v = cast_ref<DistributedVector<T> &>(other);
642 
643  std::swap(_global_size, v._global_size);
644  std::swap(_local_size, v._local_size);
645  std::swap(_first_local_index, v._first_local_index);
646  std::swap(_last_local_index, v._last_local_index);
647 
648  // This should be O(1) with any reasonable STL implementation
649  std::swap(_values, v._values);
650 }
651 
652 } // namespace libMesh
653 
654 
655 #endif // LIBMESH_DISTRIBUTED_VECTOR_H
T libmesh_real(T a)
virtual Real max() const override
virtual T dot(const NumericVector< T > &V) const override
numeric_index_type _last_local_index
virtual void scale(const T factor) override
virtual void add_vector(const NumericVector< T > &, const SparseMatrix< T > &) override
virtual numeric_index_type last_local_index() const override
virtual void close() override
virtual numeric_index_type size() const =0
virtual void set(const numeric_index_type i, const T value) override
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
const Parallel::Communicator & comm() const
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
virtual std::unique_ptr< NumericVector< T > > zero_clone() const override
const Number zero
Definition: libmesh.h:239
long double max(long double a, double b)
virtual numeric_index_type local_size() const override
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) override
virtual void reciprocal() override
virtual void zero() override
virtual Real l1_norm() const override
dof_id_type numeric_index_type
Definition: id_types.h:92
virtual ~DistributedVector()=default
virtual numeric_index_type size() const override
void init(triangulateio &t)
virtual std::unique_ptr< NumericVector< T > > clone() const override
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
virtual void conjugate() override
ParallelType type() const
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
void swap(Iterator &lhs, Iterator &rhs)
virtual Real l2_norm() const override
virtual numeric_index_type local_size() const =0
virtual void abs() override
static const bool value
Definition: xdr_io.C:109
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 void clear() override
virtual Real min() const override
long double min(long double a, double b)
virtual T sum() const override
virtual void swap(NumericVector< T > &v) override
virtual void add_vector_transpose(const NumericVector< T > &, const SparseMatrix< T > &) override
DistributedVector(const Parallel::Communicator &comm, const ParallelType=AUTOMATIC)
numeric_index_type _first_local_index
virtual T operator()(const numeric_index_type i) const override