Reference documentation for deal.II version Git 6d63218887 2020-10-30 17:17:53 -0400
\(\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 - 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>
23 #include <deal.II/base/mpi.h>
25 #include <deal.II/base/types.h>
26 #include <deal.II/base/utilities.h>
27 
29 
30 #include <cstdlib>
31 #include <cstring>
32 #include <iomanip>
33 
34 #ifdef DEAL_II_WITH_TRILINOS
38 
39 # include <Epetra_MultiVector.h>
40 
41 #endif
42 
44 
45 // Forward declarations
46 #ifndef DOXYGEN
47 namespace LinearAlgebra
48 {
49  class CommunicationPatternBase;
50  namespace distributed
51  {
52  template <typename, typename>
53  class Vector;
54  } // namespace distributed
55 } // namespace LinearAlgebra
56 
57 # ifdef DEAL_II_WITH_PETSC
58 namespace PETScWrappers
59 {
60  namespace MPI
61  {
62  class Vector;
63  }
64 } // namespace PETScWrappers
65 # endif
66 
67 # ifdef DEAL_II_WITH_TRILINOS
68 namespace TrilinosWrappers
69 {
70  namespace MPI
71  {
72  class Vector;
73  }
74 } // namespace TrilinosWrappers
75 # endif
76 
77 # ifdef DEAL_II_WITH_CUDA
78 namespace LinearAlgebra
79 {
80  namespace CUDAWrappers
81  {
82  template <typename>
83  class Vector;
84  }
85 } // namespace LinearAlgebra
86 # endif
87 #endif
88 
89 namespace LinearAlgebra
90 {
128  template <typename Number>
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 *;
143  using const_reference = const value_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 
425  begin() const;
426 
431  iterator
432  end();
433 
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
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 
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(std::free) *> values;
702 
707  mutable std::shared_ptr<::parallel::internal::TBBPartitioner>
709 
710  // Make all other ReadWriteVector types friends.
711  template <typename Number2>
712  friend class ReadWriteVector;
713 
714  private:
720  template <typename Functor>
722  {
723  public:
727  FunctorTemplate(ReadWriteVector<Number> &parent, const Functor &functor);
728 
732  virtual void
733  operator()(const size_type begin, const size_type end);
734 
735  private:
740 
744  const Functor &functor;
745  };
746  };
747 
751  /*---------------------------- Inline functions ---------------------------*/
752 
753 #ifndef DOXYGEN
754 
755  template <typename Number>
757  : Subscriptor()
758  , values(nullptr, free)
759  {
760  // virtual functions called in constructors and destructors never use the
761  // override in a derived class
762  // for clarity be explicit on which function is called
764  }
765 
766 
767 
768  template <typename Number>
770  const ReadWriteVector<Number> &v)
771  : Subscriptor()
772  , values(nullptr, free)
773  {
774  this->operator=(v);
775  }
776 
777 
778 
779  template <typename Number>
781  : Subscriptor()
782  , values(nullptr, free)
783  {
784  // virtual functions called in constructors and destructors never use the
785  // override in a derived class
786  // for clarity be explicit on which function is called
787  ReadWriteVector<Number>::reinit(size, false);
788  }
789 
790 
791 
792  template <typename Number>
794  const IndexSet &locally_stored_indices)
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
801  ReadWriteVector<Number>::reinit(locally_stored_indices);
802  }
803 
804 
805 
806  template <typename Number>
807  inline typename ReadWriteVector<Number>::size_type
809  {
810  return stored_elements.size();
811  }
812 
813 
814 
815  template <typename Number>
816  inline typename ReadWriteVector<Number>::size_type
818  {
819  return stored_elements.n_elements();
820  }
821 
822 
823 
824  template <typename Number>
825  inline const IndexSet &
827  {
828  return stored_elements;
829  }
830 
831 
832 
833  template <typename Number>
834  inline typename ReadWriteVector<Number>::iterator
836  {
837  return values.get();
838  }
839 
840 
841 
842  template <typename Number>
845  {
846  return values.get();
847  }
848 
849 
850 
851  template <typename Number>
852  inline typename ReadWriteVector<Number>::iterator
854  {
855  return values.get() + this->n_elements();
856  }
857 
858 
859 
860  template <typename Number>
863  {
864  return values.get() + this->n_elements();
865  }
866 
867 
868 
869  template <typename Number>
870  inline Number
872  {
873  return values[global_to_local(global_index)];
874  }
875 
876 
877 
878  template <typename Number>
879  inline Number &
881  {
882  return values[global_to_local(global_index)];
883  }
884 
885 
886 
887  template <typename Number>
888  inline Number ReadWriteVector<Number>::
889  operator[](const size_type global_index) const
890  {
891  return operator()(global_index);
892  }
893 
894 
895 
896  template <typename Number>
897  inline Number &ReadWriteVector<Number>::
898  operator[](const size_type global_index)
899  {
900  return operator()(global_index);
901  }
902 
903 
904 
905  template <typename Number>
906  template <typename Number2>
907  inline void
909  const std::vector<size_type> &indices,
910  std::vector<Number2> & extracted_values) const
911  {
912  for (size_type i = 0; i < indices.size(); ++i)
913  extracted_values[i] = operator()(indices[i]);
914  }
915 
916 
917 
918  template <typename Number>
919  template <typename ForwardIterator, typename OutputIterator>
920  inline void
922  ForwardIterator indices_begin,
923  const ForwardIterator indices_end,
924  OutputIterator values_begin) const
925  {
926  while (indices_begin != indices_end)
927  {
928  *values_begin = operator()(*indices_begin);
929  indices_begin++;
930  values_begin++;
931  }
932  }
933 
934 
935 
936  template <typename Number>
937  inline Number
938  ReadWriteVector<Number>::local_element(const size_type local_index) const
939  {
940  AssertIndexRange(local_index, this->n_elements());
941 
942  return values[local_index];
943  }
944 
945 
946 
947  template <typename Number>
948  inline Number &
950  {
951  AssertIndexRange(local_index, this->n_elements());
952 
953  return values[local_index];
954  }
955 
956 
957 
958  template <typename Number>
959  template <typename Number2>
960  inline void
961  ReadWriteVector<Number>::add(const std::vector<size_type> &indices,
962  const std::vector<Number2> & values)
963  {
964  AssertDimension(indices.size(), values.size());
965  add(indices.size(), indices.data(), values.data());
966  }
967 
968 
969 
970  template <typename Number>
971  template <typename Number2>
972  inline void
973  ReadWriteVector<Number>::add(const std::vector<size_type> & indices,
974  const ReadWriteVector<Number2> &values)
975  {
976  const size_type size = indices.size();
977  for (size_type i = 0; i < size; ++i)
978  {
979  Assert(
980  numbers::is_finite(values[i]),
981  ExcMessage(
982  "The given value is not finite but either infinite or Not A Number (NaN)"));
983  this->operator()(indices[i]) += values[indices[i]];
984  }
985  }
986 
987 
988 
989  template <typename Number>
990  template <typename Number2>
991  inline void
992  ReadWriteVector<Number>::add(const size_type n_indices,
993  const size_type *indices,
994  const Number2 * values_to_add)
995  {
996  for (size_type i = 0; i < n_indices; ++i)
997  {
998  Assert(
999  numbers::is_finite(values[i]),
1000  ExcMessage(
1001  "The given value is not finite but either infinite or Not A Number (NaN)"));
1002  this->operator()(indices[i]) += values_to_add[i];
1003  }
1004  }
1005 
1006 
1007 
1008  template <typename Number>
1009  template <typename Functor>
1012  const Functor & functor)
1013  : parent(parent)
1014  , functor(functor)
1015  {}
1016 
1017 
1018 
1019  template <typename Number>
1020  template <typename Functor>
1021  void
1023  operator()(const size_type begin, const size_type end)
1024  {
1025  for (size_type i = begin; i < end; ++i)
1026  functor(parent.values[i]);
1027  }
1028 
1029 #endif // ifndef DOXYGEN
1030 
1031 } // end of namespace LinearAlgebra
1032 
1033 
1034 
1042 template <typename Number>
1043 inline void
1046 {
1047  u.swap(v);
1048 }
1049 
1050 
1052 
1053 #endif
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1580
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
#define AssertIndexRange(index, range)
Definition: exceptions.h:1648
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:1632
static ::ExceptionBase & ExcMessage(std::string arg1)
void swap(LinearAlgebra::ReadWriteVector< Number > &u, LinearAlgebra::ReadWriteVector< Number > &v)
#define Assert(cond, exc)
Definition: exceptions.h:1423
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
Number operator()(const size_type global_index) const
ReadWriteVector< Number > & operator=(const ReadWriteVector< Number > &in_vector)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
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)
std::shared_ptr< CommunicationPatternBase > comm_pattern
unsigned int global_dof_index
Definition: types.h:76
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
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:1830
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