Reference documentation for deal.II version GIT e22cb6a53e 2023-09-21 14:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
read_write_vector.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2023 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 
22 #include <deal.II/base/index_set.h>
24 #include <deal.II/base/mpi_stub.h>
25 #include <deal.II/base/parallel.h>
28 #include <deal.II/base/types.h>
29 
32 
33 #include <cstdlib>
34 #include <cstring>
35 #include <iomanip>
36 
37 #ifdef DEAL_II_WITH_TRILINOS
41 
42 # include <Epetra_MultiVector.h>
43 
44 #endif
45 
47 
48 // Forward declarations
49 #ifndef DOXYGEN
50 template <typename>
51 class Vector;
52 
53 namespace LinearAlgebra
54 {
55  template <typename>
56  class Vector;
57  namespace distributed
58  {
59  template <typename, typename>
60  class Vector;
61  } // namespace distributed
62 } // namespace LinearAlgebra
63 
64 # ifdef DEAL_II_WITH_PETSC
65 namespace PETScWrappers
66 {
67  namespace MPI
68  {
69  class Vector;
70  }
71 } // namespace PETScWrappers
72 # endif
73 
74 # ifdef DEAL_II_WITH_TRILINOS
75 namespace TrilinosWrappers
76 {
77  namespace MPI
78  {
79  class Vector;
80  }
81 } // namespace TrilinosWrappers
82 # endif
83 
84 # ifdef DEAL_II_WITH_CUDA
85 namespace LinearAlgebra
86 {
87  namespace CUDAWrappers
88  {
89  template <typename>
90  class Vector;
91  }
92 } // namespace LinearAlgebra
93 # endif
94 #endif
95 
96 namespace LinearAlgebra
97 {
135  template <typename Number>
136  class ReadWriteVector : public Subscriptor, public ReadVector<Number>
137  {
138  public:
144  using value_type = Number;
145  using pointer = value_type *;
146  using const_pointer = const value_type *;
147  using iterator = value_type *;
148  using const_iterator = const value_type *;
150  using const_reference = const value_type &;
153 
162 
167 
172  explicit ReadWriteVector(const size_type size);
173 
178  explicit ReadWriteVector(const IndexSet &locally_stored_indices);
179 
183  ~ReadWriteVector() override = default;
184 
193  virtual void
194  reinit(const size_type size, const bool omit_zeroing_entries = false);
195 
204  template <typename Number2>
205  void
207  const bool omit_zeroing_entries = false);
208 
218  virtual void
219  reinit(const IndexSet &locally_stored_indices,
220  const bool omit_zeroing_entries = false);
221 
222 
223 #ifdef DEAL_II_WITH_TRILINOS
235  void
236  reinit(const TrilinosWrappers::MPI::Vector &trilinos_vec);
237 #endif
238 
252  template <typename Functor>
253  void
254  apply(const Functor &func);
255 
268  void
270 
276 
280  template <typename Number2>
283 
289  operator=(const Number s);
290 
300  void
302  const ::Vector<Number> &vec,
303  VectorOperation::values operation,
304  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
305  &communication_pattern = {});
306 
311  void
312  import(const ::Vector<Number> &V,
313  VectorOperation::values operation,
314  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
315  &communication_pattern = {})
316  {
317  import_elements(V, operation, communication_pattern);
318  }
319 
328  template <typename MemorySpace>
329  void
332  VectorOperation::values operation,
333  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
334  &communication_pattern = {});
335 
339  template <typename MemorySpace>
340  DEAL_II_DEPRECATED void
342  VectorOperation::values operation,
343  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
344  &communication_pattern = {})
345  {
346  import_elements(V, operation, communication_pattern);
347  }
348 
349 
350 #ifdef DEAL_II_WITH_PETSC
359  void
361  const PETScWrappers::MPI::Vector &petsc_vec,
362  VectorOperation::values operation,
363  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
364  &communication_pattern = {});
365 
370  void
372  VectorOperation::values operation,
373  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
374  &communication_pattern = {})
375  {
376  import_elements(V, operation, communication_pattern);
377  }
378 #endif
379 
380 #ifdef DEAL_II_WITH_TRILINOS
391  void
393  const TrilinosWrappers::MPI::Vector &trilinos_vec,
394  VectorOperation::values operation,
395  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
396  &communication_pattern = {});
397 
402  void
404  VectorOperation::values operation,
405  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
406  &communication_pattern = {})
407  {
408  import_elements(V, operation, communication_pattern);
409  }
410 
411 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
420  template <typename Dummy = Number>
421  std::enable_if_t<std::is_same_v<Dummy, Number> &&
424  const TpetraWrappers::Vector<Number> &tpetra_vec,
425  VectorOperation::values operation,
426  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
427  &communication_pattern = {});
428 
432  template <typename Dummy = Number>
433  DEAL_II_DEPRECATED std::enable_if_t<std::is_same_v<Dummy, Number> &&
436  VectorOperation::values operation,
437  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
438  &communication_pattern = {})
439  {
440  import_elements(V, operation, communication_pattern);
441  }
442 # endif
443 
452  void
454  const EpetraWrappers::Vector &epetra_vec,
455  VectorOperation::values operation,
456  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
457  &communication_pattern = {});
458 
463  void
464  import(const EpetraWrappers::Vector &V,
465  VectorOperation::values operation,
466  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
467  &communication_pattern = {})
468  {
469  import_elements(V, operation, communication_pattern);
470  }
471 #endif
472 
473 #ifdef DEAL_II_WITH_CUDA
480  void
482  const CUDAWrappers::Vector<Number> &cuda_vec,
483  VectorOperation::values operation,
484  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
485  &communication_pattern = {});
486 
487 
492  void
494  VectorOperation::values operation,
495  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
496  &communication_pattern = {})
497  {
498  import_elements(V, operation, communication_pattern);
499  }
500 #endif
501 
510  size_type
511  size() const override;
512 
517  size_type
519 
523  const IndexSet &
525 
531  iterator
532  begin();
533 
539  begin() const;
540 
545  iterator
546  end();
547 
553  end() const;
567  Number
569 
575  Number &
577 
585  Number
587 
595  Number &
597 
613  template <typename Number2>
614  void
615  extract_subvector_to(const std::vector<size_type> &indices,
616  std::vector<Number2> &values) const;
617 
621  virtual void
624  ArrayView<Number> &entries) const override;
625 
653  template <typename ForwardIterator, typename OutputIterator>
654  void
655  extract_subvector_to(ForwardIterator indices_begin,
656  const ForwardIterator indices_end,
657  OutputIterator values_begin) const;
658 
669  Number
670  local_element(const size_type local_index) const;
671 
682  Number &
683  local_element(const size_type local_index);
696  template <typename Number2>
697  void
698  add(const std::vector<size_type> &indices,
699  const std::vector<Number2> &values);
700 
705  template <typename Number2>
706  void
707  add(const std::vector<size_type> &indices,
709 
715  template <typename Number2>
716  void
717  add(const size_type n_elements,
718  const size_type *indices,
719  const Number2 *values);
720 
724  void
725  print(std::ostream &out,
726  const unsigned int precision = 3,
727  const bool scientific = true) const;
728 
732  std::size_t
736  protected:
737 #ifdef DEAL_II_WITH_TRILINOS
738 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
744  template <typename Dummy = Number>
745  std::enable_if_t<std::is_same_v<Dummy, Number> &&
748  const Tpetra::Vector<Number, int, types::signed_global_dof_index>
749  &tpetra_vector,
750  const IndexSet &locally_owned_elements,
751  VectorOperation::values operation,
752  const MPI_Comm mpi_comm,
753  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
754  &communication_pattern);
755 # endif
756 
762  void
764  const Epetra_MultiVector &multivector,
765  const IndexSet &locally_owned_elements,
766  VectorOperation::values operation,
767  const MPI_Comm mpi_comm,
768  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
769  &communication_pattern);
770 #endif
771 
775  unsigned int
777 
781  void
782  resize_val(const size_type new_allocated_size);
783 
784 #ifdef DEAL_II_WITH_TRILINOS
785 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
791  create_tpetra_comm_pattern(const IndexSet &source_index_set,
792  const MPI_Comm mpi_comm);
793 # endif
794 
800  create_epetra_comm_pattern(const IndexSet &source_index_set,
801  const MPI_Comm mpi_comm);
802 #endif
803 
808 
813 
818  std::shared_ptr<Utilities::MPI::CommunicationPatternBase> comm_pattern;
819 
823  std::unique_ptr<Number[], /*decltype(std::free) * */ void (*)(void *)>
825 
830  mutable std::shared_ptr<::parallel::internal::TBBPartitioner>
832 
833  // Make all other ReadWriteVector types friends.
834  template <typename Number2>
835  friend class ReadWriteVector;
836 
837  private:
843  template <typename Functor>
845  {
846  public:
851 
855  virtual void
857 
858  private:
863 
867  const Functor &functor;
868  };
869  };
870 
874  /*---------------------------- Inline functions ---------------------------*/
875 
876 #ifndef DOXYGEN
877 
878 
879  template <typename Number>
881  : Subscriptor()
882  , values(nullptr, free)
883  {
884  // virtual functions called in constructors and destructors never use the
885  // override in a derived class
886  // for clarity be explicit on which function is called
888  }
889 
890 
891 
892  template <typename Number>
894  const ReadWriteVector<Number> &v)
895  : Subscriptor()
896  , values(nullptr, free)
897  {
898  this->operator=(v);
899  }
900 
901 
902 
903  template <typename Number>
905  : Subscriptor()
906  , values(nullptr, free)
907  {
908  // virtual functions called in constructors and destructors never use the
909  // override in a derived class
910  // for clarity be explicit on which function is called
912  }
913 
914 
915 
916  template <typename Number>
918  const IndexSet &locally_stored_indices)
919  : Subscriptor()
920  , values(nullptr, free)
921  {
922  // virtual functions called in constructors and destructors never use the
923  // override in a derived class
924  // for clarity be explicit on which function is called
925  ReadWriteVector<Number>::reinit(locally_stored_indices);
926  }
927 
928 
929 
930  template <typename Number>
931  inline typename ReadWriteVector<Number>::size_type
933  {
934  return stored_elements.size();
935  }
936 
937 
938 
939  template <typename Number>
940  inline typename ReadWriteVector<Number>::size_type
942  {
943  return stored_elements.n_elements();
944  }
945 
946 
947 
948  template <typename Number>
949  inline const IndexSet &
951  {
952  return stored_elements;
953  }
954 
955 
956 
957  template <typename Number>
958  inline typename ReadWriteVector<Number>::iterator
960  {
961  return values.get();
962  }
963 
964 
965 
966  template <typename Number>
969  {
970  return values.get();
971  }
972 
973 
974 
975  template <typename Number>
976  inline typename ReadWriteVector<Number>::iterator
978  {
979  return values.get() + this->locally_owned_size();
980  }
981 
982 
983 
984  template <typename Number>
987  {
988  return values.get() + this->locally_owned_size();
989  }
990 
991 
992 
993  template <typename Number>
994  inline Number
996  {
998  }
999 
1000 
1001 
1002  template <typename Number>
1003  inline Number &
1005  {
1007  }
1008 
1009 
1010 
1011  template <typename Number>
1012  inline Number
1014  {
1015  return operator()(global_index);
1016  }
1017 
1018 
1019 
1020  template <typename Number>
1021  inline Number &
1023  {
1024  return operator()(global_index);
1025  }
1026 
1027 
1028 
1029  template <typename Number>
1030  template <typename Number2>
1031  inline void
1033  const std::vector<size_type> &indices,
1034  std::vector<Number2> &extracted_values) const
1035  {
1036  for (size_type i = 0; i < indices.size(); ++i)
1037  extracted_values[i] = operator()(indices[i]);
1038  }
1039 
1040 
1041 
1042  template <typename Number>
1043  void
1046  ArrayView<Number> &entries) const
1047  {
1048  AssertDimension(indices.size(), entries.size());
1049  for (unsigned int i = 0; i < indices.size(); ++i)
1050  entries[i] = (*this)[indices[i]];
1051  }
1052 
1053 
1054 
1055  template <typename Number>
1056  template <typename ForwardIterator, typename OutputIterator>
1057  inline void
1059  ForwardIterator indices_begin,
1060  const ForwardIterator indices_end,
1061  OutputIterator values_begin) const
1062  {
1063  while (indices_begin != indices_end)
1064  {
1065  *values_begin = operator()(*indices_begin);
1066  indices_begin++;
1067  values_begin++;
1068  }
1069  }
1070 
1071 
1072 
1073  template <typename Number>
1074  inline Number
1075  ReadWriteVector<Number>::local_element(const size_type local_index) const
1076  {
1077  AssertIndexRange(local_index, this->locally_owned_size());
1078 
1079  return values[local_index];
1080  }
1081 
1082 
1083 
1084  template <typename Number>
1085  inline Number &
1087  {
1088  AssertIndexRange(local_index, this->locally_owned_size());
1089 
1090  return values[local_index];
1091  }
1092 
1093 
1094 
1095  template <typename Number>
1096  template <typename Number2>
1097  inline void
1098  ReadWriteVector<Number>::add(const std::vector<size_type> &indices,
1099  const std::vector<Number2> &values)
1100  {
1101  AssertDimension(indices.size(), values.size());
1102  add(indices.size(), indices.data(), values.data());
1103  }
1104 
1105 
1106 
1107  template <typename Number>
1108  template <typename Number2>
1109  inline void
1110  ReadWriteVector<Number>::add(const std::vector<size_type> &indices,
1111  const ReadWriteVector<Number2> &values)
1112  {
1113  const size_type size = indices.size();
1114  for (size_type i = 0; i < size; ++i)
1115  {
1116  Assert(
1118  ExcMessage(
1119  "The given value is not finite but either infinite or Not A Number (NaN)"));
1120  this->operator()(indices[i]) += values[indices[i]];
1121  }
1122  }
1123 
1124 
1125 
1126  template <typename Number>
1127  template <typename Number2>
1128  inline void
1129  ReadWriteVector<Number>::add(const size_type n_indices,
1130  const size_type *indices,
1131  const Number2 *values_to_add)
1132  {
1133  for (size_type i = 0; i < n_indices; ++i)
1134  {
1135  Assert(
1137  ExcMessage(
1138  "The given value is not finite but either infinite or Not A Number (NaN)"));
1139  this->operator()(indices[i]) += values_to_add[i];
1140  }
1141  }
1142 
1143 
1144 
1145  template <typename Number>
1146  inline unsigned int
1149  {
1150  // the following will throw an exception if the global_index is not
1151  // in the remaining_elements
1152  return static_cast<unsigned int>(
1154  }
1155 
1156 
1157 
1158  template <typename Number>
1159  template <typename Functor>
1161  ReadWriteVector<Number> &parent,
1162  const Functor &functor)
1163  : parent(parent)
1164  , functor(functor)
1165  {}
1166 
1167 
1168 
1169  template <typename Number>
1170  template <typename Functor>
1171  void
1173  const size_type begin,
1174  const size_type end)
1175  {
1176  for (size_type i = begin; i < end; ++i)
1177  functor(parent.values[i]);
1178  }
1179 
1180 #endif // ifndef DOXYGEN
1181 
1182 } // end of namespace LinearAlgebra
1183 
1184 
1185 
1193 template <typename Number>
1194 inline void
1197 {
1198  u.swap(v);
1199 }
1200 
1201 
1203 
1204 #endif
std::size_t size() const
Definition: array_view.h:573
size_type size() const
Definition: index_set.h:1696
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
size_type n_elements() const
Definition: index_set.h:1854
FunctorTemplate(ReadWriteVector< Number > &parent, const Functor &functor)
virtual void operator()(const size_type begin, const size_type end)
void import_elements(const ::Vector< Number > &vec, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern={})
ReadWriteVector< Number > & operator=(const ReadWriteVector< Number2 > &in_vector)
void reinit(const ReadWriteVector< Number2 > &in_vector, const bool omit_zeroing_entries=false)
ReadWriteVector(const ReadWriteVector< Number > &in_vector)
ReadWriteVector(const size_type size)
void import_elements(const TrilinosWrappers::MPI::Vector &trilinos_vec, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern={})
size_type size() const override
void add(const std::vector< size_type > &indices, const std::vector< Number2 > &values)
Number operator()(const size_type global_index) const
void apply(const Functor &func)
std::size_t memory_consumption() const
void import_elements(const CUDAWrappers::Vector< Number > &cuda_vec, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern={})
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true) const
Number & operator()(const size_type global_index)
void add(const std::vector< size_type > &indices, const ReadWriteVector< Number2 > &values)
void import_elements(const EpetraWrappers::Vector &epetra_vec, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern={})
ReadWriteVector< Number > & operator=(const ReadWriteVector< Number > &in_vector)
ReadWriteVector< Number > & operator=(const Number s)
unsigned int global_to_local(const types::global_dof_index global_index) const
void resize_val(const size_type new_allocated_size)
Number operator[](const size_type global_index) const
Number local_element(const size_type local_index) const
void add(const size_type n_elements, const size_type *indices, const Number2 *values)
virtual void extract_subvector_to(const ArrayView< const types::global_dof_index > &indices, ArrayView< Number > &entries) const override
void import_elements(const PETScWrappers::MPI::Vector &petsc_vec, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern={})
std::enable_if_t< std::is_same_v< Dummy, Number > &&::is_tpetra_type< Number >::value > import_elements(const TpetraWrappers::Vector< Number > &tpetra_vec, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern={})
TpetraWrappers::CommunicationPattern create_tpetra_comm_pattern(const IndexSet &source_index_set, const MPI_Comm mpi_comm)
Number & local_element(const size_type local_index)
void reinit(const TrilinosWrappers::MPI::Vector &trilinos_vec)
virtual void reinit(const IndexSet &locally_stored_indices, const bool omit_zeroing_entries=false)
void extract_subvector_to(ForwardIterator indices_begin, const ForwardIterator indices_end, OutputIterator values_begin) const
virtual void reinit(const size_type size, const bool omit_zeroing_entries=false)
const IndexSet & get_stored_elements() const
size_type locally_owned_size() const
const_iterator end() const
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< Number2 > &values) const
const_iterator begin() const
~ReadWriteVector() override=default
ReadWriteVector(const IndexSet &locally_stored_indices)
std::enable_if_t< std::is_same_v< Dummy, Number > &&::is_tpetra_type< Number >::value > import_elements(const Tpetra::Vector< Number, int, types::signed_global_dof_index > &tpetra_vector, const IndexSet &locally_owned_elements, VectorOperation::values operation, const MPI_Comm mpi_comm, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern)
EpetraWrappers::CommunicationPattern create_epetra_comm_pattern(const IndexSet &source_index_set, const MPI_Comm mpi_comm)
std::unique_ptr< Number[], void(*)(void *)> values
types::global_dof_index size_type
void import_elements(const distributed::Vector< Number, MemorySpace > &vec, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern={})
void swap(ReadWriteVector< Number > &v)
std::shared_ptr< Utilities::MPI::CommunicationPatternBase > comm_pattern
std::shared_ptr<::parallel::internal::TBBPartitioner > thread_loop_partitioner
typename numbers::NumberTraits< Number >::real_type real_type
void import_elements(const Epetra_MultiVector &multivector, const IndexSet &locally_owned_elements, VectorOperation::values operation, const MPI_Comm mpi_comm, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern)
Number & operator[](const size_type global_index)
types::global_dof_index size_type
Definition: read_vector.h:44
Definition: vector.h:110
void swap(LinearAlgebra::ReadWriteVector< Number > &u, LinearAlgebra::ReadWriteVector< Number > &v)
#define DEAL_II_DEPRECATED
Definition: config.h:178
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcMessage(std::string arg1)
static const char V
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
void free(T *&pointer)
Definition: cuda.h:97
bool is_finite(const double x)
Definition: numbers.h:539
unsigned int global_dof_index
Definition: types.h:82