Home
Publications
Presentations
Support
Download
Examples
Developers
Class Docs
Namespaces
Classes
Files
Mailing Lists
Gallery
Wiki

petsc_vector.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2012 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_PETSC_VECTOR_H
#define LIBMESH_PETSC_VECTOR_H


#include "libmesh/libmesh_config.h"

#ifdef LIBMESH_HAVE_PETSC

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

EXTERN_C_FOR_PETSC_BEGIN
# include <petscvec.h>
EXTERN_C_FOR_PETSC_END

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

namespace libMesh
{



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

template <typename T>
class PetscVector : public NumericVector<T>
{
public:

  explicit
  PetscVector (const ParallelType type = AUTOMATIC);

  explicit
  PetscVector (const numeric_index_type n,
               const ParallelType type = AUTOMATIC);

  PetscVector (const numeric_index_type n,
               const numeric_index_type n_local,
               const ParallelType type = AUTOMATIC);

  PetscVector (const numeric_index_type N,
               const numeric_index_type n_local,
               const std::vector<numeric_index_type>& ghost,
               const ParallelType type = AUTOMATIC);

  PetscVector(Vec v);

  ~PetscVector ();

  void close ();

  void clear ();

  void zero ();

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

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

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

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

  PetscVector<T> & operator= (const PetscVector<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;

  numeric_index_type map_global_to_local_index(const numeric_index_type i) const;

  T operator() (const numeric_index_type i) const;

  virtual void get(const std::vector<numeric_index_type>& index, std::vector<T>& values) const;

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

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

  virtual void reciprocal();

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

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

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


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

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

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

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

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

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

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

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


  void scale (const T factor);

  virtual void abs();

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

  virtual T indefinite_dot(const NumericVector<T>&) 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);

  void print_matlab(const std::string name="NULL") const;

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

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

  Vec vec () { libmesh_assert (_vec); return _vec; }



private:

  Vec _vec;

  mutable bool _array_is_present;

#ifndef NDEBUG

  mutable numeric_index_type _local_size;
#endif

  mutable Vec _local_form;

  mutable PetscScalar* _values;

  void _get_array(void) const;

  void _restore_array(void) const;

  typedef std::map<numeric_index_type,numeric_index_type> GlobalToLocalMap;

  GlobalToLocalMap _global_to_local_map;

