<commit_before>// The libMesh Finite Element Library.
// Copyright (C) 2002-2012 <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 "libmesh/petsc_vector.h"
#include "libmesh/petsc_matrix.h"

#ifdef LIBMESH_HAVE_PETSC

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

namespace libMesh
{



//-----------------------------------------------------------------------
// 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
{
  this->_restore_array();
  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
{
  this->_restore_array();
  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
{
  this->_restore_array();
  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
{
  this->_restore_array();
  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)
{
  this->_restore_array();
  libmesh_assert(this->closed());

  this->add(1., v);

  return *this;
}



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

  this->add(-1., v);

  return *this;
}



template <typename T>
void PetscVector<T>::set (const unsigned int i, const T value)
{
  this->_restore_array();
  libmesh_assert_less (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>::reciprocal()
{
  int ierr = 0;

  // VecReciprocal has been in PETSc since at least 2.3.3 days
  ierr = VecReciprocal(_vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
}



template <typename T>
void PetscVector<T>::add (const unsigned int i, const T value)
{
  this->_restore_array();
  libmesh_assert_less (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)
{
  this->_restore_array();
  libmesh_assert_equal_to (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_equal_to (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)
{
  this->_restore_array();
  // 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_equal_to (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_transpose (const NumericVector<T>& V_in,
				           const SparseMatrix<T>& A_in)
{
  this->_restore_array();
  // 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 = MatMultTransposeAdd(const_cast<PetscMatrix<T>*>(A)->mat(), V->_vec, _vec, _vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
}

template <typename T>
void PetscVector<T>::add_vector_conjugate_transpose (const NumericVector<T>& V_in,
				                     const SparseMatrix<T>& A_in)
{
  this->_restore_array();
  // 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();

  // Store a temporary copy since MatMultHermitianTransposeAdd doesn't seem to work
  // TODO: Find out why MatMultHermitianTransposeAdd doesn't work, might be a PETSc bug?
  AutoPtr< NumericVector<Number> > this_clone = this->clone();

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

  // Add the temporary copy to the matvec result
  this->add(1., *this_clone);
}


template <typename T>
void PetscVector<T>::add (const T v_in)
{
  this->_restore_array();

  int ierr=0;
  PetscScalar* values;
  const PetscScalar v = static_cast<PetscScalar>(v_in);

  if(this->type() != GHOSTED)
    {
      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);<commit_msg>fix for ancient PETScs<commit_after>// The libMesh Finite Element Library.
// Copyright (C) 2002-2012 <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 "libmesh/petsc_vector.h"
#include "libmesh/petsc_matrix.h"

#ifdef LIBMESH_HAVE_PETSC

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

namespace libMesh
{



//-----------------------------------------------------------------------
// 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
{
  this->_restore_array();
  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
{
  this->_restore_array();
  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
{
  this->_restore_array();
  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
{
  this->_restore_array();
  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)
{
  this->_restore_array();
  libmesh_assert(this->closed());

  this->add(1., v);

  return *this;
}



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

  this->add(-1., v);

  return *this;
}



template <typename T>
void PetscVector<T>::set (const unsigned int i, const T value)
{
  this->_restore_array();
  libmesh_assert_less (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>::reciprocal()
{
  int ierr = 0;

  // VecReciprocal has been in PETSc since at least 2.3.3 days
  ierr = VecReciprocal(_vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
}



template <typename T>
void PetscVector<T>::add (const unsigned int i, const T value)
{
  this->_restore_array();
  libmesh_assert_less (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)
{
  this->_restore_array();
  libmesh_assert_equal_to (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_equal_to (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)
{
  this->_restore_array();
  // 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_equal_to (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_transpose (const NumericVector<T>& V_in,
				           const SparseMatrix<T>& A_in)
{
  this->_restore_array();
  // 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 = MatMultTransposeAdd(const_cast<PetscMatrix<T>*>(A)->mat(), V->_vec, _vec, _vec);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
}

template <typename T>
void PetscVector<T>::add_vector_conjugate_transpose (const NumericVector<T>& V_in,
				                     const SparseMatrix<T>& A_in)
{
  this->_restore_array();
  // 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();

  // Store a temporary copy since MatMultHermitianTransposeAdd doesn't seem to work
  // TODO: Find out why MatMultHermitianTransposeAdd doesn't work, might be a PETSc bug?
  AutoPtr< NumericVector<Number> > this_clone = this->clone();

  // The const_cast<> is not elegant, but it is required since PETSc
  // is not const-correct.
#if PETSC_VERSION_LESS_THAN(3,1,0)
  
  libMesh::out << "MatMultHermitianTranspose was introduced in PETSc 3.1.0,"
	       << "No one has made it backwards compatible with older "
	       << "versions of PETSc so far." << std::endl;
  libmesh_error();

#else
  
  ierr = MatMultHermitianTranspose(const_cast<PetscMatrix<T>*>(A)->mat(), V->_vec, _vec);
  CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
#endif
  
  // Add the temporary copy to the matvec result
  this->add(1., *this_clone);
}


template <typename T>
void PetscVector<T>::add (const T v_in)
{
  this->_restore_array();

  int ierr=0;
  PetscScalar* values;
  const PetscScalar v = static_cast<PetscScalar>(v_in);

  if(this->type() != GHOSTED)
    {
      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);
{'fraction_non_alphanumeric': 0.1509469816937987, 'fraction_numerical': 0.007347817824792551, 'mean_word_length': 4.048928685641189, 'pattern_counts': {'":': 0, '<': 209, '<?xml version=': 0, '>': 273, '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': '6701174', 'n_tokens_mistral': 5841, 'n_tokens_neox': 5507, 'n_words': 1407}