trilinos_epetra_vector.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2025 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