trilinos_epetra_vector.h

Go to the documentation of this file.

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

// Local includes
#include "libmesh/numeric_vector.h"
#include "libmesh/parallel.h"

// Trilinos includes
#include <Epetra_CombineMode.h>
#include <Epetra_Map.h>
#include <Epetra_MultiVector.h>
#include <Epetra_Vector.h>
#include <Epetra_MpiComm.h>

// C++ includes
#include <cstddef>
#include <vector>

// Forward declarations
class Epetra_IntSerialDenseVector;
class Epetra_SerialDenseVector;

namespace libMesh
{

// forward declarations
template <typename T> class SparseMatrix;

template <typename T>
class EpetraVector : 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
               LIBMESH_CAN_DEFAULT_TO_COMMWORLD);

  ~EpetraVector ();

  void close ();

  void clear ();

  void zero ();

  virtual UniquePtr<NumericVector<T> > zero_clone () const;

  UniquePtr<NumericVector<T> > clone () const;

  void init (const numeric_index_type N,
             const numeric_index_type n_local,
             const bool         fast=false,
             const ParallelType type=AUTOMATIC);

  void init (const numeric_index_type N,
             const bool         fast=false,
             const ParallelType type=AUTOMATIC);

  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);

  virtual void init (const NumericVector<T>& other,
                     const bool fast = false);

  //   /**
  //    * Change the dimension to that of the
  //    * vector \p V. The same applies as for
  //    * the other \p init function.
  //    *
  //    * The elements of \p V are not copied, i.e.
  //    * this function is the same as calling
  //    * \p init(V.size(),fast).
  //    */
  //   void init (const NumericVector<T>& V,
  //      const bool fast=false);

  NumericVector<T> & operator= (const T s);

  NumericVector<T> & operator= (const NumericVector<T> &V);

  EpetraVector<T> & operator= (const EpetraVector<T> &V);

  NumericVector<T> & operator= (const std::vector<T> &v);

  Real min () const;

  Real max () const;

  T sum () const;

  Real l1_norm () const;

  Real l2_norm () const;

  Real linfty_norm () const;

  numeric_index_type size () const;

  numeric_index_type local_size() const;

  numeric_index_type first_local_index() const;

  numeric_index_type last_local_index() const;

  T operator() (const numeric_index_type i) const;

  NumericVector<T> & operator += (const NumericVector<T> &V);

  NumericVector<T> & operator -= (const NumericVector<T> &V);

  virtual NumericVector<T> & operator /= (NumericVector<T> & v);

  virtual void reciprocal();

  virtual void conjugate();

  void set (const numeric_index_type i, const T value);

  void add (const numeric_index_type i, const T value);

  void add (const T s);

  void add (const NumericVector<T>& V);

  void add (const T a, const NumericVector<T>& v);

  using NumericVector<T>::add_vector;

  void add_vector (const T* v,
                   const std::vector<numeric_index_type>& dof_indices);

  void add_vector (const NumericVector<T> &V,
                   const SparseMatrix<T> &A);

  void add_vector_transpose (const NumericVector<T> &V,
                             const SparseMatrix<T> &A_trans);

  using NumericVector<T>::insert;

  virtual void insert (const T* v,
                       const std::vector<numeric_index_type>& dof_indices);

  void scale (const T factor);

  virtual void abs();


  virtual T dot(const NumericVector<T>& V) const;

  void localize (std::vector<T>& v_local) const;

  void localize (NumericVector<T>& v_local) const;

  void localize (NumericVector<T>& v_local,
                 const std::vector<numeric_index_type>& send_list) const;

  void localize (const numeric_index_type first_local_idx,
                 const numeric_index_type last_local_idx,
                 const std::vector<numeric_index_type>& send_list);

  void localize_to_one (std::vector<T>& v_local,
                        const processor_id_type proc_id=0) const;

  virtual void pointwise_mult (const NumericVector<T>& vec1,
                               const NumericVector<T>& vec2);

  virtual void create_subvector (NumericVector<T>& subvector,
                                 const std::vector<numeric_index_type>& rows) const;

  virtual void swap (NumericVector<T> &v);

  Epetra_Vector * vec () { libmesh_assert(_vec); return _vec; }

private:

  Epetra_Vector * _vec;

  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_(NULL),
    nonlocalIDs_(NULL),
    nonlocalElementSize_(NULL),
    numNonlocalIDs_(0),
    allocatedNonlocalLength_(0),
    nonlocalCoefs_(NULL),
    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_(NULL),
    nonlocalIDs_(NULL),
    nonlocalElementSize_(NULL),
    numNonlocalIDs_(0),
    allocatedNonlocalLength_(0),
    nonlocalCoefs_(NULL),
    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_(NULL),
    nonlocalIDs_(NULL),
    nonlocalElementSize_(NULL),
    numNonlocalIDs_(0),
    allocatedNonlocalLength_(0),
    nonlocalCoefs_(NULL),
    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_(NULL),
    nonlocalIDs_(NULL),
    nonlocalElementSize_(NULL),
    numNonlocalIDs_(0),
    allocatedNonlocalLength_(0),
    nonlocalCoefs_(NULL),
    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 = new Epetra_Map(_vec->GlobalLength(),
                        _vec->MyLength(),
                        0,
                        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_(NULL),
    nonlocalIDs_(NULL),
    nonlocalElementSize_(NULL),
    numNonlocalIDs_(0),
    allocatedNonlocalLength_(0),
    nonlocalCoefs_(NULL),
    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 = new 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 ()
{
  if (this->initialized())
    {
      // We might just be an interface to a user-provided _vec
      if (this->_destroy_vec_on_exit)
        {
          delete _vec;
          _vec = NULL;
        }

      // But we currently always own our own _map
      delete _map;
      _map = NULL;
    }

  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
UniquePtr<NumericVector<T> > EpetraVector<T>::zero_clone () const
{
  UniquePtr<NumericVector<T> > cloned_vector
    (new EpetraVector<T>(this->comm(), AUTOMATIC));

  cloned_vector->init(*this);

  return cloned_vector;
}



template <typename T>
inline
UniquePtr<NumericVector<T> > EpetraVector<T>::clone () const
{
  UniquePtr<NumericVector<T> > cloned_vector
    (new 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);
  std::swap(_map, 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_);
}

} // namespace libMesh


#endif // #ifdef HAVE_EPETRA
#endif // LIBMESH_TRILINOS_EPETRA_VECTOR_H