petsc_vector.c

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

// $Id: petsc_vector.C 2789 2008-04-13 02:24:40Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// 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



// C++ includes

// Local Includes
#include "petsc_vector.h"
#include "petsc_matrix.h"

#ifdef HAVE_PETSC

#include "parallel.h"
#include "utility.h"
#include "dense_vector.h"
#include "petsc_macro.h"



//-----------------------------------------------------------------------
// PetscVector members

// void PetscVector<T>::init (const NumericVector<T>& v, const bool fast)
// {
//   libmesh_error();
  
//   init (v.local_size(), v.size(), fast);

//   vec = dynamic_cast<const PetscVector<T>&>(v).vec;
// }

template <typename T>
T PetscVector<T>::sum () const
{
  libmesh_assert(this->closed());
  
  int ierr=0;
  PetscScalar value=0.;
  
  ierr = VecSum (_vec, &value);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
  return static_cast<T>(value);
}


template <typename T>
Real PetscVector<T>::l1_norm () const
{
  libmesh_assert(this->closed());
  
  int ierr=0;
  PetscReal value=0.;
  
  ierr = VecNorm (_vec, NORM_1, &value);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
  return static_cast<Real>(value);
}



template <typename T>
Real PetscVector<T>::l2_norm () const
{
  libmesh_assert(this->closed());
  
  int ierr=0;
  PetscReal value=0.;
  
  ierr = VecNorm (_vec, NORM_2, &value);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
  return static_cast<Real>(value);
}




template <typename T>
Real PetscVector<T>::linfty_norm () const
{
  libmesh_assert(this->closed());
  
  int ierr=0;
  PetscReal value=0.;
  
  ierr = VecNorm (_vec, NORM_INFINITY, &value);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
  return static_cast<Real>(value);
}




template <typename T>
NumericVector<T>&
PetscVector<T>::operator += (const NumericVector<T>& v)
{
  libmesh_assert(this->closed());
  
  this->add(1., v);
  
  return *this;
}



template <typename T>
NumericVector<T>&
PetscVector<T>::operator -= (const NumericVector<T>& v)
{
  libmesh_assert(this->closed());
  
  this->add(-1., v);
  
  return *this;
}



template <typename T>
void PetscVector<T>::set (const unsigned int i, const T value)
{
  libmesh_assert(i<size());
  
  int ierr=0;
  int i_val = static_cast<int>(i);
  PetscScalar petsc_value = static_cast<PetscScalar>(value);

  ierr = VecSetValues (_vec, 1, &i_val, &petsc_value, INSERT_VALUES);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  this->_is_closed = false;
}



template <typename T>
void PetscVector<T>::add (const unsigned int i, const T value)
{
  libmesh_assert(i<size());
  
  int ierr=0;
  int i_val = static_cast<int>(i);
  PetscScalar petsc_value = static_cast<PetscScalar>(value);

  ierr = VecSetValues (_vec, 1, &i_val, &petsc_value, ADD_VALUES);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  this->_is_closed = false;
}



template <typename T>
void PetscVector<T>::add_vector (const std::vector<T>& v,
				 const std::vector<unsigned int>& dof_indices)
{
  libmesh_assert (v.size() == dof_indices.size());

  for (unsigned int i=0; i<v.size(); i++)
    this->add (dof_indices[i], v[i]);
}



template <typename T>
void PetscVector<T>::add_vector (const NumericVector<T>& V,
				 const std::vector<unsigned int>& dof_indices)
{
  libmesh_assert (V.size() == dof_indices.size());

  for (unsigned int i=0; i<V.size(); i++)
    this->add (dof_indices[i], V(i));
}



template <typename T>
void PetscVector<T>::add_vector (const NumericVector<T>& V_in,
				 const SparseMatrix<T>& A_in)
{
  const PetscVector<T>* V = dynamic_cast<const PetscVector<T>*>(&V_in);
  const PetscMatrix<T>* A = dynamic_cast<const PetscMatrix<T>*>(&A_in);

  libmesh_assert (V != NULL);
  libmesh_assert (A != NULL);
  
  int ierr=0;

  A->close();

  // The const_cast<> is not elegant, but it is required since PETSc
  // is not const-correct.  
  ierr = MatMultAdd(const_cast<PetscMatrix<T>*>(A)->mat(), V->_vec, _vec, _vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr); 
}



template <typename T>
void PetscVector<T>::add_vector (const DenseVector<T>& V,
				 const std::vector<unsigned int>& dof_indices)
{
  libmesh_assert (V.size() == dof_indices.size());

  for (unsigned int i=0; i<V.size(); i++)
    this->add (dof_indices[i], V(i));
}



template <typename T>
void PetscVector<T>::add (const T v_in)
{
  int ierr=0;
  PetscScalar* values;
  const PetscScalar v = static_cast<PetscScalar>(v_in);  
  const int n   = static_cast<int>(this->local_size());
  const int fli = static_cast<int>(this->first_local_index());
  
  for (int i=0; i<n; i++)
    {
      ierr = VecGetArray (_vec, &values);
  	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      int ig = fli + i;      
      
      PetscScalar value = (values[ig] + v);
      
      ierr = VecRestoreArray (_vec, &values);
  	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      ierr = VecSetValues (_vec, 1, &ig, &value, INSERT_VALUES);
 	     CHKERRABORT(libMesh::COMM_WORLD,ierr); 
    }

  this->_is_closed = false;
}



template <typename T>
void PetscVector<T>::add (const NumericVector<T>& v)
{
  this->add (1., v);
}