  bool _destroy_vec_on_exit;
};


/*----------------------- Inline functions ----------------------------------*/



template <typename T>
inline
PetscVector<T>::PetscVector (const ParallelType ptype)
  : _array_is_present(false),
    _local_form(NULL),
    _values(NULL),
    _global_to_local_map(),
    _destroy_vec_on_exit(true)
{
  this->_type = ptype;
}



template <typename T>
inline
PetscVector<T>::PetscVector (const numeric_index_type n,
                             const ParallelType ptype)
  : _array_is_present(false),
    _local_form(NULL),
    _values(NULL),
    _global_to_local_map(),
    _destroy_vec_on_exit(true)
{
  this->init(n, n, false, ptype);
}



template <typename T>
inline
PetscVector<T>::PetscVector (const numeric_index_type n,
                             const numeric_index_type n_local,
                             const ParallelType ptype)
  : _array_is_present(false),
    _local_form(NULL),
    _values(NULL),
    _global_to_local_map(),
    _destroy_vec_on_exit(true)
{
  this->init(n, n_local, false, ptype);
}



template <typename T>
inline
PetscVector<T>::PetscVector (const numeric_index_type n,
                             const numeric_index_type n_local,
                             const std::vector<numeric_index_type>& ghost,
                             const ParallelType ptype)
  : _array_is_present(false),
    _local_form(NULL),
    _values(NULL),
    _global_to_local_map(),
    _destroy_vec_on_exit(true)
{
  this->init(n, n_local, ghost, false, ptype);
}





template <typename T>
inline
PetscVector<T>::PetscVector (Vec v)
  : _array_is_present(false),
    _local_form(NULL),
    _values(NULL),
    _global_to_local_map(),
    _destroy_vec_on_exit(false)
{
  this->_vec = v;
  this->_is_closed = true;
  this->_is_initialized = true;

  /* We need to ask PETSc about the (local to global) ghost value
     mapping and create the inverse mapping out of it.  */
  PetscErrorCode ierr=0;
  PetscInt petsc_local_size=0;
  ierr = VecGetLocalSize(_vec, &petsc_local_size);
  CHKERRABORT(libMesh::COMM_WORLD,ierr);

  // Get the vector type from PETSc.
  // As of Petsc 3.0.0, the VecType #define lost its const-ness, so we
  // need to have it in the code
#if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_VERSION_RELEASE
  // Pre-3.0 and petsc-dev (as of October 2012) use non-const versions
  VecType ptype;
#else
  const VecType ptype;
#endif
  ierr = VecGetType(_vec, &ptype);
  CHKERRABORT(libMesh::COMM_WORLD,ierr);

  if((std::strcmp(ptype,VECSHARED) == 0) || (std::strcmp(ptype,VECMPI) == 0))
  {
#if PETSC_VERSION_RELEASE && PETSC_VERSION_LESS_THAN(3,1,1)
    ISLocalToGlobalMapping mapping = _vec->mapping;
#else
    ISLocalToGlobalMapping mapping;
    ierr = VecGetLocalToGlobalMapping(_vec, &mapping);
    CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif

    // If is a sparsely stored vector, set up our new mapping
    if (mapping)
    {
      const numeric_index_type my_local_size = static_cast<numeric_index_type>(petsc_local_size);
      const numeric_index_type ghost_begin = static_cast<numeric_index_type>(petsc_local_size);
      const numeric_index_type ghost_end = static_cast<numeric_index_type>(mapping->n);
#if PETSC_VERSION_RELEASE && PETSC_VERSION_LESS_THAN(3,1,1)
      const PetscInt *indices = mapping->indices;
#else
      const PetscInt *indices = mapping->indices;
      ierr = ISLocalToGlobalMappingGetIndices(mapping,&indices);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif
      for(numeric_index_type i=ghost_begin; i<ghost_end; i++)
        _global_to_local_map[indices[i]] = i-my_local_size;
      this->_type = GHOSTED;
    }
    else
      this->_type = PARALLEL;
  }
  else
    this->_type = SERIAL;

  this->close();
}




template <typename T>
inline
PetscVector<T>::~PetscVector ()
{
  this->clear ();
}



template <typename T>
inline
void PetscVector<T>::init (const numeric_index_type n,
                           const numeric_index_type n_local,
                           const bool fast,
                           const ParallelType ptype)
{
  PetscErrorCode ierr=0;
  PetscInt petsc_n=static_cast<PetscInt>(n);
  PetscInt petsc_n_local=static_cast<PetscInt>(n_local);


  // Clear initialized vectors
  if (this->initialized())
    this->clear();

  if (ptype == AUTOMATIC)
    {
      if (n == n_local)
        this->_type = SERIAL;
      else
        this->_type = PARALLEL;
    }
  else
    this->_type = ptype;

  libmesh_assert ((this->_type==SERIAL && n==n_local) ||
                  this->_type==PARALLEL);

  // create a sequential vector if on only 1 processor
  if (this->_type == SERIAL)
    {
      ierr = VecCreateSeq (PETSC_COMM_SELF, petsc_n, &_vec);
             CHKERRABORT(PETSC_COMM_SELF,ierr);

      ierr = VecSetFromOptions (_vec);
             CHKERRABORT(PETSC_COMM_SELF,ierr);
    }
  // otherwise create an MPI-enabled vector
  else if (this->_type == PARALLEL)
    {
#ifdef LIBMESH_HAVE_MPI
      libmesh_assert_less_equal (n_local, n);
      ierr = VecCreateMPI (libMesh::COMM_WORLD, petsc_n_local, petsc_n,
                           &_vec);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
#else
      libmesh_assert_equal_to (n_local, n);
      ierr = VecCreateSeq (PETSC_COMM_SELF, petsc_n, &_vec);
             CHKERRABORT(PETSC_COMM_SELF,ierr);
#endif

      ierr = VecSetFromOptions (_vec);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }
  else
    libmesh_error();

  this->_is_initialized = true;
  this->_is_closed = true;


  if (fast == false)
    this->zero ();
}



template <typename T>
inline
void PetscVector<T>::init (const numeric_index_type n,
                           const bool fast,
                           const ParallelType ptype)
{
  this->init(n,n,fast,ptype);
}



template <typename T>
inline
void PetscVector<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 libmesh_dbg_var(ptype))
{
  PetscErrorCode ierr=0;
  PetscInt petsc_n=static_cast<PetscInt>(n);
  PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
  PetscInt petsc_n_ghost=static_cast<PetscInt>(ghost.size());

  // If the mesh is not disjoint, every processor will either have
  // all the dofs, none of the dofs, or some non-zero dofs at the
  // boundary between processors.
  //
  // However we can't assert this, because someone might want to
  // construct a GHOSTED vector which doesn't include neighbor element
  // dofs.  Boyce tried to do so in user code, and we're going to want
  // to do so in System::project_vector().
  //
  // libmesh_assert(n_local == 0 || n_local == n || !ghost.empty());

  libmesh_assert(sizeof(PetscInt) == sizeof(numeric_index_type));

  PetscInt* petsc_ghost = ghost.empty() ? PETSC_NULL :
    const_cast<PetscInt*>(reinterpret_cast<const PetscInt*>(&ghost[0]));

  // Clear initialized vectors
  if (this->initialized())
    this->clear();

  libmesh_assert(ptype == AUTOMATIC || ptype == GHOSTED);
  this->_type = GHOSTED;

  /* Make the global-to-local ghost cell map.  */
  for(numeric_index_type i=0; i<ghost.size(); i++)
    {
      _global_to_local_map[ghost[i]] = i;
    }

  /* Create vector.  */
  ierr = VecCreateGhost (libMesh::COMM_WORLD, petsc_n_local, petsc_n,
                         petsc_n_ghost, petsc_ghost,
                         &_vec);
  CHKERRABORT(libMesh::COMM_WORLD,ierr);

  ierr = VecSetFromOptions (_vec);
  CHKERRABORT(libMesh::COMM_WORLD,ierr);

  this->_is_initialized = true;
  this->_is_closed = true;

  if (fast == false)
    this->zero ();
}



