Reference documentation for deal.II version Git 7e76c9f1da 2019-09-19 09:07:25 +0200
\(\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/thread_management.h>
26 #include <deal.II/base/types.h>
27 #include <deal.II/base/utilities.h>
28 
29 #include <deal.II/lac/vector_operation.h>
30 
31 #include <cstdlib>
32 #include <cstring>
33 #include <iomanip>
34 
35 #ifdef DEAL_II_WITH_TRILINOS
36 # include <deal.II/lac/trilinos_epetra_communication_pattern.h>
37 # include <deal.II/lac/trilinos_epetra_vector.h>
38 # include <deal.II/lac/trilinos_tpetra_vector.h>
39 
40 # include <Epetra_MultiVector.h>
41 
42 #endif
43 
44 DEAL_II_NAMESPACE_OPEN
45 
46 // Forward declarations
47 #ifndef DOXYGEN
48 namespace LinearAlgebra
49 {
50  class CommunicationPatternBase;
51  namespace distributed
52  {
53  template <typename, typename>
54  class Vector;
55  } // namespace distributed
56 } // namespace LinearAlgebra
57 
58 # ifdef DEAL_II_WITH_PETSC
59 namespace PETScWrappers
60 {
61  namespace MPI
62  {
63  class Vector;
64  }
65 } // namespace PETScWrappers
66 # endif
67 
68 # ifdef DEAL_II_WITH_TRILINOS
69 namespace TrilinosWrappers
70 {
71  namespace MPI
72  {
73  class Vector;
74  }
75 } // namespace TrilinosWrappers
76 # endif
77 
78 # ifdef DEAL_II_WITH_CUDA
79 namespace LinearAlgebra
80 {
81  namespace CUDAWrappers
82  {
83  template <typename>
84  class Vector;
85  }
86 } // namespace LinearAlgebra
87 # endif
88 #endif
89 
90 namespace LinearAlgebra
91 {
131  template <typename Number>
133  {
134  public:
140  using value_type = Number;
141  using pointer = value_type *;
142  using const_pointer = const value_type *;
143  using iterator = value_type *;
144  using const_iterator = const value_type *;
145  using reference = value_type &;
146  using const_reference = const value_type &;
147  using size_type = types::global_dof_index;
148  using real_type = typename numbers::NumberTraits<Number>::real_type;
149 
157  ReadWriteVector();
158 
162  ReadWriteVector(const ReadWriteVector<Number> &in_vector);
163 
168  explicit ReadWriteVector(const size_type size);
169 
174  explicit ReadWriteVector(const IndexSet &locally_stored_indices);
175 
179  ~ReadWriteVector() override = default;
180 
189  virtual void
190  reinit(const size_type size, const bool omit_zeroing_entries = false);
191 
200  template <typename Number2>
201  void
202  reinit(const ReadWriteVector<Number2> &in_vector,
203  const bool omit_zeroing_entries = false);
204 
214  virtual void
215  reinit(const IndexSet &locally_stored_indices,
216  const bool omit_zeroing_entries = false);
217 
218 
219 #ifdef DEAL_II_WITH_TRILINOS
220 # ifdef DEAL_II_WITH_MPI
221 
232  void
233  reinit(const TrilinosWrappers::MPI::Vector &trilinos_vec);
234 # endif
235 #endif
236 
250  template <typename Functor>
251  void
252  apply(const Functor &func);
253 
266  void
268 
273  operator=(const ReadWriteVector<Number> &in_vector);
274 
278  template <typename Number2>
280  operator=(const ReadWriteVector<Number2> &in_vector);
281 
287  operator=(const Number s);
288 
297  template <typename MemorySpace>
298  void
300  VectorOperation::values operation,
301  const std::shared_ptr<const CommunicationPatternBase>
302  &communication_pattern =
303  std::shared_ptr<const CommunicationPatternBase>());
304 
305 #ifdef DEAL_II_WITH_PETSC
306 
314  void
315  import(const PETScWrappers::MPI::Vector &petsc_vec,
316  VectorOperation::values operation,
317  const std::shared_ptr<const CommunicationPatternBase>
318  &communication_pattern =
319  std::shared_ptr<const CommunicationPatternBase>());
320 #endif
321 
322 #ifdef DEAL_II_WITH_TRILINOS
323 
333  void
334  import(const TrilinosWrappers::MPI::Vector &trilinos_vec,
335  VectorOperation::values operation,
336  const std::shared_ptr<const CommunicationPatternBase>
337  &communication_pattern =
338  std::shared_ptr<const CommunicationPatternBase>());
339 
340 # ifdef DEAL_II_WITH_MPI
341 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
342 
350  void
351  import(const TpetraWrappers::Vector<Number> &tpetra_vec,
352  VectorOperation::values operation,
353  const std::shared_ptr<const CommunicationPatternBase>
354  &communication_pattern =
355  std::shared_ptr<const CommunicationPatternBase>());
356 # endif
357 
366  void
367  import(const EpetraWrappers::Vector &epetra_vec,
368  VectorOperation::values operation,
369  const std::shared_ptr<const CommunicationPatternBase>
370  &communication_pattern =
371  std::shared_ptr<const CommunicationPatternBase>());
372 # endif
373 #endif
374 
375 #ifdef DEAL_II_WITH_CUDA
376 
382  void
383  import(const CUDAWrappers::Vector<Number> &cuda_vec,
384  VectorOperation::values operation,
385  const std::shared_ptr<const CommunicationPatternBase>
386  &communication_pattern =
387  std::shared_ptr<const CommunicationPatternBase>());
388 #endif
389 
398  size_type
399  size() const;
400 
406  size_type
407  n_elements() const;
408 
412  const IndexSet &
413  get_stored_elements() const;
414 
420  iterator
421  begin();
422 
427  const_iterator
428  begin() const;
429 
434  iterator
435  end();
436 
441  const_iterator
442  end() const;
444 
445 
450 
456  Number
457  operator()(const size_type global_index) const;
458 
464  Number &
465  operator()(const size_type global_index);
466 
474  Number operator[](const size_type global_index) const;
475 
483  Number &operator[](const size_type global_index);
484 
500  template <typename Number2>
501  void
502  extract_subvector_to(const std::vector<size_type> &indices,
503  std::vector<Number2> & values) const;
504 
532  template <typename ForwardIterator, typename OutputIterator>
533  void
534  extract_subvector_to(ForwardIterator indices_begin,
535  const ForwardIterator indices_end,
536  OutputIterator values_begin) const;
537 
548  Number
549  local_element(const size_type local_index) const;
550 
561  Number &
562  local_element(const size_type local_index);
564 
565 
570 
575  template <typename Number2>
576  void
577  add(const std::vector<size_type> &indices,
578  const std::vector<Number2> & values);
579 
584  template <typename Number2>
585  void
586  add(const std::vector<size_type> & indices,
587  const ReadWriteVector<Number2> &values);
588 
594  template <typename Number2>
595  void
596  add(const size_type n_elements,
597  const size_type *indices,
598  const Number2 * values);
599 
603  void
604  print(std::ostream & out,
605  const unsigned int precision = 3,
606  const bool scientific = true) const;
607 
611  std::size_t
612  memory_consumption() const;
614 
615  protected:
616 #ifdef DEAL_II_WITH_TRILINOS
617 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
618 
623  void
624  import(
625  const Tpetra::Vector<Number, int, types::global_dof_index> &tpetra_vector,
626  const IndexSet & locally_owned_elements,
627  VectorOperation::values operation,
628  const MPI_Comm & mpi_comm,
629  const std::shared_ptr<const CommunicationPatternBase>
630  &communication_pattern);
631 # endif
632 
638  void
639  import(const Epetra_MultiVector &multivector,
640  const IndexSet & locally_owned_elements,
641  VectorOperation::values operation,
642  const MPI_Comm & mpi_comm,
643  const std::shared_ptr<const CommunicationPatternBase>
644  &communication_pattern);
645 #endif
646 
650  unsigned int
651  global_to_local(const types::global_dof_index global_index) const
652  {
653  // the following will throw an exception if the global_index is not
654  // in the remaining_elements
655  return static_cast<unsigned int>(
656  stored_elements.index_within_set(global_index));
657  }
658 
662  void
663  resize_val(const size_type new_allocated_size);
664 
665 #if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI)
666 # ifdef DEAL_II_TRILINOS_WITH_TPETRA
667 
671  TpetraWrappers::CommunicationPattern
672  create_tpetra_comm_pattern(const IndexSet &source_index_set,
673  const MPI_Comm &mpi_comm);
674 # endif
675 
681  create_epetra_comm_pattern(const IndexSet &source_index_set,
682  const MPI_Comm &mpi_comm);
683 #endif
684 
689 
694 
699  std::shared_ptr<CommunicationPatternBase> comm_pattern;
700 
704  std::unique_ptr<Number[], decltype(std::free) *> values;
705 
710  mutable std::shared_ptr<::parallel::internal::TBBPartitioner>
712 
713  // Make all other ReadWriteVector types friends.
714  template <typename Number2>
715  friend class ReadWriteVector;
716 
717  private:
723  template <typename Functor>
725  {
726  public:
730  FunctorTemplate(ReadWriteVector<Number> &parent, const Functor &functor);
731 
735  virtual void
736  operator()(const size_type begin, const size_type end);
737 
738  private:
743 
747  const Functor &functor;
748  };
749  };
750 
754  /*---------------------------- Inline functions ---------------------------*/
755 
756 #ifndef DOXYGEN
757 
758  template <typename Number>
760  : Subscriptor()
761  , values(nullptr, free)
762  {
763  // virtual functions called in constructors and destructors never use the
764  // override in a derived class
765  // for clarity be explicit on which function is called
767  }
768 
769 
770 
771  template <typename Number>
773  const ReadWriteVector<Number> &v)
774  : Subscriptor()
775  , values(nullptr, free)
776  {
777  this->operator=(v);
778  }
779 
780 
781 
782  template <typename Number>
783  inline ReadWriteVector<Number>::ReadWriteVector(const size_type size)
784  : Subscriptor()
785  , values(nullptr, free)
786  {
787  // virtual functions called in constructors and destructors never use the
788  // override in a derived class
789  // for clarity be explicit on which function is called
791  }
792 
793 
794 
795  template <typename Number>
797  const IndexSet &locally_stored_indices)
798  : Subscriptor()
799  , values(nullptr, free)
800  {
801  // virtual functions called in constructors and destructors never use the
802  // override in a derived class
803  // for clarity be explicit on which function is called
804  ReadWriteVector<Number>::reinit(locally_stored_indices);
805  }
806 
807 
808 
809  template <typename Number>
810  inline typename ReadWriteVector<Number>::size_type
812  {
813  return stored_elements.size();
814  }
815 
816 
817 
818  template <typename Number>
819  inline typename ReadWriteVector<Number>::size_type
821  {
822  return stored_elements.n_elements();
823  }
824 
825 
826 
827  template <typename Number>
828  inline const IndexSet &
830  {
831  return stored_elements;
832  }
833 
834 
835 
836  template <typename Number>
837  inline typename ReadWriteVector<Number>::iterator
839  {
840  return values.get();
841  }
842 
843 
844 
845  template <typename Number>
846  inline typename ReadWriteVector<Number>::const_iterator
848  {
849  return values.get();
850  }
851 
852 
853 
854  template <typename Number>
855  inline typename ReadWriteVector<Number>::iterator
857  {
858  return values.get() + this->n_elements();
859  }
860 
861 
862 
863  template <typename Number>
864  inline typename ReadWriteVector<Number>::const_iterator
866  {
867  return values.get() + this->n_elements();
868  }
869 
870 
871 
872  template <typename Number>
873  inline Number
874  ReadWriteVector<Number>::operator()(const size_type global_index) const
875  {
876  return values[global_to_local(global_index)];
877  }
878 
879 
880 
881  template <typename Number>
882  inline Number &
883  ReadWriteVector<Number>::operator()(const size_type global_index)
884  {
885  return values[global_to_local(global_index)];
886  }
887 
888 
889 
890  template <typename Number>
891  inline Number ReadWriteVector<Number>::
892  operator[](const size_type global_index) const
893  {
894  return operator()(global_index);
895  }
896 
897 
898 
899  template <typename Number>
900  inline Number &ReadWriteVector<Number>::
901  operator[](const size_type global_index)
902  {
903  return operator()(global_index);
904  }
905 
906 
907 
908  template <typename Number>
909  template <typename Number2>
910  inline void
912  const std::vector<size_type> &indices,
913  std::vector<Number2> & extracted_values) const
914  {
915  for (size_type i = 0; i < indices.size(); ++i)
916  extracted_values[i] = operator()(indices[i]);
917  }
918 
919 
920 
921  template <typename Number>
922  template <typename ForwardIterator, typename OutputIterator>
923  inline void
925  ForwardIterator indices_begin,
926  const ForwardIterator indices_end,
927  OutputIterator values_begin) const
928  {
929  while (indices_begin != indices_end)
930  {
931  *values_begin = operator()(*indices_begin);
932  indices_begin++;
933  values_begin++;
934  }
935  }
936 
937 
938 
939  template <typename Number>
940  inline Number
941  ReadWriteVector<Number>::local_element(const size_type local_index) const
942  {
943  AssertIndexRange(local_index, this->n_elements());
944 
945  return values[local_index];
946  }
947 
948 
949 
950  template <typename Number>
951  inline Number &
952  ReadWriteVector<Number>::local_element(const size_type local_index)
953  {
954  AssertIndexRange(local_index, this->n_elements());
955 
956  return values[local_index];
957  }
958 
959 
960 
961  template <typename Number>
962  template <typename Number2>
963  inline void
964  ReadWriteVector<Number>::add(const std::vector<size_type> &indices,
965  const std::vector<Number2> & values)
966  {
967  AssertDimension(indices.size(), values.size());
968  add(indices.size(), indices.data(), values.data());
969  }
970 
971 
972 
973  template <typename Number>
974  template <typename Number2>
975  inline void
976  ReadWriteVector<Number>::add(const std::vector<size_type> & indices,
977  const ReadWriteVector<Number2> &values)
978  {
979  const size_type size = indices.size();
980  for (size_type i = 0; i < size; ++i)
981  {
982  Assert(
983  numbers::is_finite(values[i]),
984  ExcMessage(
985  "The given value is not finite but either infinite or Not A Number (NaN)"));
986  this->operator()(indices[i]) += values[indices[i]];
987  }
988  }
989 
990 
991 
992  template <typename Number>
993  template <typename Number2>
994  inline void
995  ReadWriteVector<Number>::add(const size_type n_indices,
996  const size_type *indices,
997  const Number2 * values_to_add)
998  {
999  for (size_type i = 0; i < n_indices; ++i)
1000  {
1001  Assert(
1002  numbers::is_finite(values[i]),
1003  ExcMessage(
1004  "The given value is not finite but either infinite or Not A Number (NaN)"));
1005  this->operator()(indices[i]) += values_to_add[i];
1006  }
1007  }
1008 
1009 
1010 
1011  template <typename Number>
1012  template <typename Functor>
1015  const Functor & functor)
1016  : parent(parent)
1017  , functor(functor)
1018  {}
1019 
1020 
1021 
1022  template <typename Number>
1023  template <typename Functor>
1024  void
1026  operator()(const size_type begin, const size_type end)
1027  {
1028  for (size_type i = begin; i < end; ++i)
1029  functor(parent.values[i]);
1030  }
1031 
1032 #endif // ifndef DOXYGEN
1033 
1034 } // end of namespace LinearAlgebra
1035 
1036 
1037 
1045 template <typename Number>
1046 inline void
1049 {
1050  u.swap(v);
1051 }
1052 
1053 
1054 DEAL_II_NAMESPACE_CLOSE
1055 
1056 #endif
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
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:1637
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:1600
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
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:1732
const IndexSet & get_stored_elements() const
std::shared_ptr<::parallel::internal::TBBPartitioner > thread_loop_partitioner