trilinos_epetra_vector.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2016 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_trilinos_epetra_vector_h
#define dealii_trilinos_epetra_vector_h
 
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_WITH_TRILINOS
 
#  include <deal.II/base/index_set.h>
#  include <deal.II/base/mpi_stub.h>
#  include <deal.II/base/subscriptor.h>
 
#  include <deal.II/lac/read_vector.h>
#  include <deal.II/lac/trilinos_epetra_communication_pattern.h>
#  include <deal.II/lac/vector_operation.h>
#  include <deal.II/lac/vector_type_traits.h>
 
#  include <Epetra_FEVector.h>
 
#  include <memory>
 
DEAL_II_NAMESPACE_OPEN
 
namespace LinearAlgebra
{
  // Forward declaration
  template <typename Number>
  class ReadWriteVector;
 
  namespace EpetraWrappers
  {
    class VectorTraits
    {
    public:
      using value_type = double;
      using size_type  = ::types::global_dof_index;
    };
 
    namespace internal
    {
      class VectorReference
      {
      private:
        using value_type = VectorTraits::value_type;
        using size_type  = VectorTraits::size_type;
 
        VectorReference(Vector &vector, const size_type index);
 
      public:
        VectorReference(const VectorReference &) = default;
 
        const VectorReference &
        operator=(const VectorReference &r) const;
 
        VectorReference &
        operator=(const VectorReference &r);
 
        const VectorReference &
        operator=(const value_type &s) const;
 
        const VectorReference &
        operator+=(const value_type &s) const;
 
        const VectorReference &
        operator-=(const value_type &s) const;
 
        const VectorReference &
        operator*=(const value_type &s) const;
 
        const VectorReference &
        operator/=(const value_type &s) const;
 
        operator value_type() const;
 
        DeclException1(ExcTrilinosError,
                       int,
                       << "An error with error number " << arg1
                       << " occurred while calling a Trilinos function");
 
        /*
         * Access to a an element that is not (locally-)owned.
         *
         * @ingroup Exceptions
         */
        DeclException4(
          ExcAccessToNonLocalElement,
          size_type,
          size_type,
          size_type,
          size_type,
          << "You are trying to access element " << arg1
          << " of a distributed vector, but this element is not stored "
          << "on the current processor. Note: There are " << arg2
          << " elements stored "
          << "on the current processor from within the range [" << arg3 << ','
          << arg4 << "] but Trilinos vectors need not store contiguous "
          << "ranges on each processor, and not every element in "
          << "this range may in fact be stored locally."
          << "\n\n"
          << "A common source for this kind of problem is that you "
          << "are passing a 'fully distributed' vector into a function "
          << "that needs read access to vector elements that correspond "
          << "to degrees of freedom on ghost cells (or at least to "
          << "'locally active' degrees of freedom that are not also "
          << "'locally owned'). You need to pass a vector that has these "
          << "elements as ghost entries.");
 
      private:
        Vector &vector;
 
        const size_type index;
 