template <typename T>
inline
void PetscVector<T>::init (const NumericVector<T>& other,
                           const bool fast)
{
  // Clear initialized vectors
  if (this->initialized())
    this->clear();

  const PetscVector<T>& v = libmesh_cast_ref<const PetscVector<T>&>(other);

  // Other vector should restore array.
  if(v.initialized())
    {
      v._restore_array();
    }

  this->_global_to_local_map = v._global_to_local_map;
  this->_is_closed      = v._is_closed;
  this->_is_initialized = v._is_initialized;
  this->_type = v._type;

  if (v.size() != 0)
    {
      PetscErrorCode ierr = 0;

      ierr = VecDuplicate (v._vec, &this->_vec);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }

  if (fast == false)
    this->zero ();
}



template <typename T>
inline
void PetscVector<T>::close ()
{
  this->_restore_array();

  PetscErrorCode ierr=0;

  ierr = VecAssemblyBegin(_vec);
  CHKERRABORT(libMesh::COMM_WORLD,ierr);
  ierr = VecAssemblyEnd(_vec);
  CHKERRABORT(libMesh::COMM_WORLD,ierr);

  if(this->type() == GHOSTED)
    {
      ierr = VecGhostUpdateBegin(_vec,INSERT_VALUES,SCATTER_FORWARD);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = VecGhostUpdateEnd(_vec,INSERT_VALUES,SCATTER_FORWARD);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }

  this->_is_closed = true;
}



template <typename T>
inline
void PetscVector<T>::clear ()
{
  if (this->initialized())
    this->_restore_array();

  if ((this->initialized()) && (this->_destroy_vec_on_exit))
    {
      PetscErrorCode ierr=0;

      ierr = LibMeshVecDestroy(&_vec);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }

  this->_is_closed = this->_is_initialized = false;

  _global_to_local_map.clear();
}