template <typename T>
void PetscVector<T>::add (const T a_in, const NumericVector<T>& v_in)
{
  int ierr = 0;
  PetscScalar a = static_cast<PetscScalar>(a_in);

  const PetscVector<T>* v = dynamic_cast<const PetscVector<T>*>(&v_in);

  libmesh_assert (v != NULL);
  libmesh_assert(this->size() == v->size());
  

#if PETSC_VERSION_LESS_THAN(2,3,0)
	 
  // 2.2.x & earlier style
  ierr = VecAXPY(&a, v->_vec, _vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

#else
	 
  // 2.3.x & later style
  ierr = VecAXPY(_vec, a, v->_vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
#endif
}



template <typename T>
void PetscVector<T>::insert (const std::vector<T>& v,
			     const std::vector<unsigned int>& dof_indices)
{
  libmesh_assert (v.size() == dof_indices.size());

  for (unsigned int i=0; i<v.size(); i++)
    this->set (dof_indices[i], v[i]);
}



template <typename T>
void PetscVector<T>::insert (const NumericVector<T>& V,
			     const std::vector<unsigned int>& dof_indices)
{
  libmesh_assert (V.size() == dof_indices.size());

  for (unsigned int i=0; i<V.size(); i++)
    this->set (dof_indices[i], V(i));
}



template <typename T>
void PetscVector<T>::insert (const DenseVector<T>& V,
			     const std::vector<unsigned int>& dof_indices)
{
  libmesh_assert (V.size() == dof_indices.size());

  for (unsigned int i=0; i<V.size(); i++)
    this->set (dof_indices[i], V(i));
}



template <typename T>
void PetscVector<T>::scale (const T factor_in)
{
  int ierr = 0;
  PetscScalar factor = static_cast<PetscScalar>(factor_in);
  
#if PETSC_VERSION_LESS_THAN(2,3,0)

  // 2.2.x & earlier style
  ierr = VecScale(&factor, _vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

#else
  
  // 2.3.x & later style	 
  ierr = VecScale(_vec, factor);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

#endif
}




template <typename T>
T PetscVector<T>::dot (const NumericVector<T>& V) const
{
  // Error flag
  int ierr = 0;
  
  // Return value
  PetscScalar value=0.;
  
  // Make sure the NumericVector passed in is really a PetscVector
  const PetscVector<T>* v = dynamic_cast<const PetscVector<T>*>(&V);
  libmesh_assert (v != NULL);

  // 2.3.x (at least) style.  Untested for previous versions.
  ierr = VecDot(this->_vec, v->_vec, &value);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  return static_cast<T>(value);
}




template <typename T>
NumericVector<T>& 
PetscVector<T>::operator = (const T s_in)
{
  libmesh_assert(this->closed());

  int ierr = 0;
  PetscScalar s = static_cast<PetscScalar>(s_in);

  if (this->size() != 0)
    {
#if PETSC_VERSION_LESS_THAN(2,3,0)
      
  // 2.2.x & earlier style
  ierr = VecSet(&s, _vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	     
#else

  // 2.3.x & later style	 
  ierr = VecSet(_vec, s);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	     
#endif
    }
  
  return *this;
}



template <typename T>
NumericVector<T>&
PetscVector<T>::operator = (const NumericVector<T>& v_in)
{
  const PetscVector<T>* v = dynamic_cast<const PetscVector<T>*>(&v_in);

  libmesh_assert (v != NULL);
  
  *this = *v;
  
  return *this;
}



template <typename T>
PetscVector<T>&
PetscVector<T>::operator = (const PetscVector<T>& v)
{
  if (v.initialized())
    {
      this->init (v.size(), v.local_size());
      this->_is_closed      = v._is_closed;
      this->_is_initialized = v._is_initialized;
  
      if (v.size() != 0)
	{
	  int ierr = 0;

	  ierr = VecCopy (v._vec, this->_vec);
       	         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	}
    }
  
  return *this;
}



template <typename T>
NumericVector<T>&
PetscVector<T>::operator = (const std::vector<T>& v)
{
  const unsigned int nl   = this->local_size();
  const unsigned int ioff = this->first_local_index();
  int ierr=0;
  PetscScalar* values;
      
  /**
   * Case 1:  The vector is the same size of
   * The global vector.  Only add the local components.
   */
  if (this->size() == v.size())
    {
      ierr = VecGetArray (_vec, &values);
 	     CHKERRABORT(libMesh::COMM_WORLD,ierr);

      for (unsigned int i=0; i<nl; i++)
	values[i] =  static_cast<PetscScalar>(v[i+ioff]);
      
      ierr = VecRestoreArray (_vec, &values);
	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }

  /**
   * Case 2: The vector is the same size as our local
   * piece.  Insert directly to the local piece.
   */
  else
    {
      libmesh_assert (this->local_size() == v.size());

      ierr = VecGetArray (_vec, &values);
	     CHKERRABORT(libMesh::COMM_WORLD,ierr);

      for (unsigned int i=0; i<nl; i++)
	values[i] = static_cast<PetscScalar>(v[i]);
      
      ierr = VecRestoreArray (_vec, &values);
	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }

  return *this;
}



template <typename T>
void PetscVector<T>::localize (NumericVector<T>& v_local_in) const
{
  PetscVector<T>* v_local = dynamic_cast<PetscVector<T>*>(&v_local_in);

  libmesh_assert (v_local != NULL);
  libmesh_assert (v_local->local_size() == this->size());

  int ierr = 0;
  const int n = this->size();

  IS is;
  VecScatter scatter;

  // Create idx, idx[i] = i;
  std::vector<int> idx(n); Utility::iota (idx.begin(), idx.end(), 0);

  // Create the index set & scatter object
  ierr = ISCreateGeneral(libMesh::COMM_WORLD, n, &idx[0], &is);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  ierr = VecScatterCreate(_vec,          is,