Reference documentation for deal.II version Git 497f915867 2021-09-17 22:46:48 +0200
\(\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 - 2021 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.h>
25 #include <deal.II/base/parallel.h>
28 #include <deal.II/base/types.h>
29 #include <deal.II/base/utilities.h>
30 
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 namespace LinearAlgebra
51 {
52  template <typename>
53  class Vector;
54  namespace distributed
55  {
56  template <typename, typename>
57  class Vector;
58  } // namespace distributed
59 } // namespace LinearAlgebra
60 
61 # ifdef DEAL_II_WITH_PETSC
62 namespace PETScWrappers
63 {
64  namespace MPI
65  {
66  class Vector;
67  }
68 } // namespace PETScWrappers
69 # endif
70 
71 # ifdef DEAL_II_WITH_TRILINOS
72 namespace TrilinosWrappers
73 {
74  namespace MPI
75  {
76  class Vector;
77  }
78 } // namespace TrilinosWrappers
79 # endif
80 
81 # ifdef DEAL_II_WITH_CUDA
82 namespace LinearAlgebra
83 {
84  namespace CUDAWrappers
85  {
86  template <typename>
87  class Vector;
88  }
89 } // namespace LinearAlgebra
90 # endif
91 #endif
92 
93 namespace LinearAlgebra
94 {
132  template <typename Number>
134  {
135  public:
141  using value_type = Number;
142  using pointer = value_type *;
143  using const_pointer = const value_type *;
144  using iterator = value_type *;
145  using const_iterator = const value_type *;
147  using const_reference = const value_type &;
150 
158  ReadWriteVector();
159 
163  ReadWriteVector(const ReadWriteVector<Number> &in_vector);
164 
169  explicit ReadWriteVector(const size_type size);
170 
175  explicit ReadWriteVector(const IndexSet &locally_stored_indices);
176 
180  ~ReadWriteVector() override = default;
181 
190  virtual void
191  reinit(const size_type size, const bool omit_zeroing_entries = false);
192 
201  template <typename Number2>
202  void
203  reinit(const ReadWriteVector<Number2> &in_vector,
204  const bool omit_zeroing_entries = false);
205 
215  virtual void
216  reinit(const IndexSet &locally_stored_indices,
217  const bool omit_zeroing_entries = false);
218 
219 
220 #ifdef DEAL_II_WITH_TRILINOS
221 
232  void
233  reinit(const TrilinosWrappers::MPI::Vector &trilinos_vec);
234 #endif
235 
249  template <typename Functor>
250  void
251  apply(const Functor &func);
252 
265  void
267 
272  operator=(const ReadWriteVector<Number> &in_vector);
273 
277  template <typename Number2>
279  operator=(const ReadWriteVector<Number2> &in_vector);
280 
286  operator=(const Number s);
287 
297  void
298  import(const ::Vector<Number> &vec,
299  VectorOperation::values operation,
300  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
301  &communication_pattern = {});
302 
312  void
313  import(const LinearAlgebra::Vector<Number> &vec,
314  VectorOperation::values operation,
315  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
316  &communication_pattern = {});
317 
326  template <typename MemorySpace>
327  void
329  VectorOperation::values operation,
330  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
331  &communication_pattern = {});
332 
333 #ifdef DEAL_II_WITH_PETSC
334 
342  void
343  import(const PETScWrappers::MPI::Vector &petsc_vec,
344  VectorOperation::values operation,
345  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
346  &communication_pattern = {});
347 #endif
348 
349 #ifdef DEAL_II_WITH_TRILINOS
350 
360  void
361  import(const TrilinosWrappers::MPI::Vector &trilinos_vec,
362  VectorOperation::values operation,
363  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
364  &communication_pattern = {});
365 
366 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
367 
375  void
376  import(const TpetraWrappers::Vector<Number> &tpetra_vec,
377  VectorOperation::values operation,
378  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
379  &communication_pattern = {});
380 # endif
381 
390  void
391  import(const EpetraWrappers::Vector &epetra_vec,
392  VectorOperation::values operation,
393  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
394  &communication_pattern = {});
395 #endif
396 
397 #ifdef DEAL_II_WITH_CUDA
398 
404  void
405  import(const CUDAWrappers::Vector<Number> &cuda_vec,
406  VectorOperation::values operation,
407  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
408  &communication_pattern = {});
409 #endif
410 
419  size_type
420  size() const;
421 
430  size_type
431  n_elements() const;
432 
437  size_type
438  locally_owned_size() const;
439 
443  const IndexSet &
444  get_stored_elements() const;
445 
451  iterator
452  begin();
453 
459  begin() const;
460 
465  iterator
466  end();
467 
473  end() const;
475 
476 
481 
487  Number
488  operator()(const size_type global_index) const;
489 
495  Number &
496  operator()(const size_type global_index);
497 
505  Number
506  operator[](const size_type global_index) const;
507 
515  Number &
516  operator[](const size_type global_index);
517 
533  template <typename Number2>
534  void
535  extract_subvector_to(const std::vector<size_type> &indices,
536  std::vector<Number2> & values) const;
537 
565  template <typename ForwardIterator, typename OutputIterator>
566  void
567  extract_subvector_to(ForwardIterator indices_begin,
568  const ForwardIterator indices_end,
569  OutputIterator values_begin) const;
570 
581  Number
582  local_element(const size_type local_index) const;
583 
594  Number &
595  local_element(const size_type local_index);
597 
598 
603 
608  template <typename Number2>
609  void
610  add(const std::vector<size_type> &indices,
611  const std::vector<Number2> & values);
612 
617  template <typename Number2>
618  void
619  add(const std::vector<size_type> & indices,
621 
627  template <typename Number2>
628  void
629  add(const size_type n_elements,
630  const size_type *indices,
631  const Number2 * values);
632 
636  void
637  print(std::ostream & out,
638  const unsigned int precision = 3,
639  const bool scientific = true) const;
640 
644  std::size_t
645  memory_consumption() const;
647 
648  protected:
649 #ifdef DEAL_II_WITH_TRILINOS
650 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
651 
656  void
657  import(
658  const Tpetra::Vector<Number, int, types::global_dof_index> &tpetra_vector,
659  const IndexSet & locally_owned_elements,
660  VectorOperation::values operation,
661  const MPI_Comm & mpi_comm,
662  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
663  &communication_pattern);
664 # endif
665 
671  void
672  import(const Epetra_MultiVector &multivector,
673  const IndexSet & locally_owned_elements,
674  VectorOperation::values operation,
675  const MPI_Comm & mpi_comm,
676  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
677  &communication_pattern);
678 #endif
679 
683  unsigned int
685  {
686  // the following will throw an exception if the global_index is not
687  // in the remaining_elements
688  return static_cast<unsigned int>(
689  stored_elements.index_within_set(global_index));
690  }
691 
695  void
696  resize_val(const size_type new_allocated_size);
697 
698 #ifdef DEAL_II_WITH_TRILINOS
699 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
700 
705  create_tpetra_comm_pattern(const IndexSet &source_index_set,
706  const MPI_Comm &mpi_comm);
707 # endif
708 
714  create_epetra_comm_pattern(const IndexSet &source_index_set,
715  const MPI_Comm &mpi_comm);
716 #endif
717 
722 
727 
732  std::shared_ptr<Utilities::MPI::CommunicationPatternBase> comm_pattern;
733 
737  std::unique_ptr<Number[], decltype(std::free) *> values;
738 
743  mutable std::shared_ptr<::parallel::internal::TBBPartitioner>
745 
746  // Make all other ReadWriteVector types friends.
747  template <typename Number2>
748  friend class ReadWriteVector;
749 
750  private:
756  template <typename Functor>
758  {
759  public:
763  FunctorTemplate(ReadWriteVector<Number> &parent, const Functor &functor);
764 
768  virtual void
769  operator()(const size_type begin, const size_type end);
770 
771  private:
776 
780  const Functor &functor;
781  };
782  };
783 
787  /*---------------------------- Inline functions ---------------------------*/
788 
789 #ifndef DOXYGEN
790 
791  template <typename Number>
793  : Subscriptor()
794  , values(nullptr, free)
795  {
796  // virtual functions called in constructors and destructors never use the
797  // override in a derived class
798  // for clarity be explicit on which function is called
800  }
801 
802 
803 
804  template <typename Number>
806  const ReadWriteVector<Number> &v)
807  : Subscriptor()
808  , values(nullptr, free)
809  {
810  this->operator=(v);
811  }
812 
813 
814 
815  template <typename Number>
817  : Subscriptor()
818  , values(nullptr, free)
819  {
820  // virtual functions called in constructors and destructors never use the
821  // override in a derived class
822  // for clarity be explicit on which function is called
823  ReadWriteVector<Number>::reinit(size, false);
824  }
825 
826 
827 
828  template <typename Number>
830  const IndexSet &locally_stored_indices)
831  : Subscriptor()
832  , values(nullptr, free)
833  {
834  // virtual functions called in constructors and destructors never use the
835  // override in a derived class
836  // for clarity be explicit on which function is called
837  ReadWriteVector<Number>::reinit(locally_stored_indices);
838  }
839 
840 
841 
842  template <typename Number>
843  inline typename ReadWriteVector<Number>::size_type
845  {
846  return stored_elements.size();
847  }
848 
849 
850 
851  template <typename Number>
852  inline typename ReadWriteVector<Number>::size_type
854  {
855  return stored_elements.n_elements();
856  }
857 
858 
859 
860  template <typename Number>
861  inline typename ReadWriteVector<Number>::size_type
863  {
864  return stored_elements.n_elements();
865  }
866 
867 
868 
869  template <typename Number>
870  inline const IndexSet &
872  {
873  return stored_elements;
874  }
875 
876 
877 
878  template <typename Number>
879  inline typename ReadWriteVector<Number>::iterator
881  {
882  return values.get();
883  }
884 
885 
886 
887  template <typename Number>
890  {
891  return values.get();
892  }
893 
894 
895 
896  template <typename Number>
897  inline typename ReadWriteVector<Number>::iterator
899  {
900  return values.get() + this->n_elements();
901  }
902 
903 
904 
905  template <typename Number>
908  {
909  return values.get() + this->n_elements();
910  }
911 
912 
913 
914  template <typename Number>
915  inline Number
917  {
918  return values[global_to_local(global_index)];
919  }
920 
921 
922 
923  template <typename Number>
924  inline Number &
926  {
927  return values[global_to_local(global_index)];
928  }
929 
930 
931 
932  template <typename Number>
933  inline Number
934  ReadWriteVector<Number>::operator[](const size_type global_index) const
935  {
936  return operator()(global_index);
937  }
938 
939 
940 
941  template <typename Number>
942  inline Number &
944  {
945  return operator()(global_index);
946  }
947 
948 
949 
950  template <typename Number>
951  template <typename Number2>
952  inline void
954  const std::vector<size_type> &indices,
955  std::vector<Number2> & extracted_values) const
956  {
957  for (size_type i = 0; i < indices.size(); ++i)
958  extracted_values[i] = operator()(indices[i]);
959  }
960 
961 
962 
963  template <typename Number>
964  template <typename ForwardIterator, typename OutputIterator>
965  inline void
967  ForwardIterator indices_begin,
968  const ForwardIterator indices_end,
969  OutputIterator values_begin) const
970  {
971  while (indices_begin != indices_end)
972  {
973  *values_begin = operator()(*indices_begin);
974  indices_begin++;
975  values_begin++;
976  }
977  }
978 
979 
980 
981  template <typename Number>
982  inline Number
983  ReadWriteVector<Number>::local_element(const size_type local_index) const
984  {
985  AssertIndexRange(local_index, this->n_elements());
986 
987  return values[local_index];
988  }
989 
990 
991 
992  template <typename Number>
993  inline Number &
995  {
996  AssertIndexRange(local_index, this->n_elements());
997 
998  return values[local_index];
999  }
1000 
1001 
1002 
1003  template <typename Number>
1004  template <typename Number2>
1005  inline void
1006  ReadWriteVector<Number>::add(const std::vector<size_type> &indices,
1007  const std::vector<Number2> & values)
1008  {
1009  AssertDimension(indices.size(), values.size());
1010  add(indices.size(), indices.data(), values.data());
1011  }
1012 
1013 
1014 
1015  template <typename Number>
1016  template <typename Number2>
1017  inline void
1018  ReadWriteVector<Number>::add(const std::vector<size_type> & indices,
1019  const ReadWriteVector<Number2> &values)
1020  {
1021  const size_type size = indices.size();
1022  for (size_type i = 0; i < size; ++i)
1023  {
1024  Assert(
1025  numbers::is_finite(values[i]),
1026  ExcMessage(
1027  "The given value is not finite but either infinite or Not A Number (NaN)"));
1028  this->operator()(indices[i]) += values[indices[i]];
1029  }
1030  }
1031 
1032 
1033 
1034  template <typename Number>
1035  template <typename Number2>
1036  inline void
1037  ReadWriteVector<Number>::add(const size_type n_indices,
1038  const size_type *indices,
1039  const Number2 * values_to_add)
1040  {
1041  for (size_type i = 0; i < n_indices; ++i)
1042  {
1043  Assert(
1044  numbers::is_finite(values[i]),
1045  ExcMessage(
1046  "The given value is not finite but either infinite or Not A Number (NaN)"));
1047  this->operator()(indices[i]) += values_to_add[i];
1048  }
1049  }
1050 
1051 
1052 
1053  template <typename Number>
1054  template <typename Functor>
1057  const Functor & functor)
1058  : parent(parent)
1059  , functor(functor)
1060  {}
1061 
1062 
1063 
1064  template <typename Number>
1065  template <typename Functor>
1066  void
1068  const size_type begin,
1069  const size_type end)
1070  {
1071  for (size_type i = begin; i < end; ++i)
1072  functor(parent.values[i]);
1073  }
1074 
1075 #endif // ifndef DOXYGEN
1076 
1077 } // end of namespace LinearAlgebra
1078 
1079 
1080 
1088 template <typename Number>
1089 inline void
1092 {
1093  u.swap(v);
1094 }
1095 
1096 
1098 
1099 #endif
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1653
Number operator[](const size_type global_index) const
void swap(ReadWriteVector< Number > &v)
Number local_element(const size_type local_index) const
std::unique_ptr< Number[], decltype(std::free) * > values
std::shared_ptr< Utilities::MPI::CommunicationPatternBase > comm_pattern
#define AssertIndexRange(index, range)
Definition: exceptions.h:1718
bool is_finite(const double x)
Definition: numbers.h:537
virtual void reinit(const size_type size, const bool omit_zeroing_entries=false)
size_type size() const
Definition: index_set.h:1639
static ::ExceptionBase & ExcMessage(std::string arg1)
void swap(LinearAlgebra::ReadWriteVector< Number > &u, LinearAlgebra::ReadWriteVector< Number > &v)
#define Assert(cond, exc)
Definition: exceptions.h:1461
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:36
size_type locally_owned_size() const
Number operator()(const size_type global_index) const
ReadWriteVector< Number > & operator=(const ReadWriteVector< Number > &in_vector)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
VectorType::value_type * end(VectorType &V)
void add(const std::vector< size_type > &indices, const std::vector< Number2 > &values)
virtual void operator()(const size_type begin, const size_type end)
unsigned int global_dof_index
Definition: types.h:76
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< Number2 > &values) const
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
void free(T *&pointer)
Definition: cuda.h:97
unsigned int global_to_local(const types::global_dof_index global_index) const
size_type n_elements() const
#define DEAL_II_DEPRECATED
Definition: config.h:160
size_type n_elements() const
Definition: index_set.h:1837
const IndexSet & get_stored_elements() const
std::shared_ptr<::parallel::internal::TBBPartitioner > thread_loop_partitioner
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
typename numbers::NumberTraits< Number >::real_type real_type