vector/PetscVector.cxx

// Copyright (C) 2011-2012, INRIA
// Author(s): Marc Fragu
//
// This file is part of the linear-algebra library Seldon,
// http://seldon.sourceforge.net/.
//
// Seldon 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.
//
// Seldon 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 Seldon. If not, see http://www.gnu.org/licenses/.


#ifndef SELDON_FILE_VECTOR_PETSCVECTOR_CXX


#include "PetscVector.hxx"


namespace Seldon
{


  /****************
   * CONSTRUCTORS *
   ****************/



  template <class T, class Allocator>
  PETScVector<T, Allocator>::PETScVector():
    Vector_Base<T, Allocator>()
  {
    this->m_ = 0;
    petsc_vector_deallocated_ = false;
  }



  template <class T, class Allocator>
  PETScVector<T, Allocator>::PETScVector(int i, MPI_Comm mpi_communicator)
    :Vector_Base<T, Allocator>(i)
  {
    this->m_ = i;
    petsc_vector_deallocated_ = false;
  }



  template <class T, class Allocator>
  PETScVector<T, Allocator>::
  PETScVector(Vec& petsc_vector)
  {
    petsc_vector_deallocated_ = true;
    Copy(petsc_vector);
  }



  template <class T, class Allocator>
  PETScVector<T, Allocator>::
  PETScVector(const PETScVector<T, Allocator>& V)
  {
    Copy(V);
  }


  /**************
   * DESTRUCTOR *
   **************/


  template <class T, class Allocator>
  PETScVector<T, Allocator>::~PETScVector()
  {
    Clear();
  }



  template <class T, class Allocator>
  Vec& PETScVector<T, Allocator>::GetPetscVector()
  {
    return petsc_vector_;
  }



  template <class T, class Allocator>
  const Vec& PETScVector<T, Allocator>::GetPetscVector() const
  {
    return petsc_vector_;
  }



  template <class T, class Allocator>
  void PETScVector<T, Allocator>::SetCommunicator(MPI_Comm mpi_communicator)
  {
    mpi_communicator_ = mpi_communicator;
  }


  /*********************
   * MEMORY MANAGEMENT *
   *********************/



  template <class T, class Allocator>
  inline void PETScVector<T, Allocator>::Clear()
  {
    if (petsc_vector_deallocated_)
      return;
    int ierr;
    ierr = VecDestroy(&petsc_vector_);
    CHKERRABORT(mpi_communicator_, ierr);
    petsc_vector_deallocated_ = true;
  }



  template <class T, class Allocator>
  inline void PETScVector<T, Allocator>::Resize(int n)
  {
    throw Undefined("PETScVector<T, Allocator>::Resize(int n)");
  }


  template <class T, class Allocator>
  inline void PETScVector<T, Allocator>
  ::SetData(int i, typename PETScVector<T, Allocator>::pointer data)
  {
    throw Undefined("PETScVector<T, Allocator>::SetData(int i, "
                    "typename PETScVector<T, Allocator>::pointer data)");
  }



  template <class T, class Allocator>
  void PETScVector<T, Allocator>::Nullify()
  {
    throw Undefined("PETScVector<T, Allocator>::Nullify()");
  }


  /**********************************
   * ELEMENT ACCESS AND AFFECTATION *
   **********************************/



