trilinos_epetra_vector.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
#ifndef LIBMESH_TRILINOS_EPETRA_VECTOR_H
#define LIBMESH_TRILINOS_EPETRA_VECTOR_H
 
 
#include "libmesh/libmesh_common.h"
 
#ifdef LIBMESH_TRILINOS_HAVE_EPETRA
 
// Local includes
#include "libmesh/numeric_vector.h"
#include "libmesh/parallel.h"
 
// Trilinos includes
#include "libmesh/ignore_warnings.h"
#include <Epetra_CombineMode.h>
#include <Epetra_Map.h>
#include <Epetra_MultiVector.h>
#include <Epetra_Vector.h>
#include <Epetra_MpiComm.h>
#include "libmesh/restore_warnings.h"
 
// C++ includes
#include <cstddef>
#include <limits>
#include <memory>
#include <mutex>
#include <vector>
 
// Forward declarations
class Epetra_IntSerialDenseVector;
class Epetra_SerialDenseVector;
 
namespace libMesh
{
 
// forward declarations
template <typename T> class SparseMatrix;
 
template <typename T>
class EpetraVector final : public NumericVector<T>
{
public:
 
  explicit
  EpetraVector (const Parallel::Communicator & comm,
                const ParallelType type = AUTOMATIC);
 
  explicit
  EpetraVector (const Parallel::Communicator & comm,
                const numeric_index_type n,
                const ParallelType type = AUTOMATIC);
 
  EpetraVector (const Parallel::Communicator & comm,
                const numeric_index_type n,
                const numeric_index_type n_local,
                const ParallelType type = AUTOMATIC);
 
  EpetraVector (const Parallel::Communicator & comm,
                const numeric_index_type N,
                const numeric_index_type n_local,
                const std::vector<numeric_index_type> & ghost,
                const ParallelType type = AUTOMATIC);
 
  EpetraVector(Epetra_Vector & v,
               const Parallel::Communicator & comm);
 
  EpetraVector (EpetraVector &&) = delete;
  EpetraVector (const EpetraVector &) = delete;
  EpetraVector & operator= (const EpetraVector &) = delete;
  EpetraVector & operator= (EpetraVector &&) = delete;
  virtual ~EpetraVector ();
 
  virtual void close () override;
 
  virtual void clear () noexcept override;
 
  virtual void zero () override;
 
  virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
 
  virtual std::unique_ptr<NumericVector<T>> clone () const override;
 
  virtual void init (const numeric_index_type N,
                     const numeric_index_type n_local,
                     const bool fast=false,
                     const ParallelType type=AUTOMATIC) override;
 
  virtual void init (const numeric_index_type N,
                     const bool fast=false,
                     const ParallelType type=AUTOMATIC) override;
 
  virtual void init (const numeric_index_type N,
                     const numeric_index_type n_local,
                     const std::vector<numeric_index_type> & ghost,
                     const bool fast = false,
                     const ParallelType = AUTOMATIC) override;
 
  virtual void init (const NumericVector<T> & other,
                     const bool fast = false) override;
 
  virtual NumericVector<T> & operator= (const T s) override;
 
  virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
 
  virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
 
  virtual Real min () const override;
 
  virtual Real max () const override;
 
  virtual T sum () const override;
 
  virtual Real l1_norm () const override;
 
  virtual Real l2_norm () const override;
 
  virtual Real linfty_norm () const override;
 
  virtual numeric_index_type size () const override;
 
  virtual numeric_index_type local_size() const override;
 
  virtual numeric_index_type first_local_index() const override;
 
  virtual numeric_index_type last_local_index() const override;
 
  virtual T operator() (const numeric_index_type i) const override;
 
  virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
 
  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
 
  virtual NumericVector<T> & operator *= (const NumericVector<T> & v) override;
 
  virtual NumericVector<T> & operator /= (const NumericVector<T> & v) override;
 
  virtual void reciprocal() override;
 
  virtual void conjugate() override;
 
  virtual void set (const numeric_index_type i, const T value) override;
 
  virtual void add (const numeric_index_type i, const T value) override;
 
