Reference documentation for deal.II version Git 8236ec9b21 2019-11-20 14:24:58 -0500
\(\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/types.h>
26 #include <deal.II/base/utilities.h>
27 
28 #include <deal.II/lac/vector_operation.h>
29 
30 #include <cstdlib>
31 #include <cstring>
32 #include <iomanip>
33 
34 #ifdef DEAL_II_WITH_TRILINOS
35 # include <deal.II/lac/trilinos_epetra_communication_pattern.h>
36 # include <deal.II/lac/trilinos_epetra_vector.h>
37 # include <deal.II/lac/trilinos_tpetra_vector.h>
38 
39 # include <Epetra_MultiVector.h>
40 
41 #endif
42 
43 DEAL_II_NAMESPACE_OPEN
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 {
130  template <typename Number>
132  {
133  public:
139  using value_type = Number;
140  using pointer = value_type *;
141  using const_pointer = const value_type *;
142  using iterator = value_type *;
143  using const_iterator = const value_type *;
144  using reference = value_type &;
145  using const_reference = const value_type &;
146  using size_type = types::global_dof_index;
147  using real_type = typename numbers::NumberTraits<Number>::real_type;
148 
156  ReadWriteVector();
157 
161  ReadWriteVector(const ReadWriteVector<Number> &in_vector);
162 
167  explicit ReadWriteVector(const size_type size);
168 
173  explicit ReadWriteVector(const IndexSet &locally_stored_indices);
174 
178  ~ReadWriteVector() override = default;
179 
188  virtual void
189  reinit(const size_type size, const bool omit_zeroing_entries = false);
190 
199  template <typename Number2>
200  void
201  reinit(const ReadWriteVector<Number2> &in_vector,
202  const bool omit_zeroing_entries = false);
203 
213  virtual void
214  reinit(const IndexSet &locally_stored_indices,
215  const bool omit_zeroing_entries = false);
216 
217 
218 #ifdef DEAL_II_WITH_TRILINOS
219 # ifdef DEAL_II_WITH_MPI
220 
231  void
232  reinit(const TrilinosWrappers::MPI::Vector &trilinos_vec);
233 # endif
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 
296  template <typename MemorySpace>
297  void
299  VectorOperation::values operation,
300  const std::shared_ptr<const CommunicationPatternBase>
301  &communication_pattern =
302  std::shared_ptr<const CommunicationPatternBase>());
303 
304 #ifdef DEAL_II_WITH_PETSC
305 
313  void
314  import(const PETScWrappers::MPI::Vector &petsc_vec,
315  VectorOperation::values operation,
316  const std::shared_ptr<const CommunicationPatternBase>
317  &communication_pattern =
318  std::shared_ptr<const CommunicationPatternBase>());
319 #endif
320 
321 #ifdef DEAL_II_WITH_TRILINOS
322 
332  void
333  import(const TrilinosWrappers::MPI::Vector &trilinos_vec,
334  VectorOperation::values operation,
335  const std::shared_ptr<const CommunicationPatternBase>
336  &communication_pattern =
337  std::shared_ptr<const CommunicationPatternBase>());
338 
339 # ifdef DEAL_II_WITH_MPI
340 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
341 
349  void
350  import(const TpetraWrappers::Vector<Number> &tpetra_vec,
351  VectorOperation::values operation,
352  const std::shared_ptr<const CommunicationPatternBase>
353  &communication_pattern =
354  std::shared_ptr<const CommunicationPatternBase>());
355 # endif
356 
365  void
366  import(const EpetraWrappers::Vector &epetra_vec,
367  VectorOperation::values operation,
368  const std::shared_ptr<const CommunicationPatternBase>
369  &communication_pattern =
370  std::shared_ptr<const CommunicationPatternBase>());
371 # endif
372 #endif
373 
374 #ifdef DEAL_II_WITH_CUDA
375 
381  void
382  import(const CUDAWrappers::Vector<Number> &cuda_vec,
383  VectorOperation::values operation,
384  const std::shared_ptr<const CommunicationPatternBase>
385  &communication_pattern =
386  std::shared_ptr<const CommunicationPatternBase>());
387 #endif
388 
397  size_type
398  size() const;
399 
405  size_type
406  n_elements() const;
407 
411  const IndexSet &
412  get_stored_elements() const;
413 
419  iterator
420  begin();
421 
426  const_iterator
427  begin() const;
428 
433  iterator
434  end();
435 
440  const_iterator
441  end() const;
443 
444 
449 
455  Number
456  operator()(const size_type global_index) const;
457 
463  Number &
464  operator()(const size_type global_index);
465 
473  Number operator[](const size_type global_index) const;
474 
482  Number &operator[](const size_type global_index);
483 
499  template <typename Number2>
500  void
501  extract_subvector_to(const std::vector<size_type> &indices,
502  std::vector<Number2> & values) const;
503 
531  template <typename ForwardIterator, typename OutputIterator>
532  void
533  extract_subvector_to(ForwardIterator indices_begin,
534  const ForwardIterator indices_end,
535  OutputIterator values_begin) const;
536 
547  Number
548  local_element(const size_type local_index) const;
549 
560  Number &
561  local_element(const size_type local_index);
563 
564 
569 
574  template <typename Number2>
575  void
576  add(const std::vector<size_type> &indices,
577  const std::vector<Number2> & values);
578 
583  template <typename Number2>
584  void
585  add(const std::vector<size_type> & indices,
586  const ReadWriteVector<Number2> &values);
587 
593  template <typename Number2>
594  void
595  add(const size_type n_elements,
596  const size_type *indices,
597  const Number2 * values);
598 
602  void
603  print(std::ostream & out,
604  const unsigned int precision = 3,
605  const bool scientific = true) const;
606 
610  std::size_t
611  memory_consumption() const;
613 
614  protected:
615 #ifdef DEAL_II_WITH_TRILINOS
616 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
617 
622  void
623  import(
624  const Tpetra::Vector<Number, int, types::global_dof_index> &tpetra_vector,
625  const IndexSet & locally_owned_elements,
626  VectorOperation::values operation,
627  const MPI_Comm & mpi_comm,
628  const std::shared_ptr<const CommunicationPatternBase>
629  &communication_pattern);
630 # endif
631 
637  void
638  import(const Epetra_MultiVector &multivector,
639  const IndexSet & locally_owned_elements,
640  VectorOperation::values operation,
641  const MPI_Comm & mpi_comm,
642  const std::shared_ptr<const CommunicationPatternBase>
643  &communication_pattern);
644 #endif
645 
649  unsigned int
650  global_to_local(const types::global_dof_index global_index) const
651  {
652  // the following will throw an exception if the global_index is not
653  // in the remaining_elements
654  return static_cast<unsigned int>(
655  stored_elements.index_within_set(global_index));
656  }
657 
661  void
662  resize_val(const size_type new_allocated_size);
663 
664 #if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI)
665 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
666 
670  TpetraWrappers::CommunicationPattern
671  create_tpetra_comm_pattern(const IndexSet &source_index_set,
672  const MPI_Comm &mpi_comm);
673 # endif
674 
680  create_epetra_comm_pattern(const IndexSet &source_index_set,
681  const MPI_Comm &mpi_comm);
682 #endif
683 
688 
693 
698  std::shared_ptr<CommunicationPatternBase> comm_pattern;
699 
703  std::unique_ptr<Number[], decltype(std::free) *> values;
704 
709  mutable std::shared_ptr<::parallel::internal::TBBPartitioner>
711 
712  // Make all other ReadWriteVector types friends.
713  template <typename Number2>
714  friend class ReadWriteVector;
715 
716  private:
722  template <typename Functor>
724  {
725  public:
729  FunctorTemplate(ReadWriteVector<Number> &parent, const Functor &functor);
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(
982  numbers::is_finite(values[i]),
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(
1001  numbers::is_finite(values[i]),
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>
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:1571
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:1641
void swap(LinearAlgebra::ReadWriteVector< Number > &u, LinearAlgebra::ReadWriteVector< Number > &v)
LinearAlgebra::distributed::Vector< Number > Vector
bool is_finite(const double x)
Definition: numbers.h:555
virtual void reinit(const size_type size, const bool omit_zeroing_entries=false)
size_type size() const
Definition: index_set.h:1607
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1411
Number operator()(const size_type global_index) const
ReadWriteVector< Number > & operator=(const ReadWriteVector< Number > &in_vector)
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
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
unsigned int global_dof_index
Definition: types.h:89
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< Number2 > &values) const
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:1806
const IndexSet & get_stored_elements() const
std::shared_ptr<::parallel::internal::TBBPartitioner > thread_loop_partitioner