template <typename T>
inline
void PetscVector<T>::zero ()
{
  libmesh_assert(this->closed());

  this->_restore_array();

  PetscErrorCode ierr=0;

  PetscScalar z=0.;

  if(this->type() != GHOSTED)
    {
#if PETSC_VERSION_LESS_THAN(2,3,0)
      // 2.2.x & earlier style
      ierr = VecSet (&z, _vec);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
#else
      // 2.3.x & newer
      ierr = VecSet (_vec, z);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif
    }
  else
    {
      /* Vectors that include ghost values require a special
         handling.  */
      Vec loc_vec;
      ierr = VecGhostGetLocalForm (_vec,&loc_vec);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
#if PETSC_VERSION_LESS_THAN(2,3,0)
      // 2.2.x & earlier style
      ierr = VecSet (&z, loc_vec);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
#else
      // 2.3.x & newer
      ierr = VecSet (loc_vec, z);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif
      ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }
}



template <typename T>
inline
AutoPtr<NumericVector<T> > PetscVector<T>::zero_clone () const
{
  AutoPtr<NumericVector<T> > cloned_vector (new PetscVector<T>);

  cloned_vector->init(*this);

  return cloned_vector;
}



template <typename T>
inline
AutoPtr<NumericVector<T> > PetscVector<T>::clone () const
{
  AutoPtr<NumericVector<T> > cloned_vector (new PetscVector<T>);

  cloned_vector->init(*this, true);

  *cloned_vector = *this;

  return cloned_vector;
}



template <typename T>
inline
numeric_index_type PetscVector<T>::size () const
{
  libmesh_assert (this->initialized());

  PetscErrorCode ierr=0, petsc_size=0;

  if (!this->initialized())
    return 0;

  ierr = VecGetSize(_vec, &petsc_size);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  return static_cast<numeric_index_type>(petsc_size);
}



template <typename T>
inline
numeric_index_type PetscVector<T>::local_size () const
{
  libmesh_assert (this->initialized());

  PetscErrorCode ierr=0, petsc_size=0;

  ierr = VecGetLocalSize(_vec, &petsc_size);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  return static_cast<numeric_index_type>(petsc_size);
}



template <typename T>
inline
numeric_index_type PetscVector<T>::first_local_index () const
{
  libmesh_assert (this->initialized());

  PetscErrorCode ierr=0, petsc_first=0, petsc_last=0;

  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  return static_cast<numeric_index_type>(petsc_first);
}



template <typename T>
inline
numeric_index_type PetscVector<T>::last_local_index () const
{
  libmesh_assert (this->initialized());

  PetscErrorCode ierr=0, petsc_first=0, petsc_last=0;

  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  return static_cast<numeric_index_type>(petsc_last);
}



template <typename T>
inline
numeric_index_type PetscVector<T>::map_global_to_local_index (const numeric_index_type i) const
{
  libmesh_assert (this->initialized());

  PetscErrorCode ierr=0, petsc_first=0, petsc_last=0;
  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
  CHKERRABORT(libMesh::COMM_WORLD,ierr);
  const numeric_index_type first = static_cast<numeric_index_type>(petsc_first);
  const numeric_index_type last = static_cast<numeric_index_type>(petsc_last);

  if((i>=first) && (i<last))
    {
      return i-first;
    }

  GlobalToLocalMap::const_iterator it = _global_to_local_map.find(i);
#ifndef NDEBUG
  const GlobalToLocalMap::const_iterator end = _global_to_local_map.end();
  if (it == end)
    {
      std::ostringstream error_message;
      error_message << "No index " << i << " in ghosted vector.\n"
                    << "Vector contains [" << first << ',' << last << ")\n";
      GlobalToLocalMap::const_iterator b = _global_to_local_map.begin();
      if (b == end)
        {
          error_message << "And empty ghost array.\n";
        }
      else
        {
          error_message << "And ghost array {" << b->first;
          for (b++; b != end; b++)
            error_message << ',' << b->first;
          error_message << "}\n";
        }

      std::cerr << error_message.str();

      libmesh_error();
    }
  libmesh_assert (it != _global_to_local_map.end());
#endif
  return it->second+last-first;
}