  virtual void add (const T s) override;
 
  virtual void add (const NumericVector<T> & v) override;
 
  virtual void add (const T a, const NumericVector<T> & v) override;
 
  using NumericVector<T>::add_vector;
 
  virtual void add_vector (const T * v,
                           const std::vector<numeric_index_type> & dof_indices) override;
 
  virtual void add_vector (const NumericVector<T> & v,
                           const SparseMatrix<T> & A) override;
 
  virtual void add_vector_transpose (const NumericVector<T> & v,
                                     const SparseMatrix<T> & A) override;
 
  using NumericVector<T>::insert;
 
  virtual void insert (const T * v,
                       const std::vector<numeric_index_type> & dof_indices) override;
 
  virtual void scale (const T factor) override;
 
  virtual void abs() override;
 
  virtual T dot(const NumericVector<T> & v) const override;
 
  virtual void localize (std::vector<T> & v_local) const override;
 
  virtual void localize (NumericVector<T> & v_local) const override;
 
  virtual void localize (NumericVector<T> & v_local,
                         const std::vector<numeric_index_type> & send_list) const override;
 
  virtual void localize (std::vector<T> & v_local,
                         const std::vector<numeric_index_type> & indices) const override;
 
  virtual void localize (const numeric_index_type first_local_idx,
                         const numeric_index_type last_local_idx,
                         const std::vector<numeric_index_type> & send_list) override;
 
  virtual void localize_to_one (std::vector<T> & v_local,
                                const processor_id_type proc_id=0) const override;
 
  virtual void pointwise_mult (const NumericVector<T> & vec1,
                               const NumericVector<T> & vec2) override;
 
  virtual void swap (NumericVector<T> & v) override;
 
  virtual std::size_t max_allowed_id() const override;
 
  Epetra_Vector * vec () { libmesh_assert(_vec); return _vec; }
 
private:
 
  Epetra_Vector * _vec;
 
  std::unique_ptr<Epetra_Map> _map;
 
  bool _destroy_vec_on_exit;
 
  // The following were copied (and slightly modified) from
  // Epetra_FEVector.h in order to allow us to use a standard
  // Epetra_Vector... which is more compatible with other Trilinos
  // packages such as NOX.  All of this code is originally under LGPL
 
  int SumIntoGlobalValues(int numIDs,
                          const int * GIDs,
                          const double * values);
 
  int SumIntoGlobalValues(const Epetra_IntSerialDenseVector & GIDs,
                          const Epetra_SerialDenseVector & values);
 
  int ReplaceGlobalValues(int numIDs,
                          const int * GIDs,
                          const double * values);
 
  int ReplaceGlobalValues(const Epetra_IntSerialDenseVector & GIDs,
                          const Epetra_SerialDenseVector & values);
 
  int SumIntoGlobalValues(int numIDs,
                          const int * GIDs,
                          const int * numValuesPerID,
                          const double * values);
 
  int ReplaceGlobalValues(int numIDs,
                          const int * GIDs,
                          const int * numValuesPerID,
                          const double * values);
 
  int GlobalAssemble(Epetra_CombineMode mode = Add);
 
  void setIgnoreNonLocalEntries(bool flag)
  {
    ignoreNonLocalEntries_ = flag;
  }
 
  void FEoperatorequals(const EpetraVector & source);
 
  int inputValues(int numIDs,
                  const int * GIDs,
                  const double * values,
                  bool accumulate);
 
  int inputValues(int numIDs,
                  const int * GIDs,
                  const int * numValuesPerID,
                  const double * values,
                  bool accumulate);
 
  int inputNonlocalValue(int GID,
                         double value,
                         bool accumulate);
 
  int inputNonlocalValues(int GID,
                          int numValues,
                          const double * values,
                          bool accumulate);
 
  void destroyNonlocalData();
 
  int myFirstID_;
  int myNumIDs_;
  double * myCoefs_;
 
  int * nonlocalIDs_;
  int * nonlocalElementSize_;
  int numNonlocalIDs_;
  int allocatedNonlocalLength_;
  double ** nonlocalCoefs_;
 
