petsc_vector.h

一个用来实现偏微分方程中网格的计算库

   * Creates a local copy of the global vector in
   * \p v_local only on processor \p proc_id.  By
   * default the data is sent to processor 0.  This method
   * is useful for outputting data from one processor.
   */
  void localize_to_one (std::vector<T>& v_local,
			const unsigned int proc_id=0) const;
  
  /**
   * Print the contents of the vector in Matlab
   * format. Optionally prints the
   * matrix to the file named \p name.  If \p name
   * is not specified it is dumped to the screen.
   */
  void print_matlab(const std::string name="NULL") const;

  /**
   * Creates a "subvector" from this vector using the rows indices
   * of the "rows" array.
   */
  virtual void create_subvector(NumericVector<T>& subvector,
				const std::vector<unsigned int>& rows) const;

  /**
   * Swaps the raw PETSc vector context pointers.
   */
  void swap (PetscVector<T> &v);

  /**
   * Returns the raw PETSc vector context pointer.  Note this is generally
   * not required in user-level code. Just don't do anything crazy like
   * calling VecDestroy()!
   */
  Vec vec () { libmesh_assert (_vec != NULL); return _vec; }


  
private:

  /**
   * Actual Petsc vector datatype
   * to hold vector entries
   */
  Vec _vec;

  /**
   * This boolean value should only be set to false
   * for the constructor which takes a PETSc Vec object. 
   */
  bool _destroy_vec_on_exit;
};


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



template <typename T>
inline
PetscVector<T>::PetscVector ()
  : _destroy_vec_on_exit(true)
{}



template <typename T>
inline
PetscVector<T>::PetscVector (const unsigned int n)
  : _destroy_vec_on_exit(true)
{
  this->init(n, n, false);
}



template <typename T>
inline
PetscVector<T>::PetscVector (const unsigned int n,
			     const unsigned int n_local)
  : _destroy_vec_on_exit(true)
{
  this->init(n, n_local, false);
}





template <typename T>
inline
PetscVector<T>::PetscVector (Vec v)
  : _destroy_vec_on_exit(false)
{
  this->_vec = v;
  this->_is_initialized = true;
}




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



template <typename T>
inline
void PetscVector<T>::init (const unsigned int n,
			   const unsigned int n_local,
			   const bool fast)
{
  int ierr=0;
  int petsc_n=static_cast<int>(n);
  int petsc_n_local=static_cast<int>(n_local);


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

  
  // create a sequential vector if on only 1 processor 
  if (n_local == n)
    {
      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
    {
      libmesh_assert (n_local < n);
      
      ierr = VecCreateMPI (libMesh::COMM_WORLD, petsc_n_local, petsc_n,
			   &_vec);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      ierr = VecSetFromOptions (_vec);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }  
  
  this->_is_initialized = true;
  
  
  if (fast == false)
    this->zero ();
}



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



template <typename T>
inline
void PetscVector<T>::close ()
{
  libmesh_assert (this->initialized());
  
  int ierr=0;
  
  ierr = VecAssemblyBegin(_vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
  ierr = VecAssemblyEnd(_vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  this->_is_closed = true;
}



template <typename T>
inline
void PetscVector<T>::clear ()
{
  if ((this->initialized()) && (this->_destroy_vec_on_exit))
    {
      int ierr=0;

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

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



template <typename T>
inline
void PetscVector<T>::zero ()
{
  libmesh_assert (this->initialized());
  
  int ierr=0;

  PetscScalar z=0.;

#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
}



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

  *cloned_vector = *this;

  return cloned_vector;
}



template <typename T>
inline
unsigned int PetscVector<T>::size () const
{
  libmesh_assert (this->initialized());
  
  int ierr=0, petsc_size=0;
  
  if (!this->initialized())
    return 0;
  
  ierr = VecGetSize(_vec, &petsc_size);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  return static_cast<unsigned int>(petsc_size);
}



template <typename T>
inline
unsigned int PetscVector<T>::local_size () const
{
  libmesh_assert (this->initialized());
  
  int ierr=0, petsc_size=0;
  
  ierr = VecGetLocalSize(_vec, &petsc_size);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
  return static_cast<unsigned int>(petsc_size);
}



template <typename T>
inline
unsigned int PetscVector<T>::first_local_index () const
{
  libmesh_assert (this->initialized());
  
  int ierr=0, petsc_first=0, petsc_last=0;
  
  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
  return static_cast<unsigned int>(petsc_first);
}



template <typename T>
inline
unsigned int PetscVector<T>::last_local_index () const
{
  libmesh_assert (this->initialized());
  
  int ierr=0, petsc_first=0, petsc_last=0;
  
  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
  return static_cast<unsigned int>(petsc_last);
}



template <typename T>
inline
T PetscVector<T>::operator() (const unsigned int i) const
{
  libmesh_assert (this->initialized());
  libmesh_assert ( ((i >= this->first_local_index()) &&
	    (i <  this->last_local_index())) );

  int ierr=0;
  PetscScalar *values, value=0.;
  

  ierr = VecGetArray(_vec, &values);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
  value = values[i - this->first_local_index()];
  
  ierr = VecRestoreArray (_vec, &values);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
  return static_cast<T>(value);
}



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

  int index=0, ierr=0;
  PetscReal min=0.;

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

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



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

  int index=0, ierr=0;
  PetscReal max=0.;

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

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



template <typename T>
inline
void PetscVector<T>::swap (PetscVector<T> &v)
{
  std::swap(_vec, v._vec);
  std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
}


#endif // #ifdef HAVE_PETSC
#endif // #ifdef __petsc_vector_h__