        // Make the vector class a friend, so that it can create objects of the
        // present type.
        friend class ::LinearAlgebra::EpetraWrappers::Vector;
      }; // class VectorReference
 
    } // namespace internal
    class Vector : public ReadVector<VectorTraits::value_type>,
                   public Subscriptor
    {
    public:
      using value_type      = VectorTraits::value_type;
      using size_type       = VectorTraits::size_type;
      using real_type       = value_type;
      using reference       = internal::VectorReference;
      using const_reference = const internal::VectorReference;
 
      Vector();
 
      Vector(const Vector &V);
 
      explicit Vector(const IndexSet &parallel_partitioner,
                      const MPI_Comm  communicator);
 
      void
      reinit(const IndexSet &parallel_partitioner,
             const MPI_Comm  communicator,
             const bool      omit_zeroing_entries = false);
 
      void
      reinit(const Vector &V, const bool omit_zeroing_entries = false);
 
      virtual void
      extract_subvector_to(
        const ArrayView<const types::global_dof_index> &indices,
        ArrayView<double> &elements) const override;
 
      Vector &
      operator=(const Vector &V);
 
      Vector &
      operator=(const double s);
 
      void
      import_elements(
        const ReadWriteVector<double> &V,
        VectorOperation::values        operation,
        const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
          &communication_pattern = {});
 
      DEAL_II_DEPRECATED
      void
      import(
        const ReadWriteVector<double> &V,
        const VectorOperation::values  operation,
        const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
          &communication_pattern = {})
      {
        import_elements(V, operation, communication_pattern);
      }
 
      reference
      operator()(const size_type index);
 
      value_type
      operator()(const size_type index) const;
 
      reference
      operator[](const size_type index);
 
      value_type
      operator[](const size_type index) const;
 
      Vector &
      operator*=(const double factor);
 
      Vector &
      operator/=(const double factor);
 
      Vector &
      operator+=(const Vector &V);
 
      Vector &
      operator-=(const Vector &V);
 
      double
      operator*(const Vector &V) const;
 
      void
      add(const double a);
 
      void
      add(const double a, const Vector &V);
 
      void
      add(const double a, const Vector &V, const double b, const Vector &W);
 
      void
      sadd(const double s, const double a, const Vector &V);
 
      void
      scale(const Vector &scaling_factors);
 
      void
      equ(const double a, const Vector &V);
 
      bool
      all_zero() const;
 
      double
      mean_value() const;
 
      double
      l1_norm() const;
 
      double
      l2_norm() const;
 
      double
      linfty_norm() const;
 
      double
      add_and_dot(const double a, const Vector &V, const Vector &W);
 
      bool
      has_ghost_elements() const;
 
      virtual size_type
      size() const override;
 
      size_type
      locally_owned_size() const;
 
      MPI_Comm
      get_mpi_communicator() const;
 
      ::IndexSet
      locally_owned_elements() const;
 
      void
      compress(const VectorOperation::values operation);
 
      const Epetra_FEVector &
      trilinos_vector() const;
 
      Epetra_FEVector &
      trilinos_vector();
 
      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;
 
      DeclException0(ExcDifferentParallelPartitioning);
 
      DeclException0(ExcVectorTypeNotCompatible);
 
      DeclException1(ExcTrilinosError,
                     int,
                     << "An error with error number " << arg1
                     << " occurred while calling a Trilinos function");
 
    private:
      void
      create_epetra_comm_pattern(const IndexSet &source_index_set,
                                 const MPI_Comm  mpi_comm);
 
      std::unique_ptr<Epetra_FEVector> vector;
 
      ::IndexSet source_stored_elements;
 
      std::shared_ptr<const CommunicationPattern> epetra_comm_pattern;
 
      // Make the reference class a friend.
      friend class internal::VectorReference;
    };
 
#  ifndef DOXYGEN
 
    // VectorReference
    namespace internal
    {
      inline VectorReference::VectorReference(Vector         &vector,
                                              const size_type index)
        : vector(vector)
        , index(index)
      {}
 
 
 
      inline const VectorReference &
      VectorReference::operator=(const VectorReference &r) const
      {
        // as explained in the class
        // documentation, this is not the copy
        // operator. so simply pass on to the
        // "correct" assignment operator
        *this = static_cast<value_type>(r);
 
        return *this;
      }
 
 
 
      inline VectorReference &
      VectorReference::operator=(const VectorReference &r)
      {
        // as above
        *this = static_cast<value_type>(r);
 
        return *this;
      }
 
 
 
      inline const VectorReference &
      VectorReference::operator=(const value_type &value) const
      {
        (*vector.vector)[0][index] = value;
 
        return *this;
      }
 
 
 
      inline const VectorReference &
      VectorReference::operator+=(const value_type &value) const
      {
        value_type new_value       = static_cast<value_type>(*this) + value;
        (*vector.vector)[0][index] = new_value;
 
        return *this;
      }
 
 
 
      inline const VectorReference &
      VectorReference::operator-=(const value_type &value) const
      {
        value_type new_value       = static_cast<value_type>(*this) - value;
        (*vector.vector)[0][index] = new_value;
 
        return *this;
      }
 
 
 
      inline const VectorReference &
      VectorReference::operator*=(const value_type &value) const
      {
        value_type new_value       = static_cast<value_type>(*this) * value;
        (*vector.vector)[0][index] = new_value;
 
        return *this;
      }
 
 
 
      inline const VectorReference &
      VectorReference::operator/=(const value_type &value) const
      {
        value_type new_value       = static_cast<value_type>(*this) / value;
        (*vector.vector)[0][index] = new_value;
 
        return *this;
      }
    } // namespace internal
 
#  endif /* DOXYGEN */
 
 
    inline internal::VectorReference
    Vector::operator()(const size_type index)
    {
      return internal::VectorReference(*this, index);
    }
 
    inline internal::VectorReference
    Vector::operator[](const size_type index)
    {
      return operator()(index);
    }
 
    inline Vector::value_type
    Vector::operator[](const size_type index) const
    {
      return operator()(index);
    }
 
 
    inline bool
    Vector::has_ghost_elements() const
    {
      return false;
    }
  } // namespace EpetraWrappers
} // namespace LinearAlgebra
 
 
template <>
struct is_serial_vector<LinearAlgebra::EpetraWrappers::Vector> : std::false_type
{};
 
DEAL_II_NAMESPACE_CLOSE
 
#endif
 
#endif