  unsigned char last_edit;
 
  bool ignoreNonLocalEntries_;
};
 
 
/*----------------------- Inline functions ----------------------------------*/
 
 
 
template <typename T>
inline
EpetraVector<T>::EpetraVector (const Parallel::Communicator & comm,
                               const ParallelType type) :
  NumericVector<T>(comm, type),
  _destroy_vec_on_exit(true),
  myFirstID_(0),
  myNumIDs_(0),
  myCoefs_(nullptr),
  nonlocalIDs_(nullptr),
  nonlocalElementSize_(nullptr),
  numNonlocalIDs_(0),
  allocatedNonlocalLength_(0),
  nonlocalCoefs_(nullptr),
  last_edit(0),
  ignoreNonLocalEntries_(false)
{
  this->_type = type;
}
 
 
 
template <typename T>
inline
EpetraVector<T>::EpetraVector (const Parallel::Communicator & comm,
                               const numeric_index_type n,
                               const ParallelType type) :
  NumericVector<T>(comm, type),
  _destroy_vec_on_exit(true),
  myFirstID_(0),
  myNumIDs_(0),
  myCoefs_(nullptr),
  nonlocalIDs_(nullptr),
  nonlocalElementSize_(nullptr),
  numNonlocalIDs_(0),
  allocatedNonlocalLength_(0),
  nonlocalCoefs_(nullptr),
  last_edit(0),
  ignoreNonLocalEntries_(false)
 
{
  this->init(n, n, false, type);
}
 
 
 
template <typename T>
inline
EpetraVector<T>::EpetraVector (const Parallel::Communicator & comm,
                               const numeric_index_type n,
                               const numeric_index_type n_local,
                               const ParallelType type) :
  NumericVector<T>(comm, type),
  _destroy_vec_on_exit(true),
  myFirstID_(0),
  myNumIDs_(0),
  myCoefs_(nullptr),
  nonlocalIDs_(nullptr),
  nonlocalElementSize_(nullptr),
  numNonlocalIDs_(0),
  allocatedNonlocalLength_(0),
  nonlocalCoefs_(nullptr),
  last_edit(0),
  ignoreNonLocalEntries_(false)
{
  this->init(n, n_local, false, type);
}
 
 
 
 
template <typename T>
inline
EpetraVector<T>::EpetraVector(Epetra_Vector & v,
                              const Parallel::Communicator & comm) :
  NumericVector<T>(comm, AUTOMATIC),
  _destroy_vec_on_exit(false),
  myFirstID_(0),
  myNumIDs_(0),
  myCoefs_(nullptr),
  nonlocalIDs_(nullptr),
  nonlocalElementSize_(nullptr),
  numNonlocalIDs_(0),
  allocatedNonlocalLength_(0),
  nonlocalCoefs_(nullptr),
  last_edit(0),
  ignoreNonLocalEntries_(false)
{
  _vec = &v;
 
  this->_type = PARALLEL; // FIXME - need to determine this from v!
 
  myFirstID_ = _vec->Map().MinMyGID();
  myNumIDs_ = _vec->Map().NumMyElements();
 
  _map = std::make_unique<Epetra_Map>
    (_vec->GlobalLength(),
     _vec->MyLength(),
     0, // IndexBase = 0 for C/C++, 1 for Fortran.
     Epetra_MpiComm (this->comm().get()));
 
  //Currently we impose the restriction that NumVectors==1, so we won't
  //need the LDA argument when calling ExtractView. Hence the "dummy" arg.
  int dummy;
  _vec->ExtractView(&myCoefs_, &dummy);
 
  this->_is_closed = true;
  this->_is_initialized = true;
}
 
 
 
template <typename T>
inline
EpetraVector<T>::EpetraVector (const Parallel::Communicator & comm,
                               const numeric_index_type n,
                               const numeric_index_type n_local,
                               const std::vector<numeric_index_type> & ghost,
                               const ParallelType type) :
  NumericVector<T>(comm, AUTOMATIC),
  _destroy_vec_on_exit(true),
  myFirstID_(0),
  myNumIDs_(0),
  myCoefs_(nullptr),
  nonlocalIDs_(nullptr),
  nonlocalElementSize_(nullptr),
  numNonlocalIDs_(0),
  allocatedNonlocalLength_(0),
  nonlocalCoefs_(nullptr),
  last_edit(0),
  ignoreNonLocalEntries_(false)
{
  this->init(n, n_local, ghost, false, type);
}
 
 
 
// Default implementation for solver packages for which ghosted
// vectors are not yet implemented.
template <class T>
void EpetraVector<T>::init (const NumericVector<T> & other,
                            const bool fast)
{
  this->init(other.size(),other.local_size(),fast,other.type());
}
 
 
 
template <typename T>
inline
EpetraVector<T>::~EpetraVector ()
{
  this->clear ();
}
 
 
 
template <typename T>
inline
void EpetraVector<T>::init (const numeric_index_type n,
                            const numeric_index_type n_local,
                            const bool fast,
                            const ParallelType type)
{
  // We default to allocating n_local local storage
  numeric_index_type my_n_local = n_local;
 
  if (type == AUTOMATIC)
    {
      if (n == n_local)
        this->_type = SERIAL;
      else
        this->_type = PARALLEL;
    }
  else if (type == GHOSTED)
    {
      // We don't yet support GHOSTED Epetra vectors, so to get the
      // same functionality we need a SERIAL vector with local
      // storage allocated for every entry.
      this->_type = SERIAL;
      my_n_local = n;
    }
  else
    this->_type = type;
 
  libmesh_assert ((this->_type==SERIAL && n==my_n_local) ||
                  this->_type==PARALLEL);
 
  _map = std::make_unique<Epetra_Map>
    (static_cast<int>(n),
     my_n_local,
     0,
     Epetra_MpiComm (this->comm().get()));
 
  _vec = new Epetra_Vector(*_map);
 
  myFirstID_ = _vec->Map().MinMyGID();
  myNumIDs_ = _vec->Map().NumMyElements();
 
  // Currently we impose the restriction that NumVectors==1, so we won't
  // need the LDA argument when calling ExtractView. Hence the "dummy" arg.
  int dummy;
  _vec->ExtractView(&myCoefs_, &dummy);
 
  this->_is_initialized = true;
  this->_is_closed = true;
  this->last_edit = 0;
 
  if (fast == false)
    this->zero ();
}
 
 
template <typename T>
inline
void EpetraVector<T>::init (const numeric_index_type n,
                            const numeric_index_type n_local,
                            const std::vector<numeric_index_type> & /*ghost*/,
                            const bool fast,
                            const ParallelType type)
{
  // TODO: we shouldn't ignore the ghost sparsity pattern
  this->init(n, n_local, fast, type);
}
 
 
 
template <typename T>
inline
void EpetraVector<T>::init (const numeric_index_type n,
                            const bool fast,
                            const ParallelType type)
{
  this->init(n,n,fast,type);
}
 
 
 
template <typename T>
inline
void EpetraVector<T>::close ()
{
  libmesh_assert (this->initialized());
 
  // Are we adding or inserting?
  unsigned char global_last_edit = last_edit;
  this->comm().max(global_last_edit);
  libmesh_assert(!last_edit || last_edit == global_last_edit);
 
  if (global_last_edit == 1)
    this->GlobalAssemble(Insert);
  else if (global_last_edit == 2)
    this->GlobalAssemble(Add);
  else
    libmesh_assert(!global_last_edit);
 
  this->_is_closed = true;
  this->last_edit = 0;
}
 
 
 
template <typename T>
inline
void EpetraVector<T>::clear () noexcept
{
  if (this->initialized())
    {
      // We might just be an interface to a user-provided _vec
      if (this->_destroy_vec_on_exit)
        {
          delete _vec;
          _vec = nullptr;
        }
 
      // But we currently always own our own _map
      _map.reset();
    }
 
  this->_is_closed = this->_is_initialized = false;
}
 
 
 
template <typename T>
inline
void EpetraVector<T>::zero ()
{
  libmesh_assert (this->initialized());
  libmesh_assert (this->closed());
 
  _vec->PutScalar(0.0);
}
 
 
 
template <typename T>
inline
std::unique_ptr<NumericVector<T>> EpetraVector<T>::zero_clone () const
{
  std::unique_ptr<NumericVector<T>> cloned_vector =
    std::make_unique<EpetraVector<T>>(this->comm(), AUTOMATIC);
  cloned_vector->init(*this);
  return cloned_vector;
}
 
 
 
template <typename T>
inline
std::unique_ptr<NumericVector<T>> EpetraVector<T>::clone () const
{
  std::unique_ptr<NumericVector<T>> cloned_vector =
    std::make_unique<EpetraVector<T>>(this->comm(), AUTOMATIC);
  cloned_vector->init(*this, true);
  *cloned_vector = *this;
  return cloned_vector;
}
 
 
 
template <typename T>
inline
numeric_index_type EpetraVector<T>::size () const
{
  libmesh_assert (this->initialized());
 
  return _vec->GlobalLength();
}
 
 
 
template <typename T>
inline
numeric_index_type EpetraVector<T>::local_size () const
{
  libmesh_assert (this->initialized());
 
  return _vec->MyLength();
}
 
template <typename T>
inline
numeric_index_type EpetraVector<T>::first_local_index () const
{
  libmesh_assert (this->initialized());
 
  return _vec->Map().MinMyGID();
}
 
 
 
template <typename T>
inline
numeric_index_type EpetraVector<T>::last_local_index () const
{
  libmesh_assert (this->initialized());
 
  return _vec->Map().MaxMyGID()+1;
}
 
 
template <typename T>
inline
T EpetraVector<T>::operator() (const numeric_index_type i) const
{
  libmesh_assert (this->initialized());
  libmesh_assert ( ((i >= this->first_local_index()) &&
                    (i <  this->last_local_index())) );
 
  return (*_vec)[i-this->first_local_index()];
}
 
 
 
template <typename T>
inline
Real EpetraVector<T>::min () const
{
  libmesh_assert (this->initialized());
 
  T value;
 
  _vec->MinValue(&value);
 
  return value;
}
 
 
 
template <typename T>
inline
Real EpetraVector<T>::max() const
{
  libmesh_assert (this->initialized());
 
  T value;
 
  _vec->MaxValue(&value);
 
  return value;
}
 
 
 
template <typename T>
inline
void EpetraVector<T>::swap (NumericVector<T> & other)
{
  NumericVector<T>::swap(other);
 
  EpetraVector<T> & v = cast_ref<EpetraVector<T> &>(other);
 
  std::swap(_vec, v._vec);
  _map.swap(v._map);
  std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
  std::swap(myFirstID_, v.myFirstID_);
  std::swap(myNumIDs_, v.myNumIDs_);
  std::swap(myCoefs_, v.myCoefs_);
  std::swap(nonlocalIDs_, v.nonlocalIDs_);
  std::swap(nonlocalElementSize_, v.nonlocalElementSize_);
  std::swap(numNonlocalIDs_, v.numNonlocalIDs_);
  std::swap(allocatedNonlocalLength_, v.allocatedNonlocalLength_);
  std::swap(nonlocalCoefs_, v.nonlocalCoefs_);
  std::swap(last_edit, v.last_edit);
  std::swap(ignoreNonLocalEntries_, v.ignoreNonLocalEntries_);
}
 
 
 
template <typename T>
inline
std::size_t EpetraVector<T>::max_allowed_id () const
{
  // Epetra_Vector seems to use hard-coded ints in its various indexing routines.
  return std::numeric_limits<int>::max();
}
 
 
 
// Trilinos only got serious about const in version 10.4
inline
int * numeric_trilinos_cast(const numeric_index_type * p)
{
  return reinterpret_cast<int *>(const_cast<numeric_index_type *>(p));
}
 
} // namespace libMesh
 
 
#endif // #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
#endif // LIBMESH_TRILINOS_EPETRA_VECTOR_H