<commit_before>// $Id$

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

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



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

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

//   vec = libmesh_cast_ref<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)
{
  // Make sure the data passed in are really of Petsc types
  const PetscVector<T>* V = libmesh_cast_ptr<const PetscVector<T>*>(&V_in);
  const PetscMatrix<T>* A = libmesh_cast_ptr<const PetscMatrix<T>*>(&A_in);

  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[i] + 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);

  // Make sure the NumericVector passed in is really a PetscVector
  const PetscVector<T>* v = libmesh_cast_ptr<const PetscVector<T>*>(&v_in);

  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>::insert (const DenseSubVector<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>
void PetscVector<T>::abs()
{
  int ierr = 0;
  ierr = VecAbs(_vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);  
}

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 = libmesh_cast_ptr<const PetscVector<T>*>(&V);

  // 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)
{
  // Make sure the NumericVector passed in is really a PetscVector
  const PetscVector<T>* v = libmesh_cast_ptr<const PetscVector<T>*>(&v_in);

  *this = *v;
  
  return *this;
}



template <typename T>
PetscVector<T>&
PetscVector<T>::operator = (const PetscVector<T>& v)
{
  libmesh_assert (this->initialized());
  libmesh_assert (v.initialized());
  libmesh_assert (this->_type == v._type);
  libmesh_assert (this->size() == v.size());
  libmesh_assert (this->local_size() == v.local_size());
  libmesh_assert (this->_global_to_local_map ==
                  v._global_to_local_map);

  if (v.size() != 0)
    {
      int 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)<commit_msg>Bugfix (or bug workaround?) for operator= on ghosted vectors<commit_after>// $Id$

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

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



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

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

//   vec = libmesh_cast_ref<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)
{
  // Make sure the data passed in are really of Petsc types
  const PetscVector<T>* V = libmesh_cast_ptr<const PetscVector<T>*>(&V_in);
  const PetscMatrix<T>* A = libmesh_cast_ptr<const PetscMatrix<T>*>(&A_in);

  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[i] + 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);

  // Make sure the NumericVector passed in is really a PetscVector
  const PetscVector<T>* v = libmesh_cast_ptr<const PetscVector<T>*>(&v_in);

  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>::insert (const DenseSubVector<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>
void PetscVector<T>::abs()
{
  int ierr = 0;
  ierr = VecAbs(_vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);  
}

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 = libmesh_cast_ptr<const PetscVector<T>*>(&V);

  // 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)
{
  // Make sure the NumericVector passed in is really a PetscVector
  const PetscVector<T>* v = libmesh_cast_ptr<const PetscVector<T>*>(&v_in);

  *this = *v;
  
  return *this;
}



template <typename T>
PetscVector<T>&
PetscVector<T>::operator = (const PetscVector<T>& v)
{
  libmesh_assert (this->initialized());
  libmesh_assert (v.initialized());
  libmesh_assert (this->_type == v._type);
  libmesh_assert (this->size() == v.size());
  libmesh_assert (this->local_size() == v.local_size());
  libmesh_assert (this->_global_to_local_map ==
                  v._global_to_local_map);

  if (v.size() != 0)
    {
      int ierr = VecCopy (v._vec, this->_vec);
      	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }
  
  // Make sure ghost dofs are copied
  if (this->type() == GHOSTED)
    this->close();

  return *this;
}



template <typename T>
NumericVector<T>&
PetscVector<T>::operator = (const std::vector<T>& v)
{'fraction_non_alphanumeric': 0.15812742545561723, 'fraction_numerical': 0.008350935796040675, 'mean_word_length': 3.6671251719394773, 'pattern_counts': {'":': 0, '<': 267, '<?xml version=': 0, '>': 325, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4811843', 'n_tokens_mistral': 7852, 'n_tokens_neox': 7477, 'n_words': 1765}