Reference documentation for deal.II version 9.3.0
\(\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 # ifdef DEAL_II_WITH_MPI
222 
233  void
234  reinit(const TrilinosWrappers::MPI::Vector &trilinos_vec);
235 # endif
236 #endif
237 
251  template <typename Functor>
252  void
253  apply(const Functor &func);
254 
267  void
269 
274  operator=(const ReadWriteVector<Number> &in_vector);
275 
279  template <typename Number2>
281  operator=(const ReadWriteVector<Number2> &in_vector);
282 
288  operator=(const Number s);
289 
299  void
300  import(const ::Vector<Number> &vec,
301  VectorOperation::values operation,
302  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
303  &communication_pattern = {});
304 
314  void
315  import(const LinearAlgebra::Vector<Number> &vec,
316  VectorOperation::values operation,
317  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
318  &communication_pattern = {});
319 
328  template <typename MemorySpace>
329  void
331  VectorOperation::values operation,
332  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
333  &communication_pattern = {});
334 
335 #ifdef DEAL_II_WITH_PETSC
336 
344  void
345  import(const PETScWrappers::MPI::Vector &petsc_vec,
346  VectorOperation::values operation,
347  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
348  &communication_pattern = {});
349 #endif
350 
351 #ifdef DEAL_II_WITH_TRILINOS
352 
362  void
363  import(const TrilinosWrappers::MPI::Vector &trilinos_vec,
364  VectorOperation::values operation,
365  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
366  &communication_pattern = {});
367 
368 # ifdef DEAL_II_WITH_MPI
369 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
370 
378  void
379  import(const TpetraWrappers::Vector<Number> &tpetra_vec,
380  VectorOperation::values operation,
381  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
382  &communication_pattern = {});
383 # endif
384 
393  void
394  import(const EpetraWrappers::Vector &epetra_vec,
395  VectorOperation::values operation,
396  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
397  &communication_pattern = {});
398 # endif
399 #endif
400 
401 #ifdef DEAL_II_WITH_CUDA
402 
408  void
409  import(const CUDAWrappers::Vector<Number> &cuda_vec,
410  VectorOperation::values operation,
411  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
412  &communication_pattern = {});
413 #endif
414 
423  size_type
424  size() const;
425 
434  size_type
435  n_elements() const;
436 
441  size_type
442  locally_owned_size() const;
443 
447  const IndexSet &
448  get_stored_elements() const;
449 
455  iterator
456  begin();
457 
463  begin() const;
464 
469  iterator
470  end();
471 
477  end() const;
479 
480 
485 
491  Number
492  operator()(const size_type global_index) const;
493 
499  Number &
500  operator()(const size_type global_index);
501 
509  Number operator[](const size_type global_index) const;
510 
518  Number &operator[](const size_type global_index);
519 
535  template <typename Number2>
536  void
537  extract_subvector_to(const std::vector<size_type> &indices,
538  std::vector<Number2> & values) const;
539 
567  template <typename ForwardIterator, typename OutputIterator>
568  void
569  extract_subvector_to(ForwardIterator indices_begin,
570  const ForwardIterator indices_end,
571  OutputIterator values_begin) const;
572 
583  Number
584  local_element(const size_type local_index) const;
585 
596  Number &
597  local_element(const size_type local_index);
599 
600 
605 
610  template <typename Number2>
611  void
612  add(const std::vector<size_type> &indices,
613  const std::vector<Number2> & values);
614 
619  template <typename Number2>
620  void
621  add(const std::vector<size_type> & indices,
623 
629  template <typename Number2>
630  void
631  add(const size_type n_elements,
632  const size_type *indices,
633  const Number2 * values);
634 
638  void
639  print(std::ostream & out,
640  const unsigned int precision = 3,
641  const bool scientific = true) const;
642 
646  std::size_t
647  memory_consumption() const;
649 
650  protected:
651 #ifdef DEAL_II_WITH_TRILINOS
652 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
653 
658  void
659  import(
660  const Tpetra::Vector<Number, int, types::global_dof_index> &tpetra_vector,
661  const IndexSet & locally_owned_elements,
662  VectorOperation::values operation,
663  const MPI_Comm & mpi_comm,
664  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
665  &communication_pattern);
666 # endif
667 
673  void
674  import(const Epetra_MultiVector &multivector,
675  const IndexSet & locally_owned_elements,
676  VectorOperation::values operation,
677  const MPI_Comm & mpi_comm,
678  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
679  &communication_pattern);
680 #endif
681 
685  unsigned int
687  {
688  // the following will throw an exception if the global_index is not
689  // in the remaining_elements
690  return static_cast<unsigned int>(
691  stored_elements.index_within_set(global_index));
692  }
693 
697  void
698  resize_val(const size_type new_allocated_size);
699 
700 #if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI)
701 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
702 
707  create_tpetra_comm_pattern(const IndexSet &source_index_set,
708  const MPI_Comm &mpi_comm);
709 # endif
710 
716  create_epetra_comm_pattern(const IndexSet &source_index_set,
717  const MPI_Comm &mpi_comm);
718 #endif
719 
724 
729 
734  std::shared_ptr<Utilities::MPI::CommunicationPatternBase> comm_pattern;
735 
739  std::unique_ptr<Number[], decltype(std::free) *> values;
740 
745  mutable std::shared_ptr<::parallel::internal::TBBPartitioner>
747 
748  // Make all other ReadWriteVector types friends.
749  template <typename Number2>
750  friend class ReadWriteVector;
751 
752  private:
758  template <typename Functor>
760  {
761  public:
765  FunctorTemplate(ReadWriteVector<Number> &parent, const Functor &functor);
766 
770  virtual void
771  operator()(const size_type begin, const size_type end);
772 
773  private:
778 
782  const Functor &functor;
783  };
784  };
785 
789  /*---------------------------- Inline functions ---------------------------*/
790 
791 #ifndef DOXYGEN
792 
793  template <typename Number>
795  : Subscriptor()
796  , values(nullptr, free)
797  {
798  // virtual functions called in constructors and destructors never use the
799  // override in a derived class
800  // for clarity be explicit on which function is called
802  }
803 
804 
805 
806  template <typename Number>
808  const ReadWriteVector<Number> &v)
809  : Subscriptor()
810  , values(nullptr, free)
811  {
812  this->operator=(v);
813  }
814 
815 
816 
817  template <typename Number>
819  : Subscriptor()
820  , values(nullptr, free)
821  {
822  // virtual functions called in constructors and destructors never use the
823  // override in a derived class
824  // for clarity be explicit on which function is called
825  ReadWriteVector<Number>::reinit(size, false);
826  }
827 
828 
829 
830  template <typename Number>
832  const IndexSet &locally_stored_indices)
833  : Subscriptor()
834  , values(nullptr, free)
835  {
836  // virtual functions called in constructors and destructors never use the
837  // override in a derived class
838  // for clarity be explicit on which function is called
839  ReadWriteVector<Number>::reinit(locally_stored_indices);
840  }
841 
842 
843 
844  template <typename Number>
845  inline typename ReadWriteVector<Number>::size_type
847  {
848  return stored_elements.size();
849  }
850 
851 
852 
853  template <typename Number>
854  inline typename ReadWriteVector<Number>::size_type
856  {
857  return stored_elements.n_elements();
858  }
859 
860 
861 
862  template <typename Number>
863  inline typename ReadWriteVector<Number>::size_type
865  {
866  return stored_elements.n_elements();
867  }
868 
869 
870 
871  template <typename Number>
872  inline const IndexSet &
874  {
875  return stored_elements;
876  }
877 
878 
879 
880  template <typename Number>
881  inline typename ReadWriteVector<Number>::iterator
883  {
884  return values.get();
885  }
886 
887 
888 
889  template <typename Number>
892  {
893  return values.get();
894  }
895 
896 
897 
898  template <typename Number>
899  inline typename ReadWriteVector<Number>::iterator
901  {
902  return values.get() + this->n_elements();
903  }
904 
905 
906 
907  template <typename Number>
910  {
911  return values.get() + this->n_elements();
912  }
913 
914 
915 
916  template <typename Number>
917  inline Number
919  {
920  return values[global_to_local(global_index)];
921  }
922 
923 
924 
925  template <typename Number>
926  inline Number &
928  {
929  return values[global_to_local(global_index)];
930  }
931 
932 
933 
934  template <typename Number>
935  inline Number ReadWriteVector<Number>::
936  operator[](const size_type global_index) const
937  {
938  return operator()(global_index);
939  }
940 
941 
942 
943  template <typename Number>
944  inline Number &ReadWriteVector<Number>::
945  operator[](const size_type global_index)
946  {
947  return operator()(global_index);
948  }
949 
950 
951 
952  template <typename Number>
953  template <typename Number2>
954  inline void
956  const std::vector<size_type> &indices,
957  std::vector<Number2> & extracted_values) const
958  {
959  for (size_type i = 0; i < indices.size(); ++i)
960  extracted_values[i] = operator()(indices[i]);
961  }
962 
963 
964 
965  template <typename Number>
966  template <typename ForwardIterator, typename OutputIterator>
967  inline void
969  ForwardIterator indices_begin,
970  const ForwardIterator indices_end,
971  OutputIterator values_begin) const
972  {
973  while (indices_begin != indices_end)
974  {
975  *values_begin = operator()(*indices_begin);
976  indices_begin++;
977  values_begin++;
978  }
979  }
980 
981 
982 
983  template <typename Number>
984  inline Number
985  ReadWriteVector<Number>::local_element(const size_type local_index) const
986  {
987  AssertIndexRange(local_index, this->n_elements());
988 
989  return values[local_index];
990  }
991 
992 
993 
994  template <typename Number>
995  inline Number &
997  {
998  AssertIndexRange(local_index, this->n_elements());
999 
1000  return values[local_index];
1001  }
1002 
1003 
1004 
1005  template <typename Number>
1006  template <typename Number2>
1007  inline void
1008  ReadWriteVector<Number>::add(const std::vector<size_type> &indices,
1009  const std::vector<Number2> & values)
1010  {
1011  AssertDimension(indices.size(), values.size());
1012  add(indices.size(), indices.data(), values.data());
1013  }
1014 
1015 
1016 
1017  template <typename Number>
1018  template <typename Number2>
1019  inline void
1020  ReadWriteVector<Number>::add(const std::vector<size_type> & indices,
1021  const ReadWriteVector<Number2> &values)
1022  {
1023  const size_type size = indices.size();
1024  for (size_type i = 0; i < size; ++i)
1025  {
1026  Assert(
1027  numbers::is_finite(values[i]),
1028  ExcMessage(
1029  "The given value is not finite but either infinite or Not A Number (NaN)"));
1030  this->operator()(indices[i]) += values[indices[i]];
1031  }
1032  }
1033 
1034 
1035 
1036  template <typename Number>
1037  template <typename Number2>
1038  inline void
1039  ReadWriteVector<Number>::add(const size_type n_indices,
1040  const size_type *indices,
1041  const Number2 * values_to_add)
1042  {
1043  for (size_type i = 0; i < n_indices; ++i)
1044  {
1045  Assert(
1046  numbers::is_finite(values[i]),
1047  ExcMessage(
1048  "The given value is not finite but either infinite or Not A Number (NaN)"));
1049  this->operator()(indices[i]) += values_to_add[i];
1050  }
1051  }
1052 
1053 
1054 
1055  template <typename Number>
1056  template <typename Functor>
1059  const Functor & functor)
1060  : parent(parent)
1061  , functor(functor)
1062  {}
1063 
1064 
1065 
1066  template <typename Number>
1067  template <typename Functor>
1068  void
1070  operator()(const size_type begin, const size_type end)
1071  {
1072  for (size_type i = begin; i < end; ++i)
1073  functor(parent.values[i]);
1074  }
1075 
1076 #endif // ifndef DOXYGEN
1077 
1078 } // end of namespace LinearAlgebra
1079 
1080 
1081 
1089 template <typename Number>
1090 inline void
1093 {
1094  u.swap(v);
1095 }
1096 
1097 
1099 
1100 #endif
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
Number operator[](const size_type global_index) const
#define DEAL_II_DEPRECATED_EARLY
Definition: config.h:165
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:1690
bool is_finite(const double x)
Definition: numbers.h:527
virtual void reinit(const size_type size, const bool omit_zeroing_entries=false)
size_type size() const
Definition: index_set.h:1634
static ::ExceptionBase & ExcMessage(std::string arg1)
void swap(LinearAlgebra::ReadWriteVector< Number > &u, LinearAlgebra::ReadWriteVector< Number > &v)
#define Assert(cond, exc)
Definition: exceptions.h:1465
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:395
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:394
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
size_type n_elements() const
Definition: index_set.h:1832
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