  template <class T, class Allocator>
  inline typename PETScVector<T, Allocator>::value_type
  PETScVector<T, Allocator>::operator() (int i) const
  {
#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongIndex("Vector<PETSc>::operator()",
                       string("Index should be in [0, ") + to_str(this->m_-1)
                       + "], but is equal to " + to_str(i) + ".");
#endif
    int ierr;
    value_type ret[1];
    int index[1];
    index[0] = i;
    ierr = VecGetValues(petsc_vector_, 1, index, ret);
    CHKERRABORT(mpi_communicator_, ierr);
    return ret[0];
  }



  template <class T, class Allocator>
  inline void PETScVector<T, Allocator>
  ::SetBuffer(int i, T value, InsertMode insert_mode)
  {
    int ierr;
    int ix[1] = {i};
    T data[1] = {value};
    ierr = VecSetValues(petsc_vector_, 1, ix,
                        data, INSERT_VALUES);
    CHKERRABORT(mpi_communicator_, ierr);
  }


  template <class T, class Allocator>
  inline void PETScVector<T, Allocator>
  ::Flush()
  {
    int ierr;
    ierr = VecAssemblyBegin(petsc_vector_);
    CHKERRABORT(mpi_communicator_, ierr);
    ierr = VecAssemblyEnd(petsc_vector_);
    CHKERRABORT(mpi_communicator_, ierr);
  }



  template <class T, class Allocator>
  inline void PETScVector<T, Allocator>
  ::GetProcessorRange(int& i, int& j) const
  {
    int ierr;
    ierr = VecGetOwnershipRange(petsc_vector_, &i, &j);
    CHKERRABORT(mpi_communicator_, ierr);
  }



  template <class T, class Allocator>
  inline void PETScVector<T, Allocator>
  ::Copy(const PETScVector<T, Allocator>& X)
  {
    Copy(X.GetPetscVector());
  }



  template <class T, class Allocator>
  inline void PETScVector<T, Allocator>
  ::Copy(const Vec& petsc_vector)
  {
    Clear();

    int ierr;
    ierr = VecGetSize(petsc_vector, &this->m_);
    CHKERRABORT(mpi_communicator_, ierr);
    PetscObjectGetComm(reinterpret_cast<PetscObject>(petsc_vector),
                       &mpi_communicator_);
    CHKERRABORT(mpi_communicator_, ierr);
    ierr = VecDuplicate(petsc_vector, &petsc_vector_);
    CHKERRABORT(mpi_communicator_, ierr);
    ierr = VecCopy(petsc_vector, petsc_vector_);
    CHKERRABORT(mpi_communicator_, ierr);
    petsc_vector_deallocated_ = false;
  }



  template <class T, class Allocator>
  inline void PETScVector<T, Allocator>::Append(const T& x)
  {
    throw Undefined("PETScVector<T, Allocator>::Append(const T& x)");
  }


  /*******************
   * BASIC FUNCTIONS *
   *******************/



  template <class T, class Allocator>
  int PETScVector<T, Allocator>::GetDataSize() const
  {
    return this->m_;
  }



  template <class T, class Allocator>
  int PETScVector<T, Allocator>::GetLocalM() const
  {
    int size;
    VecGetLocalSize(petsc_vector_, &size);
    return size;
  }


  /************************
   * CONVENIENT FUNCTIONS *
   ************************/



  template <class T, class Allocator>
  void PETScVector<T, Allocator>::Zero()
  {
    int ierr;
    ierr = VecSet(petsc_vector_, 0.);
    CHKERRABORT(mpi_communicator_, ierr);
  }


  template <class T, class Allocator>
  void PETScVector<T, Allocator>::Fill()
  {
    int ierr;
    Vector<int> index(this->m_);
    Vector<double> data(this->m_);
    index.Fill();
    data.Fill();
    ierr = VecSetValues(petsc_vector_, this->m_, index.GetData(),
                        data.GetData(), INSERT_VALUES);
    CHKERRABORT(mpi_communicator_, ierr);
    Flush();
  }



  template <class T, class Allocator>
  template <class T0>
  void PETScVector<T, Allocator>::Fill(const T0& x)
  {
    int ierr;
    ierr = VecSet(petsc_vector_, double(x));
    CHKERRABORT(mpi_communicator_, ierr);
    Flush();
  }



  template <class T, class Allocator>
  void PETScVector<T, Allocator>::FillRand()
  {
    int ierr;
    Vector<int> index(this->m_);
    Vector<double> data(this->m_);
    index.Fill();
    data.FillRand();
    ierr = VecSetValues(petsc_vector_, this->m_, index.GetData(),
                        data.GetData(), INSERT_VALUES);
    CHKERRABORT(mpi_communicator_, ierr);
    Flush();
  }


  /*********
   * NORMS *
   *********/



  template <class T, class Allocator>
  typename PETScVector<T, Allocator>::value_type
  PETScVector<T, Allocator>::GetNormInf() const
  {
    int ierr;
    value_type res;
    ierr = VecNorm(petsc_vector_, NORM_INFINITY, &res);
    CHKERRABORT(mpi_communicator_, ierr);
    return res;
  }



  template <class T, class Allocator>
  int PETScVector<T, Allocator>::GetNormInfIndex() const
  {
    throw Undefined("PETScVector<T, Allocator>::GetNormInfIndex()");
  }


  // SEQUENTIAL PETSC VECTOR //


  /****************
   * CONSTRUCTORS *
   ****************/



  template <class T, class Allocator>
  Vector<T, PETScSeq, Allocator>::Vector():
    PETScVector<T, Allocator>()
  {
    this->mpi_communicator_ = MPI_COMM_WORLD;
    this->m_ = 0;
    int ierr;
    ierr = VecCreateSeq(PETSC_COMM_SELF, 0, &this->petsc_vector_);
    this->petsc_vector_deallocated_ = false;
  }



  template <class T, class Allocator>
  Vector<T, PETScSeq, Allocator>::Vector(int i, MPI_Comm mpi_communicator)
    :PETScVector<T, Allocator>(i)
  {
    int ierr;
    this->mpi_communicator_ = mpi_communicator;
    this->m_ = i;
    ierr = VecCreateSeq(PETSC_COMM_SELF, i, &this->petsc_vector_);
    CHKERRABORT(this->mpi_communicator_, ierr);
    ierr = VecSet(this->petsc_vector_, 0.);
    CHKERRABORT(this->mpi_communicator_, ierr);
    this->Flush();
  }



  template <class T, class Allocator>
  Vector<T, PETScSeq, Allocator>::
  Vector(Vec& petsc_vector): PETScVector<T, Allocator>(petsc_vector)
  {
  }



  template <class T, class Allocator>
  Vector<T, PETScSeq, Allocator>::
  Vector(const Vector<T, PETScSeq, Allocator>& V)
  {
    Copy(V);
  }


  /**************
   * DESTRUCTOR *
   **************/


  template <class T, class Allocator>
  Vector<T, PETScSeq, Allocator>::~Vector()
  {
  }



  template <class T, class Allocator>
  void Vector<T, PETScSeq, Allocator>
  ::Copy(const Vector<T, PETScSeq, Allocator>& X)
  {
    Copy(X.GetPetscVector());
  }



  template <class T, class Allocator>
  void Vector<T, PETScSeq, Allocator>
  ::Copy(const Vec& X)
  {
    PETScVector<T, Allocator>::Copy(X);
  }



  template <class T, class Allocator>
  inline void Vector<T, PETScSeq, Allocator>
  ::Reallocate(int i)
  {
    this->Clear();
    int ierr;
    this->m_ = i;
    ierr = VecCreateSeq(PETSC_COMM_SELF, i, &this->petsc_vector_);
    CHKERRABORT(this->mpi_communicator_, ierr);
    Fill(T(0));
    this->Flush();
    this->petsc_vector_deallocated_ = false;
  }



  template <class T, class Allocator>
  inline Vector<T, PETScSeq, Allocator>& Vector<T, PETScSeq, Allocator>
  ::operator= (const Vector<T, PETScSeq, Allocator>& X)
  {
    this->Copy(X);
    return *this;
  }



  template <class T, class Allocator>
  template <class T0>
  Vector<T, PETScSeq, Allocator>&
  Vector<T, PETScSeq, Allocator>::operator= (const T0& x)
  {
    this->Fill(x);
    return *this;
  }



  template <class T, class Allocator>
  template<class T0>
  inline Vector<T, PETScSeq, Allocator>& Vector<T, PETScSeq, Allocator>
  ::operator*= (const T0& alpha)
  {
    int ierr;
    ierr = VecScale(this->petsc_vector_, alpha);
    CHKERRABORT(this->mpi_communicator_, ierr);
    return *this;
  }


  template <class T, class Allocator>
  void Vector<T, PETScSeq, Allocator>::Print() const
  {
    int ierr;
    ierr = VecView(this->petsc_vector_, PETSC_VIEWER_STDOUT_SELF);
  }


  /**************************
   * OUTPUT/INPUT FUNCTIONS *
   **************************/



  template <class T, class Allocator>
  void Vector<T, PETScSeq, Allocator>
  ::Write(string FileName, bool with_size) const
  {
    throw Undefined("Vector<T, PETScSeq, Allocator>"
                    "::Write(string FileName, bool with_size) const");
  }



  template <class T, class Allocator>
  void Vector<T, PETScSeq, Allocator>
  ::Write(ostream& FileStream, bool with_size) const
  {
    throw Undefined("Vector<T, PETScSeq, Allocator>"
                    "::Write(string FileName, bool with_size) const");
  }



  template <class T, class Allocator>
  void Vector<T, PETScSeq, Allocator>::WriteText(string FileName) const
  {
    ofstream FileStream;
    FileStream.precision(cout.precision());
    FileStream.flags(cout.flags());
    FileStream.open(FileName.c_str());
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("Vector<T, PETScSeq, Allocator>"
                    "::WriteText(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
    this->WriteText(FileStream);
    FileStream.close();
  }



  template <class T, class Allocator>
  void Vector<T, PETScSeq, Allocator>::WriteText(ostream& FileStream) const
  {
    throw Undefined("Vector<T, PETScSeq, Allocator>"
                    "::WriteText(ostream& FileStream) const");
  }



  template <class T, class Allocator>
  void Vector<T, PETScSeq, Allocator>
  ::Read(string FileName, bool with_size)
  {
    ifstream FileStream;
    FileStream.open(FileName.c_str());
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("Vector<T, PETScSeq, Allocator>::Read(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
    this->Read(FileStream, with_size);
    FileStream.close();
  }



  template <class T, class Allocator>
  void Vector<T, PETScSeq, Allocator>
  ::Read(istream& FileStream, bool with_size)
  {
    throw Undefined("Vector<T, PETScSeq, Allocator>"
                    "::Read(istream& FileStream, bool with_size)");
  }



  template <class T, class Allocator>
  void Vector<T, PETScSeq, Allocator>::ReadText(string FileName)
  {
    throw Undefined("Vector<T, PETScSeq, Allocator>"
                    "::ReadText(string FileName)");
  }



  template <class T, class Allocator>
  void Vector<T, PETScSeq, Allocator>::ReadText(istream& FileStream)
  {
    throw Undefined("Vector<T, PETScSeq, Allocator>"
                    "::ReadText(istream& FileStream)");
  }



  template <class T, class Allocator>
  ostream& operator << (ostream& out,
                        const Vector<T, PETScSeq, Allocator>& V)
  {
    out << '\t';
    for (int i = 0; i < V.GetLength() - 1; i++)
      out << V(i) << '\t';
    if (V.GetLength() != 0)
      out << V(V.GetLength() - 1);
    return out;
  }


  // PARALLEL PETSC VECTOR //


  /****************
   * CONSTRUCTORS *
   ****************/



  template <class T, class Allocator>
  Vector<T, PETScPar, Allocator>::Vector():
    PETScVector<T, Allocator>()
  {
    this->mpi_communicator_ = MPI_COMM_WORLD;
    int ierr;
    ierr = VecCreateMPI(this->mpi_communicator_,
                        PETSC_DECIDE, 0, &this->petsc_vector_);
    CHKERRABORT(this->mpi_communicator_, ierr);
  }



  template <class T, class Allocator>
  Vector<T, PETScPar, Allocator>::Vector(int i, MPI_Comm mpi_communicator):
    PETScVector<T, Allocator>(i)
  {
    int ierr;
    this->mpi_communicator_ = mpi_communicator;
    ierr = VecCreateMPI(this->mpi_communicator_, PETSC_DECIDE, i,
                        &this->petsc_vector_);
    CHKERRABORT(this->mpi_communicator_, ierr);
    ierr = VecSet(this->petsc_vector_, 0.);
    CHKERRABORT(this->mpi_communicator_, ierr);
    this->Flush();
  }



  template <class T, class Allocator>
  Vector<T, PETScPar, Allocator>::Vector(int i, int Nlocal,
                                         MPI_Comm mpi_communicator):
    PETScVector<T, Allocator>(i)
  {
    int ierr;
    this->mpi_communicator_ = mpi_communicator;
    ierr = VecCreateMPI(this->mpi_communicator_, Nlocal,
                        i, &this->petsc_vector_);
    CHKERRABORT(this->mpi_communicator_, ierr);
    ierr = VecSet(this->petsc_vector_, 0.);
    CHKERRABORT(this->mpi_communicator_, ierr);
    this->Flush();
  }



  template <class T, class Allocator>
  Vector<T, PETScPar, Allocator>::
  Vector(Vec& petsc_vector): PETScVector<T, Allocator>(petsc_vector)
  {
  }



  template <class T, class Allocator>
  Vector<T, PETScPar, Allocator>::
  Vector(const Vector<T, PETScPar, Allocator>& V)
  {
    Copy(V);
  }


  /**************
   * DESTRUCTOR *
   **************/


  template <class T, class Allocator>
  Vector<T, PETScPar, Allocator>::~Vector()
  {
  }



  template <class T, class Allocator>
  void Vector<T, PETScPar, Allocator>
  ::Copy(const Vector<T, PETScPar, Allocator>& X)
  {
    Copy(X.GetPetscVector());
    this->mpi_communicator_ = X.mpi_communicator_;
  }



  template <class T, class Allocator>
  void Vector<T, PETScPar, Allocator>
  ::Copy(const Vec& X)
  {
    PETScVector<T, Allocator>::Copy(X);
  }



  template <class T, class Allocator>
  inline Vector<T, PETScPar, Allocator>& Vector<T, PETScPar, Allocator>
  ::operator= (const Vector<T, PETScPar, Allocator>& X)
  {
    this->Copy(X);
    return *this;
  }



  template <class T, class Allocator>
  template <class T0>
  Vector<T, PETScPar, Allocator>&
  Vector<T, PETScPar, Allocator>::operator= (const T0& x)
  {
    this->Fill(x);
    return *this;
  }



  template <class T, class Allocator>
  template<class T0>
  inline Vector<T, PETScPar, Allocator>& Vector<T, PETScPar, Allocator>
  ::operator*= (const T0& alpha)
  {
    int ierr;
    ierr = VecScale(this->petsc_vector_, alpha);
    CHKERRABORT(this->mpi_communicator_, ierr);
    return *this;
  }


  template <class T, class Allocator>
  void Vector<T, PETScPar, Allocator>::Print() const
  {
    int ierr;
    ierr = VecView(this->petsc_vector_, PETSC_VIEWER_STDOUT_WORLD);
  }



  template <class T, class Allocator>
  inline void Vector<T, PETScPar, Allocator>
  ::Reallocate(int i, int local_size)
  {
    this->Clear();
    int ierr;
    this->m_ = i;
    ierr = VecCreateMPI(this->mpi_communicator_, local_size, i,
                        &this->petsc_vector_);
    CHKERRABORT(this->mpi_communicator_, ierr);
    Fill(T(0));
    this->Flush();
    this->petsc_vector_deallocated_ = false;
  }


  /**************************
   * OUTPUT/INPUT FUNCTIONS *
   **************************/



  template <class T, class Allocator>
  void Vector<T, PETScPar, Allocator>
  ::Write(string FileName, bool with_size) const
  {
    throw Undefined("Vector<T, PETScPar, Allocator>"
                    "::Write(string FileName, bool with_size) const");
  }



  template <class T, class Allocator>
  void Vector<T, PETScPar, Allocator>
  ::Write(ostream& FileStream, bool with_size) const
  {
    int local_n;
    VecGetLocalSize(this->petsc_vector_, &local_n);

    Vector<int> index(local_n);
    index.Fill();
    int i_start, i_end;
    this->GetProcessorRange(i_start, i_end);
    for (int i = 0; i < local_n; i++)
      index(i) += i_start;
    Vector<T> data(local_n);
    VecGetValues(this->petsc_vector_, local_n, index.GetData(),
                 data.GetData());
    data.Write(FileStream, with_size);
  }



  template <class T, class Allocator>
  void Vector<T, PETScPar, Allocator>::WriteText(string FileName) const
  {
    ofstream FileStream;
    FileStream.precision(cout.precision());
    FileStream.flags(cout.flags());
    FileStream.open(FileName.c_str());
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("Vector<T, PETScPar, Allocator>"
                    "::WriteText(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
    this->WriteText(FileStream);
    FileStream.close();

  }



  template <class T, class Allocator>
  void Vector<T, PETScPar, Allocator>::WriteText(ostream& FileStream) const
  {
    int local_n;
    VecGetLocalSize(this->petsc_vector_, &local_n);

    Vector<int> index(local_n);
    index.Fill();
    int i_start, i_end;
    this->GetProcessorRange(i_start, i_end);
    for (int i = 0; i < local_n; i++)
      index(i) += i_start;
    Vector<T> data(local_n);
    VecGetValues(this->petsc_vector_, local_n, index.GetData(),
                 data.GetData());
    data.WriteText(FileStream);
  }



  template <class T, class Allocator>
  void Vector<T, PETScPar, Allocator>
  ::Read(string FileName, bool with_size)
  {
    ifstream FileStream;
    FileStream.open(FileName.c_str());
#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("Vector<T, PETScPar, Allocator>::Read(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif
    this->Read(FileStream, with_size);
    FileStream.close();
  }



  template <class T, class Allocator>
  void Vector<T, PETScPar, Allocator>
  ::Read(istream& FileStream, bool with_size)
  {
    throw Undefined("PETScVector<T, PETScPar, Allocator>"
                    "::Read(istream& FileStream, bool with_size)");
  }



  template <class T, class Allocator>
  void Vector<T, PETScPar, Allocator>::ReadText(string FileName)
  {
    throw Undefined("PETScVector<T, PETScPar, Allocator>"
                    "::ReadText(string FileName)");
  }



  template <class T, class Allocator>
  void Vector<T, PETScPar, Allocator>::ReadText(istream& FileStream)
  {
    throw Undefined("PETScVector<T, PETScPar, Allocator>"
                    "::ReadText(istream& FileStream)");
  }



  template <class T, class Allocator>
  ostream& operator << (ostream& out,
                        const Vector<T, PETScPar, Allocator>& V)
  {
    out << '\t';
    for (int i = 0; i < V.GetLength() - 1; i++)
      out << V(i) << '\t';
    if (V.GetLength() != 0)
      out << V(V.GetLength() - 1);
    return out;
  }


} // namespace Seldon.


#define SELDON_FILE_PETSCVECTOR_CXX
#endif