Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
read_write_vector.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_read_write_vector_h
17 #define dealii_read_write_vector_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/index_set.h>
22 #include <deal.II/base/memory_consumption.h>
23 #include <deal.II/base/mpi.h>
24 #include <deal.II/base/template_constraints.h>
25 #include <deal.II/base/thread_management.h>
26 #include <deal.II/base/types.h>
27 #include <deal.II/base/utilities.h>
28 
29 #include <deal.II/lac/vector_operation.h>
30 
31 #include <cstdlib>
32 #include <cstring>
33 #include <iomanip>
34 
35 #ifdef DEAL_II_WITH_TRILINOS
36 # include <deal.II/lac/trilinos_epetra_communication_pattern.h>
37 # include <deal.II/lac/trilinos_epetra_vector.h>
38 # include <deal.II/lac/trilinos_tpetra_vector.h>
39 
40 # include <Epetra_MultiVector.h>
41 
42 #endif
43 
44 DEAL_II_NAMESPACE_OPEN
45 
46 namespace LinearAlgebra
47 {
48  class CommunicationPatternBase;
49  namespace distributed
50  {
51  template <typename, typename>
52  class Vector;
53  } // namespace distributed
54 } // namespace LinearAlgebra
55 
56 #ifdef DEAL_II_WITH_PETSC
57 namespace PETScWrappers
58 {
59  namespace MPI
60  {
61  class Vector;
62  }
63 } // namespace PETScWrappers
64 #endif
65 
66 #ifdef DEAL_II_WITH_TRILINOS
67 namespace TrilinosWrappers
68 {
69  namespace MPI
70  {
71  class Vector;
72  }
73 } // namespace TrilinosWrappers
74 #endif
75 
76 #ifdef DEAL_II_WITH_CUDA
77 namespace LinearAlgebra
78 {
79  namespace CUDAWrappers
80  {
81  template <typename>
82  class Vector;
83  }
84 } // namespace LinearAlgebra
85 #endif
86 
87 namespace LinearAlgebra
88 {
128  template <typename Number>
129  class ReadWriteVector : public Subscriptor
130  {
131  public:
137  using value_type = Number;
138  using pointer = value_type *;
139  using const_pointer = const value_type *;
140  using iterator = value_type *;
141  using const_iterator = const value_type *;
142  using reference = value_type &;
143  using const_reference = const value_type &;
144  using size_type = types::global_dof_index;
145  using real_type = typename numbers::NumberTraits<Number>::real_type;
146 
154  ReadWriteVector();
155 
159  ReadWriteVector(const ReadWriteVector<Number> &in_vector);
160 
165  explicit ReadWriteVector(const size_type size);
166 
171  explicit ReadWriteVector(const IndexSet &locally_stored_indices);
172 
176  ~ReadWriteVector() override = default;
177 
186  virtual void
187  reinit(const size_type size, const bool omit_zeroing_entries = false);
188 
197  template <typename Number2>
198  void
199  reinit(const ReadWriteVector<Number2> &in_vector,
200  const bool omit_zeroing_entries = false);
201 
211  virtual void
212  reinit(const IndexSet &locally_stored_indices,
213  const bool omit_zeroing_entries = false);
214 
215 
216 #ifdef DEAL_II_WITH_TRILINOS
217 # ifdef DEAL_II_WITH_MPI
218 
229  void
230  reinit(const TrilinosWrappers::MPI::Vector &trilinos_vec);
231 # endif
232 #endif
233 
247  template <typename Functor>
248  void
249  apply(const Functor &func);
250 
263  void
265 
270  operator=(const ReadWriteVector<Number> &in_vector);
271 
275  template <typename Number2>
277  operator=(const ReadWriteVector<Number2> &in_vector);
278 
284  operator=(const Number s);
285 
294  template <typename MemorySpace>
295  void
297  VectorOperation::values operation,
298  const std::shared_ptr<const CommunicationPatternBase>
299  &communication_pattern =
300  std::shared_ptr<const CommunicationPatternBase>());
301 
302 #ifdef DEAL_II_WITH_PETSC
303 
311  void
312  import(const PETScWrappers::MPI::Vector &petsc_vec,
313  VectorOperation::values operation,
314  const std::shared_ptr<const CommunicationPatternBase>
315  &communication_pattern =
316  std::shared_ptr<const CommunicationPatternBase>());
317 #endif
318 
319 #ifdef DEAL_II_WITH_TRILINOS
320 
330  void
331  import(const TrilinosWrappers::MPI::Vector &trilinos_vec,
332  VectorOperation::values operation,
333  const std::shared_ptr<const CommunicationPatternBase>
334  &communication_pattern =
335  std::shared_ptr<const CommunicationPatternBase>());
336 
337 # ifdef DEAL_II_WITH_MPI
338 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
339 
347  void
348  import(const TpetraWrappers::Vector<Number> &tpetra_vec,
349  VectorOperation::values operation,
350  const std::shared_ptr<const CommunicationPatternBase>
351  &communication_pattern =
352  std::shared_ptr<const CommunicationPatternBase>());
353 # endif
354 
363  void
364  import(const EpetraWrappers::Vector &epetra_vec,
365  VectorOperation::values operation,
366  const std::shared_ptr<const CommunicationPatternBase>
367  &communication_pattern =
368  std::shared_ptr<const CommunicationPatternBase>());
369 # endif
370 #endif
371 
372 #ifdef DEAL_II_WITH_CUDA
373 
379  void
380  import(const CUDAWrappers::Vector<Number> &cuda_vec,
381  VectorOperation::values operation,
382  const std::shared_ptr<const CommunicationPatternBase>
383  &communication_pattern =
384  std::shared_ptr<const CommunicationPatternBase>());
385 #endif
386 
395  size_type
396  size() const;
397 
403  size_type
404  n_elements() const;
405 
409  const IndexSet &
410  get_stored_elements() const;
411 
417  iterator
418  begin();
419 
424  const_iterator
425  begin() const;
426 
431  iterator
432  end();
433 
438  const_iterator
439  end() const;
441 
442 
447 
453  Number
454  operator()(const size_type global_index) const;
455 
461  Number &
462  operator()(const size_type global_index);
463 
471  Number operator[](const size_type global_index) const;
472 
480  Number &operator[](const size_type global_index);
481 
497  template <typename Number2>
498  void
499  extract_subvector_to(const std::vector<size_type> &indices,
500  std::vector<Number2> & values) const;
501 
529  template <typename ForwardIterator, typename OutputIterator>
530  void
531  extract_subvector_to(ForwardIterator indices_begin,
532  const ForwardIterator indices_end,
533  OutputIterator values_begin) const;
534 
545  Number
546  local_element(const size_type local_index) const;
547 
558  Number &
559  local_element(const size_type local_index);
561 
562 
567 
572  template <typename Number2>
573  void
574  add(const std::vector<size_type> &indices,
575  const std::vector<Number2> & values);
576 
581  template <typename Number2>
582  void
583  add(const std::vector<size_type> & indices,
585 
591  template <typename Number2>
592  void
593  add(const size_type n_elements,
594  const size_type *indices,
595  const Number2 * values);
596 
600  void
601  print(std::ostream & out,
602  const unsigned int precision = 3,
603  const bool scientific = true) const;
604 
608  std::size_t
609  memory_consumption() const;
611 
612  protected:
613 #ifdef DEAL_II_WITH_TRILINOS
614 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
615 
620  void
621  import(
622  const Tpetra::Vector<Number, int, types::global_dof_index> &tpetra_vector,
623  const IndexSet & locally_owned_elements,
624  VectorOperation::values operation,
625  const MPI_Comm & mpi_comm,
626  const std::shared_ptr<const CommunicationPatternBase>
627  &communication_pattern);
628 # endif
629 
635  void
636  import(const Epetra_MultiVector &multivector,
637  const IndexSet & locally_owned_elements,
638  VectorOperation::values operation,
639  const MPI_Comm & mpi_comm,
640  const std::shared_ptr<const CommunicationPatternBase>
641  &communication_pattern);
642 #endif
643 
647  unsigned int
648  global_to_local(const types::global_dof_index global_index) const
649  {
650  // the following will throw an exception if the global_index is not
651  // in the remaining_elements
652  return static_cast<unsigned int>(
653  stored_elements.index_within_set(global_index));
654  }
655 
659  void
660  resize_val(const size_type new_allocated_size);
661 
662 #if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI)
663 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
664 
668  TpetraWrappers::CommunicationPattern
669  create_tpetra_comm_pattern(const IndexSet &source_index_set,
670  const MPI_Comm &mpi_comm);
671 # endif
672 
678  create_epetra_comm_pattern(const IndexSet &source_index_set,
679  const MPI_Comm &mpi_comm);
680 #endif
681 
686 
691 
696  std::shared_ptr<CommunicationPatternBase> comm_pattern;
697 
701  std::unique_ptr<Number[], decltype(free) *> values;
702 
707  mutable std::shared_ptr<::parallel::internal::TBBPartitioner>
709 
713  template <typename Number2>
714  friend class ReadWriteVector;
715 
716  private:
722  template <typename Functor>
724  {
725  public:
730 
734  virtual void
735  operator()(const size_type begin, const size_type end);
736 
737  private:
742 
746  const Functor &functor;
747  };
748  };
749 
753  /*---------------------------- Inline functions ---------------------------*/
754 
755 #ifndef DOXYGEN
756 
757  template <typename Number>
759  : Subscriptor()
760  , values(nullptr, free)
761  {
762  // virtual functions called in constructors and destructors never use the
763  // override in a derived class
764  // for clarity be explicit on which function is called
766  }
767 
768 
769 
770  template <typename Number>
772  const ReadWriteVector<Number> &v)
773  : Subscriptor()
774  , values(nullptr, free)
775  {
776  this->operator=(v);
777  }
778 
779 
780 
781  template <typename Number>
782  inline ReadWriteVector<Number>::ReadWriteVector(const size_type size)
783  : Subscriptor()
784  , values(nullptr, free)
785  {
786  // virtual functions called in constructors and destructors never use the
787  // override in a derived class
788  // for clarity be explicit on which function is called
790  }
791 
792 
793 
794  template <typename Number>
796  const IndexSet &locally_stored_indices)
797  : Subscriptor()
798  , values(nullptr, free)
799  {
800  // virtual functions called in constructors and destructors never use the
801  // override in a derived class
802  // for clarity be explicit on which function is called
803  ReadWriteVector<Number>::reinit(locally_stored_indices);
804  }
805 
806 
807 
808  template <typename Number>
809  inline typename ReadWriteVector<Number>::size_type
811  {
812  return stored_elements.size();
813  }
814 
815 
816 
817  template <typename Number>
818  inline typename ReadWriteVector<Number>::size_type
820  {
821  return stored_elements.n_elements();
822  }
823 
824 
825 
826  template <typename Number>
827  inline const IndexSet &
829  {
830  return stored_elements;
831  }
832 
833 
834 
835  template <typename Number>
836  inline typename ReadWriteVector<Number>::iterator
838  {
839  return values.get();
840  }
841 
842 
843 
844  template <typename Number>
845  inline typename ReadWriteVector<Number>::const_iterator
847  {
848  return values.get();
849  }
850 
851 
852 
853  template <typename Number>
854  inline typename ReadWriteVector<Number>::iterator
856  {
857  return values.get() + this->n_elements();
858  }
859 
860 
861 
862  template <typename Number>
863  inline typename ReadWriteVector<Number>::const_iterator
865  {
866  return values.get() + this->n_elements();
867  }
868 
869 
870 
871  template <typename Number>
872  inline Number
873  ReadWriteVector<Number>::operator()(const size_type global_index) const
874  {
875  return values[global_to_local(global_index)];
876  }
877 
878 
879 
880  template <typename Number>
881  inline Number &
882  ReadWriteVector<Number>::operator()(const size_type global_index)
883  {
884  return values[global_to_local(global_index)];
885  }
886 
887 
888 
889  template <typename Number>
890  inline Number ReadWriteVector<Number>::
891  operator[](const size_type global_index) const
892  {
893  return operator()(global_index);
894  }
895 
896 
897 
898  template <typename Number>
899  inline Number &ReadWriteVector<Number>::
900  operator[](const size_type global_index)
901  {
902  return operator()(global_index);
903  }
904 
905 
906 
907  template <typename Number>
908  template <typename Number2>
909  inline void
911  const std::vector<size_type> &indices,
912  std::vector<Number2> & extracted_values) const
913  {
914  for (size_type i = 0; i < indices.size(); ++i)
915  extracted_values[i] = operator()(indices[i]);
916  }
917 
918 
919 
920  template <typename Number>
921  template <typename ForwardIterator, typename OutputIterator>
922  inline void
924  ForwardIterator indices_begin,
925  const ForwardIterator indices_end,
926  OutputIterator values_begin) const
927  {
928  while (indices_begin != indices_end)
929  {
930  *values_begin = operator()(*indices_begin);
931  indices_begin++;
932  values_begin++;
933  }
934  }
935 
936 
937 
938  template <typename Number>
939  inline Number
940  ReadWriteVector<Number>::local_element(const size_type local_index) const
941  {
942  AssertIndexRange(local_index, this->n_elements());
943 
944  return values[local_index];
945  }
946 
947 
948 
949  template <typename Number>
950  inline Number &
951  ReadWriteVector<Number>::local_element(const size_type local_index)
952  {
953  AssertIndexRange(local_index, this->n_elements());
954 
955  return values[local_index];
956  }
957 
958 
959 
960  template <typename Number>
961  template <typename Number2>
962  inline void
963  ReadWriteVector<Number>::add(const std::vector<size_type> &indices,
964  const std::vector<Number2> & values)
965  {
966  AssertDimension(indices.size(), values.size());
967  add(indices.size(), indices.data(), values.data());
968  }
969 
970 
971 
972  template <typename Number>
973  template <typename Number2>
974  inline void
975  ReadWriteVector<Number>::add(const std::vector<size_type> & indices,
976  const ReadWriteVector<Number2> &values)
977  {
978  const size_type size = indices.size();
979  for (size_type i = 0; i < size; ++i)
980  {
981  Assert(
983  ExcMessage(
984  "The given value is not finite but either infinite or Not A Number (NaN)"));
985  this->operator()(indices[i]) += values[indices[i]];
986  }
987  }
988 
989 
990 
991  template <typename Number>
992  template <typename Number2>
993  inline void
994  ReadWriteVector<Number>::add(const size_type n_indices,
995  const size_type *indices,
996  const Number2 * values_to_add)
997  {
998  for (size_type i = 0; i < n_indices; ++i)
999  {
1000  Assert(
1002  ExcMessage(
1003  "The given value is not finite but either infinite or Not A Number (NaN)"));
1004  this->operator()(indices[i]) += values_to_add[i];
1005  }
1006  }
1007 
1008 
1009 
1010  template <typename Number>
1011  template <typename Functor>
1013  ReadWriteVector<Number> &parent,
1014  const Functor & functor)
1015  : parent(parent)
1016  , functor(functor)
1017  {}
1018 
1019 
1020 
1021  template <typename Number>
1022  template <typename Functor>
1023  void
1025  operator()(const size_type begin, const size_type end)
1026  {
1027  for (size_type i = begin; i < end; ++i)
1028  functor(parent.values[i]);
1029  }
1030 
1031 #endif // ifndef DOXYGEN
1032 
1033 } // end of namespace LinearAlgebra
1034 
1035 
1036 
1044 template <typename Number>
1045 inline void
1048 {
1049  u.swap(v);
1050 }
1051 
1052 
1053 DEAL_II_NAMESPACE_CLOSE
1054 
1055 #endif
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
Number operator[](const size_type global_index) const
void swap(ReadWriteVector< Number > &v)
Number local_element(const size_type local_index) const
FunctorTemplate(ReadWriteVector< Number > &parent, const Functor &functor)
void resize_val(const size_type new_allocated_size)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
std::size_t memory_consumption() const
~ReadWriteVector() override=default
void swap(LinearAlgebra::ReadWriteVector< Number > &u, LinearAlgebra::ReadWriteVector< Number > &v)
LinearAlgebra::distributed::Vector< Number > Vector
bool is_finite(const double x)
Definition: numbers.h:463
virtual void reinit(const size_type size, const bool omit_zeroing_entries=false)
size_type size() const
Definition: index_set.h:1600
std::unique_ptr< Number[], decltype(free) * > values
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1821
Number operator()(const size_type global_index) const
ReadWriteVector< Number > & operator=(const ReadWriteVector< Number > &in_vector)
EpetraWrappers::CommunicationPattern create_epetra_comm_pattern(const IndexSet &source_index_set, const MPI_Comm &mpi_comm)
void add(const std::vector< size_type > &indices, const std::vector< Number2 > &values)
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true) const
virtual void operator()(const size_type begin, const size_type end)
std::shared_ptr< CommunicationPatternBase > comm_pattern
unsigned int global_dof_index
Definition: types.h:89
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< Number2 > &values) const
void apply(const Functor &func)
void free(T *&pointer)
Definition: cuda.h:96
unsigned int global_to_local(const types::global_dof_index global_index) const
size_type n_elements() const
size_type n_elements() const
Definition: index_set.h:1732
const IndexSet & get_stored_elements() const
std::shared_ptr<::parallel::internal::TBBPartitioner > thread_loop_partitioner