template <typename T>
inline
T PetscVector<T>::operator() (const numeric_index_type i) const
{
  this->_get_array();

  const numeric_index_type local_index = this->map_global_to_local_index(i);

#ifndef NDEBUG
  if(this->type() == GHOSTED)
    {
      libmesh_assert_less (local_index, _local_size);
    }
#endif

  return static_cast<T>(_values[local_index]);
}



template <typename T>
inline
void PetscVector<T>::get(const std::vector<numeric_index_type>& index, std::vector<T>& values) const
{
  this->_get_array();

  const std::size_t num = index.size();
  values.resize(num);

  for(std::size_t i=0; i<num; i++)
    {
      const numeric_index_type local_index = this->map_global_to_local_index(index[i]);
#ifndef NDEBUG
      if(this->type() == GHOSTED)
        {
          libmesh_assert_less (local_index, _local_size);
        }
#endif
      values[i] = static_cast<T>(_values[local_index]);
    }
}



template <typename T>
inline
Real PetscVector<T>::min () const
{
  this->_restore_array();

  PetscErrorCode ierr=0;
  PetscInt index=0;
  PetscReal returnval=0.;

  ierr = VecMin (_vec, &index, &returnval);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  // this return value is correct: VecMin returns a PetscReal
  return static_cast<Real>(returnval);
}



template <typename T>
inline
Real PetscVector<T>::max() const
{
  this->_restore_array();

  PetscErrorCode ierr=0;
  PetscInt index=0;
  PetscReal returnval=0.;

  ierr = VecMax (_vec, &index, &returnval);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  // this return value is correct: VecMax returns a PetscReal
  return static_cast<Real>(returnval);
}



template <typename T>
inline
void PetscVector<T>::swap (NumericVector<T> &other)
{
  NumericVector<T>::swap(other);

  PetscVector<T>& v = libmesh_cast_ref<PetscVector<T>&>(other);

  std::swap(_vec, v._vec);
  std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
  std::swap(_global_to_local_map, v._global_to_local_map);
  std::swap(_array_is_present, v._array_is_present);
  std::swap(_local_form, v._local_form);
  std::swap(_values, v._values);
}



template <typename T>
inline
void PetscVector<T>::_get_array(void) const
{
  libmesh_assert (this->initialized());
  if(!_array_is_present)
    {
      PetscErrorCode ierr=0;
      if(this->type() != GHOSTED)
        {
          ierr = VecGetArray(_vec, &_values);
          CHKERRABORT(libMesh::COMM_WORLD,ierr);
        }
      else
        {
          ierr = VecGhostGetLocalForm (_vec,&_local_form);
          CHKERRABORT(libMesh::COMM_WORLD,ierr);
          ierr = VecGetArray(_local_form, &_values);
          CHKERRABORT(libMesh::COMM_WORLD,ierr);
#ifndef NDEBUG
          PetscInt my_local_size = 0;
          ierr = VecGetLocalSize(_local_form, &my_local_size);
          CHKERRABORT(libMesh::COMM_WORLD,ierr);
          _local_size = static_cast<numeric_index_type>(my_local_size);
#endif
        }
      _array_is_present = true;
    }
}



template <typename T>
inline
void PetscVector<T>::_restore_array(void) const
{
  libmesh_assert (this->initialized());
  if(_array_is_present)
    {
      PetscErrorCode ierr=0;
      if(this->type() != GHOSTED)
        {
          ierr = VecRestoreArray (_vec, &_values);
          CHKERRABORT(libMesh::COMM_WORLD,ierr);
          _values = NULL;
        }
      else
        {
          ierr = VecRestoreArray (_local_form, &_values);
          CHKERRABORT(libMesh::COMM_WORLD,ierr);
          _values = NULL;
          ierr = VecGhostRestoreLocalForm (_vec,&_local_form);
          CHKERRABORT(libMesh::COMM_WORLD,ierr);
          _local_form = NULL;
#ifndef NDEBUG
          _local_size = 0;
#endif
        }
      _array_is_present = false;
    }
}


} // namespace libMesh


#endif // #ifdef LIBMESH_HAVE_PETSC
#endif // LIBMESH_PETSC_VECTOR_H

---

Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:48 UTC

Hosted By: