Reference documentation for deal.II version GIT 29f9da0a34 2023-12-07 10:00:01+00:00
\(\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\}}\)
trilinos_tpetra_vector.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2023 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_trilinos_tpetra_vector_h
17 #define dealii_trilinos_tpetra_vector_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #ifdef DEAL_II_TRILINOS_WITH_TPETRA
23 
24 # include <deal.II/base/index_set.h>
25 # include <deal.II/base/mpi_stub.h>
27 
28 # include <deal.II/lac/read_vector.h>
32 
33 # include <Teuchos_Comm.hpp>
34 # include <Teuchos_OrdinalTraits.hpp>
35 # include <Tpetra_Core.hpp>
36 # include <Tpetra_Vector.hpp>
37 # include <Tpetra_Version.hpp>
38 
39 # include <memory>
40 
42 
48 template <typename Number>
49 struct is_tpetra_type : std::false_type
50 {};
51 
52 # ifdef HAVE_TPETRA_INST_FLOAT
53 template <>
54 struct is_tpetra_type<float> : std::true_type
55 {};
56 # endif
57 
58 # ifdef HAVE_TPETRA_INST_DOUBLE
59 template <>
60 struct is_tpetra_type<double> : std::true_type
61 {};
62 # endif
63 
64 # ifdef DEAL_II_WITH_COMPLEX_VALUES
65 # ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
66 template <>
67 struct is_tpetra_type<std::complex<float>> : std::true_type
68 {};
69 # endif
70 
71 # ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE
72 template <>
73 struct is_tpetra_type<std::complex<double>> : std::true_type
74 {};
75 # endif
76 # endif
77 
78 namespace LinearAlgebra
79 {
80  // Forward declaration
81 # ifndef DOXYGEN
82  template <typename Number>
83  class ReadWriteVector;
84 # endif
85 
100  namespace TpetraWrappers
101  {
112  namespace internal
113  {
118 
128  template <typename Number>
129  class VectorReference
130  {
131  private:
136  VectorReference(Vector<Number> &vector, const size_type index);
137 
138  public:
142  VectorReference(const VectorReference &) = default;
143 
155  const VectorReference &
156  operator=(const VectorReference &r) const;
157 
161  VectorReference &
162  operator=(const VectorReference &r);
163 
167  const VectorReference &
168  operator=(const Number &s) const;
169 
173  const VectorReference &
174  operator+=(const Number &s) const;
175 
179  const VectorReference &
180  operator-=(const Number &s) const;
181 
185  const VectorReference &
186  operator*=(const Number &s) const;
187 
191  const VectorReference &
192  operator/=(const Number &s) const;
193 
198  operator Number() const;
199 
204  int,
205  << "An error with error number " << arg1
206  << " occurred while calling a Trilinos function");
207 
208  private:
212  Vector<Number> &vector;
213 
217  const size_type index;
218 
219  // Make the vector class a friend, so that it can create objects of the
220  // present type.
221  friend class Vector<Number>;
222  }; // class VectorReference
223 
224  } // namespace internal
249  template <typename Number>
250  class Vector : public ReadVector<Number>, public Subscriptor
251  {
252  public:
256  using value_type = Number;
259  using reference = internal::VectorReference<Number>;
260  using MapType = Tpetra::Map<int, ::types::signed_global_dof_index>;
261  using VectorType =
262  Tpetra::Vector<Number, int, ::types::signed_global_dof_index>;
263 
274 
279  Vector(const Vector &V);
280 
284  Vector(const Teuchos::RCP<VectorType> V);
285 
292  explicit Vector(const IndexSet &parallel_partitioner,
293  const MPI_Comm communicator);
294 
301  void
302  reinit(const IndexSet &parallel_partitioner,
303  const MPI_Comm communicator,
304  const bool omit_zeroing_entries = false);
305 
306 
321  void
322  reinit(const IndexSet &locally_owned_entries,
323  const IndexSet &locally_relevant_or_ghost_entries,
324  const MPI_Comm communicator = MPI_COMM_WORLD);
325 
330  void
331  reinit(const Vector<Number> &V, const bool omit_zeroing_entries = false);
332 
336  virtual void
339  ArrayView<Number> &elements) const override;
340 
346  Vector &
347  operator=(const Vector &V);
348 
353  Vector &
354  operator=(const Number s);
355 
364  void
366  const ReadWriteVector<Number> &V,
367  VectorOperation::values operation,
368  const Teuchos::RCP<const Utilities::MPI::CommunicationPatternBase>
369  &communication_pattern);
370 
375  void
377  const ReadWriteVector<Number> &V,
378  VectorOperation::values operation,
379  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
380  &communication_pattern);
381 
382  /*
383  * Imports all the elements present in the vector's IndexSet from the
384  * input vector @p V. VectorOperation::values @p operation is used to decide if
385  * the elements in @p V should be added to the current vector or replace the
386  * current elements.
387  */
388  void
390  VectorOperation::values operation);
391 
396  void
397  import(const ReadWriteVector<Number> &V,
398  VectorOperation::values operation,
399  std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
400  communication_pattern = {})
401  {
402  import_elements(V, operation, communication_pattern);
403  }
404 
420  reference
421  operator()(const size_type index);
422 
434  Vector &
435  operator*=(const Number factor);
436 
440  Vector &
441  operator/=(const Number factor);
442 
446  Vector &
448 
452  Vector &
454 
459  Number
460  operator*(const Vector<Number> &V) const;
461 
465  void
466  add(const Number a);
467 
472  void
473  add(const Number a, const Vector<Number> &V);
474 
479  void
480  add(const Number a,
481  const Vector<Number> &V,
482  const Number b,
483  const Vector<Number> &W);
484 
489  void
490  add(const std::vector<size_type> &indices,
491  const std::vector<TrilinosScalar> &values);
492 
493 
499  void
500  add(const size_type n_elements,
501  const size_type *indices,
502  const Number *values);
503 
508  void
509  sadd(const Number s, const Number a, const Vector<Number> &V);
510 
517  void
518  set(const size_type n_elements,
519  const size_type *indices,
520  const Number *values);
521 
528  void
529  scale(const Vector<Number> &scaling_factors);
530 
534  void
535  equ(const Number a, const Vector<Number> &V);
536 
540  bool
541  all_zero() const;
542 
554  Number
555  mean_value() const;
556 
561  real_type
562  l1_norm() const;
563 
568  real_type
569  l2_norm() const;
570 
575  real_type
576  linfty_norm() const;
577 
600  Number
601  add_and_dot(const Number a,
602  const Vector<Number> &V,
603  const Vector<Number> &W);
604 
617  bool
618  has_ghost_elements() const;
619 
624  virtual size_type
625  size() const override;
626 
631  size_type
633 
637  MPI_Comm
639 
651  ::IndexSet
653 
671  void
673 
678  const Tpetra::Vector<Number, int, types::signed_global_dof_index> &
680 
685  Tpetra::Vector<Number, int, types::signed_global_dof_index> &
687 
692  Teuchos::RCP<
693  const Tpetra::Vector<Number, int, types::signed_global_dof_index>>
694  trilinos_rcp() const;
695 
700  Teuchos::RCP<Tpetra::Vector<Number, int, types::signed_global_dof_index>>
702 
706  void
707  print(std::ostream &out,
708  const unsigned int precision = 3,
709  const bool scientific = true,
710  const bool across = true) const;
711 
715  std::size_t
717 
721  MPI_Comm
722  mpi_comm() const;
723 
734 
741 
748  int,
749  << "An error with error number " << arg1
750  << " occurred while calling a Trilinos function");
751 
752  private:
758  void
759  create_tpetra_comm_pattern(const IndexSet &source_index_set,
760  const MPI_Comm mpi_comm);
761 
765  Teuchos::RCP<VectorType> vector;
766 
771 
776  Teuchos::RCP<const TpetraWrappers::CommunicationPattern>
778 
779  // Make the reference class a friend.
780  friend class internal::VectorReference<Number>;
781  };
782 
783 
784  /* ------------------------- Inline functions ---------------------- */
785 
786  template <typename Number>
787  inline bool
789  {
790  return false;
791  }
792 
793 
794 
795  template <typename Number>
796  inline void
797  Vector<Number>::add(const std::vector<size_type> &indices,
798  const std::vector<TrilinosScalar> &values)
799  {
800  // if we have ghost values, do not allow
801  // writing to this vector at all.
802  AssertDimension(indices.size(), values.size());
803 
804  add(indices.size(), indices.data(), values.data());
805  }
806 
807 
808 
809  template <typename Number>
810  inline void
811  Vector<Number>::add(const size_type n_elements,
812  const size_type *indices,
813  const Number *values)
814  {
815 # if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
816  auto vector_2d = vector->template getLocalView<Kokkos::HostSpace>(
817  Tpetra::Access::ReadWrite);
818 # else
819  vector->template sync<Kokkos::HostSpace>();
820  auto vector_2d = vector->template getLocalView<Kokkos::HostSpace>();
821 # endif
822  auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0);
823 # if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
824  vector->template modify<Kokkos::HostSpace>();
825 # endif
826 
827  for (size_type i = 0; i < n_elements; ++i)
828  {
829  const size_type row = indices[i];
830  const TrilinosWrappers::types::int_type local_row =
831  vector->getMap()->getLocalElement(row);
832 
833  if (local_row != Teuchos::OrdinalTraits<int>::invalid())
834  vector_1d(local_row) += values[i];
835  }
836 
837 # if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
838  vector->template sync<
839  typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
840  device_type::memory_space>();
841 # endif
842  }
843 
844 
845 
846  template <typename Number>
847  inline void
848  Vector<Number>::set(const size_type n_elements,
849  const size_type *indices,
850  const Number *values)
851  {
852 # if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
853  auto vector_2d = vector->template getLocalView<Kokkos::HostSpace>(
854  Tpetra::Access::ReadWrite);
855 # else
856  vector->template sync<Kokkos::HostSpace>();
857  auto vector_2d = vector->template getLocalView<Kokkos::HostSpace>();
858 # endif
859  auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0);
860 # if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
861  vector->template modify<Kokkos::HostSpace>();
862 # endif
863 
864  for (size_type i = 0; i < n_elements; ++i)
865  {
866  const size_type row = indices[i];
867  const TrilinosWrappers::types::int_type local_row =
868  vector->getMap()->getLocalElement(row);
869 
870  if (local_row != Teuchos::OrdinalTraits<int>::invalid())
871  vector_1d(local_row) = values[i];
872  }
873 
874 # if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
875  vector->template sync<
876  typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
877  device_type::memory_space>();
878 # endif
879  }
880 
881 
882 
883  template <typename Number>
884  inline internal::VectorReference<Number>
886  {
887  return internal::VectorReference(*this, index);
888  }
889 
890 # ifndef DOXYGEN
891 
892  // VectorReference
893  namespace internal
894  {
895  template <typename Number>
896  inline VectorReference<Number>::VectorReference(Vector<Number> &vector,
897  const size_type index)
898  : vector(vector)
899  , index(index)
900  {}
901 
902 
903 
904  template <typename Number>
905  inline const VectorReference<Number> &
906  VectorReference<Number>::operator=(const VectorReference<Number> &r) const
907  {
908  // as explained in the class
909  // documentation, this is not the copy
910  // operator. so simply pass on to the
911  // "correct" assignment operator
912  *this = static_cast<Number>(r);
913 
914  return *this;
915  }
916 
917 
918 
919  template <typename Number>
920  inline VectorReference<Number> &
921  VectorReference<Number>::operator=(const VectorReference<Number> &r)
922  {
923  // as above
924  *this = static_cast<Number>(r);
925 
926  return *this;
927  }
928 
929 
930 
931  template <typename Number>
932  inline const VectorReference<Number> &
933  VectorReference<Number>::operator=(const Number &value) const
934  {
935  vector.set(1, &index, &value);
936  return *this;
937  }
938 
939 
940 
941  template <typename Number>
942  inline const VectorReference<Number> &
943  VectorReference<Number>::operator+=(const Number &value) const
944  {
945  vector.add(1, &index, &value);
946  return *this;
947  }
948 
949 
950 
951  template <typename Number>
952  inline const VectorReference<Number> &
953  VectorReference<Number>::operator-=(const Number &value) const
954  {
955  Number new_value = -value;
956  vector.add(1, &index, &new_value);
957  return *this;
958  }
959 
960 
961 
962  template <typename Number>
963  inline const VectorReference<Number> &
964  VectorReference<Number>::operator*=(const Number &value) const
965  {
966  Number new_value = static_cast<Number>(*this) * value;
967  vector.set(1, &index, &new_value);
968  return *this;
969  }
970 
971 
972 
973  template <typename Number>
974  inline const VectorReference<Number> &
975  VectorReference<Number>::operator/=(const Number &value) const
976  {
977  Number new_value = static_cast<Number>(*this) / value;
978  vector.set(1, &index, &new_value);
979  return *this;
980  }
981  } // namespace internal
982 
983 # endif /* DOXYGEN */
984 
985  } // namespace TpetraWrappers
986 
989 } // namespace LinearAlgebra
990 
994 template <typename Number>
996  : std::false_type
997 {};
998 
1000 
1001 #endif
1002 
1003 #endif
Vector & operator-=(const Vector< Number > &V)
Number operator*(const Vector< Number > &V) const
void create_tpetra_comm_pattern(const IndexSet &source_index_set, const MPI_Comm mpi_comm)
Teuchos::RCP< const TpetraWrappers::CommunicationPattern > tpetra_comm_pattern
Vector & operator=(const Vector &V)
virtual size_type size() const override
Teuchos::RCP< const Tpetra::Vector< Number, int, types::signed_global_dof_index > > trilinos_rcp() const
void reinit(const IndexSet &locally_owned_entries, const IndexSet &locally_relevant_or_ghost_entries, const MPI_Comm communicator=MPI_COMM_WORLD)
void add(const Number a, const Vector< Number > &V)
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
std::size_t memory_consumption() const
Vector & operator/=(const Number factor)
Vector & operator=(const Number s)
Vector & operator*=(const Number factor)
void compress(const VectorOperation::values operation)
void import_elements(const ReadWriteVector< Number > &V, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern)
Tpetra::Vector< Number, int, types::signed_global_dof_index > & trilinos_vector()
Tpetra::Vector< Number, int, ::types::signed_global_dof_index > VectorType
void equ(const Number a, const Vector< Number > &V)
Vector & operator+=(const Vector< Number > &V)
void import_elements(const ReadWriteVector< Number > &V, VectorOperation::values operation, const Teuchos::RCP< const Utilities::MPI::CommunicationPatternBase > &communication_pattern)
void sadd(const Number s, const Number a, const Vector< Number > &V)
Tpetra::Map< int, ::types::signed_global_dof_index > MapType
typename numbers::NumberTraits< Number >::real_type real_type
void reinit(const IndexSet &parallel_partitioner, const MPI_Comm communicator, const bool omit_zeroing_entries=false)
reference operator()(const size_type index)
virtual void extract_subvector_to(const ArrayView< const types::global_dof_index > &indices, ArrayView< Number > &elements) const override
Teuchos::RCP< Tpetra::Vector< Number, int, types::signed_global_dof_index > > trilinos_rcp()
void scale(const Vector< Number > &scaling_factors)
void import_elements(const ReadWriteVector< Number > &V, VectorOperation::values operation)
::IndexSet locally_owned_elements() const
const Tpetra::Vector< Number, int, types::signed_global_dof_index > & trilinos_vector() const
Number add_and_dot(const Number a, const Vector< Number > &V, const Vector< Number > &W)
void reinit(const Vector< Number > &V, const bool omit_zeroing_entries=false)
Vector(const IndexSet &parallel_partitioner, const MPI_Comm communicator)
void add(const Number a, const Vector< Number > &V, const Number b, const Vector< Number > &W)
void set(const size_type n_elements, const size_type *indices, const Number *values)
internal::VectorReference< Number > reference
Vector(const Teuchos::RCP< VectorType > V)
types::global_dof_index size_type
Definition: read_vector.h:44
Definition: vector.h:110
#define DEAL_II_DEPRECATED
Definition: config.h:178
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define DeclException0(Exception0)
Definition: exceptions.h:472
static ::ExceptionBase & ExcDifferentParallelPartitioning()
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:517
static ::ExceptionBase & ExcVectorTypeNotCompatible()
static const char V
types::global_dof_index size_type
Definition: cuda_kernels.h:45
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
unsigned int global_dof_index
Definition: types.h:82