/home/vivien/public_html/.src_seldon/matrix/PetscMatrix.cxx

// Copyright (C) 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_MATRIX_PETSCMATRIX_CXX

#include "PetscMatrix.hxx"


namespace Seldon
{



  template <class T, class Prop, class Storage, class Allocator>
  inline PetscMatrix<T, Prop, Storage, Allocator>::PetscMatrix():
    Matrix_Base<T, Allocator>()
  {
    mpi_communicator_ = MPI_COMM_WORLD;
    MatCreate(PETSC_COMM_WORLD, &petsc_matrix_);
    petsc_matrix_deallocated_ = false;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline PetscMatrix<T, Prop, Storage, Allocator>::PetscMatrix(int i, int j):
    Matrix_Base<T, Allocator>(i, j)
  {
    int ierr;
    mpi_communicator_ = MPI_COMM_WORLD;
    MatCreate(mpi_communicator_, &petsc_matrix_);
    petsc_matrix_deallocated_ = false;
  }


  template <class T, class Prop, class Storage, class Allocator>
  inline PetscMatrix<T, Prop, Storage, Allocator>
  ::PetscMatrix(Mat& A): Matrix_Base<T, Allocator>()
  {
    Copy(A);
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void PetscMatrix<T, Prop, Storage, Allocator>
  ::SetCommunicator(MPI_Comm mpi_communicator)
  {
    mpi_communicator_ = mpi_communicator;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline MPI_Comm PetscMatrix<T, Prop, Storage, Allocator>
  ::GetCommunicator() const
  {
    return mpi_communicator_;
  }


  template <class T, class Prop, class Storage, class Allocator>
  inline PetscMatrix<T, Prop, Storage, Allocator>
  ::~PetscMatrix()
  {
    Clear();
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void PetscMatrix<T, Prop, Storage, Allocator>
  ::Clear()
  {
    if (petsc_matrix_deallocated_)
      return;
    int ierr;
    ierr = MatDestroy(&petsc_matrix_);
    CHKERRABORT(mpi_communicator_, ierr);
    petsc_matrix_deallocated_ = true;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void PetscMatrix<T, Prop, Storage, Allocator>
  ::Nullify()
  {
    throw Undefined("PetscMatrix<T, Prop, Storage, Allocator>"
                    "::Nullify()");
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline Mat& PetscMatrix<T, Prop, Storage, Allocator>
  ::GetPetscMatrix()
  {
    return petsc_matrix_;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline const Mat& PetscMatrix<T, Prop, Storage, Allocator>
  ::GetPetscMatrix() const
  {
    return petsc_matrix_;
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void PetscMatrix<T, Prop, Storage, Allocator>
  ::Resize(int i, int j)
  {
    throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
                    "::Resize(int i, int j)");
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void PetscMatrix<T, Prop, Storage, Allocator>
  ::SetData(int i, int j, pointer data)
  {
    throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
                    "::SetData(int i, int j, pointer data)");
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline typename PetscMatrix<T, Prop, Storage, Allocator>::const_reference
  PetscMatrix<T, Prop, Storage, Allocator>::Val(int i, int j) const
  {
    throw Undefined("PetscMatrix::Val(int i, int j) const");
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline typename PetscMatrix<T, Prop, Storage, Allocator>::reference
  PetscMatrix<T, Prop, Storage, Allocator>::Val(int i, int j)
  {
    throw Undefined("PetscMatrix::Val(int i, int j)");
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline typename PetscMatrix<T, Prop, Storage, Allocator>::reference
  PetscMatrix<T, Prop, Storage, Allocator>::operator[] (int i)
  {
    throw Undefined("PetscMatrix::operator[] (int i)");
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline typename PetscMatrix<T, Prop, Storage, Allocator>::const_reference
  PetscMatrix<T, Prop, Storage, Allocator>::operator[] (int i) const
  {
    throw Undefined("PetscMatrix::operator[] (int i) const");
  }


  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>::Set(int i, int j, T value)
  {
    throw Undefined("PetscMatrix::Set(int i, int j, T value)");
  }



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>
  ::SetBuffer(int i, int j, T value, InsertMode insert_mode = INSERT_VALUES)
  {
    int ierr;
    ierr = MatSetValues(petsc_matrix_, 1, &i, 1,
                        &j, &value, insert_mode);
    CHKERRABORT(mpi_communicator_, ierr);
  }


  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>::Flush() const
  {
    int ierr;
    ierr = MatAssemblyBegin(petsc_matrix_, MAT_FINAL_ASSEMBLY);
    CHKERRABORT(mpi_communicator_, ierr);
    ierr = MatAssemblyEnd(petsc_matrix_, MAT_FINAL_ASSEMBLY);
    CHKERRABORT(mpi_communicator_, ierr);
  }



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>
  ::GetProcessorRowRange(int& i, int& j) const
  {
    int ierr;
    ierr = MatGetOwnershipRange(petsc_matrix_, &i, &j);
    CHKERRABORT(mpi_communicator_, ierr);
  }



  template <class T, class Prop, class Storage, class Allocator>
  inline void PetscMatrix<T, Prop, Storage, Allocator>
  ::Copy(const Mat& A)
  {
    Clear();
    int ierr;
    ierr = MatDuplicate(A, MAT_COPY_VALUES, &petsc_matrix_);
    CHKERRABORT(mpi_communicator_, ierr);
    MatGetSize(A, &this->m_, &this->n_);
    mpi_communicator_ = MPI_COMM_WORLD;
    petsc_matrix_deallocated_ = false;
  }



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>::Zero()
  {
    MatScale(petsc_matrix_, T(0));
  }


  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>::SetIdentity()
  {
    throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
                    "::SetIdentity()");
  }



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>::Fill()
  {
    throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
                    "::Fill()");
  }



  template <class T, class Prop, class Storage, class Allocator>
  template <class T0>
  void PetscMatrix<T, Prop, Storage, Allocator>::Fill(const T0& x)
  {
    throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
                    "::Fill(const T0& x)");
  }



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>::FillRand()
  {
    throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
                    "::FillRand()");
  }



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>
  ::Print(int a, int b, int m, int n) const
  {
    throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
                    "::Print(int a, int b, int m, int n) const");
  }



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>::Print(int l) const
  {
    throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
                    "::Print(int l) const");
  }


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



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>
  ::Write(string FileName, bool with_size) const
  {
    ofstream FileStream;
    FileStream.open(FileName.c_str());

#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("PetscMatrix::Write(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif

    this->Write(FileStream, with_size);

    FileStream.close();
  }



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>
  ::Write(ostream& FileStream, bool with_size) const
  {
    throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
                    "::Write(ostream& FileStream, bool with_size) const");
  }



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, 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("PetscMatrix::WriteText(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif

    this->WriteText(FileStream);

    FileStream.close();
  }



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



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>::Read(string FileName)
  {
    ifstream FileStream;
    FileStream.open(FileName.c_str());

#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("PetscMatrix::Read(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif

    this->Read(FileStream);

    FileStream.close();
  }



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>
  ::Read(istream& FileStream)
  {
    throw Undefined("void PetscMatrix<T, Prop, Storage, Allocator>"
                    "::Read(istream& FileStream) const");
  }



  template <class T, class Prop, class Storage, class Allocator>
  void PetscMatrix<T, Prop, Storage, Allocator>::ReadText(string FileName)
  {
    ifstream FileStream;
    FileStream.open(FileName.c_str());

#ifdef SELDON_CHECK_IO
    // Checks if the file was opened.
    if (!FileStream.is_open())
      throw IOError("Matrix_Pointers::ReadText(string FileName)",
                    string("Unable to open file \"") + FileName + "\".");
#endif

    this->ReadText(FileStream);

    FileStream.close();
  }



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


  // Matrix<PETScSeqDense> //


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



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScSeqDense, Allocator>::Matrix():
    PetscMatrix<T, Prop, RowMajor, Allocator>()
  {
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScSeqDense, Allocator>::Matrix(int i, int j):
    PetscMatrix<T, Prop, RowMajor, Allocator>(i, j)
  {
    MatSetType(this->petsc_matrix_, MATSEQDENSE);
    MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScSeqDense, Allocator>::Matrix(Mat& A):
  PetscMatrix<T, Prop, RowMajor, Allocator>(A)
  {
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScSeqDense, Allocator>
  ::Matrix(const Matrix<T, Prop, PETScSeqDense, Allocator>& A)
  {
    this->petsc_matrix_deallocated_ = true;
    Copy(A);
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, PETScSeqDense, Allocator>
  ::Reallocate(int i, int j)
  {
    this->Clear();
    int ierr;
    MatCreate(MPI_COMM_SELF, &this->petsc_matrix_);
    MatSetSizes(this->petsc_matrix_, i, j, i, j);
    MatSetType(this->petsc_matrix_, MATSEQDENSE);
    this->petsc_matrix_deallocated_ = false;
    this->Flush();
  }



  template <class T, class Prop, class Allocator>
  inline typename Matrix<T, Prop, PETScSeqDense, Allocator>::value_type
  Matrix<T, Prop, PETScSeqDense, Allocator>::operator() (int i, int j)
  {
#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("PetscMatrix<PETScSeqDense>::operator()",
                     string("Index should be in [0, ")
                     + to_str(this->m_-1)
                     + "], but is equal to " + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("PetscMatrix<PETScSeqDense>::operator()",
                     string("Index should be in [0, ") + to_str(this->n_-1)
                     + "], but is equal to " + to_str(j) + ".");
#endif
    PetscInt idxm[1] = {i};
    PetscInt idxn[1] = {j};
    PetscScalar v[1];
    MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
    return v[0];
  }



  template <class T, class Prop, class Allocator>
  inline typename Matrix<T, Prop, PETScSeqDense, Allocator>::value_type
  Matrix<T, Prop, PETScSeqDense, Allocator>::operator() (int i, int j) const
  {
#ifdef SELDON_CHECK_BOUNDS
    if (i < 0 || i >= this->m_)
      throw WrongRow("PetscMatrix<PETScSeqDense>::operator()",
                     string("Index should be in [0, ")
                     + to_str(this->m_-1)
                     + "], but is equal to " + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("PetscMatrix<PETScSeqDense>::operator()",
                     string("Index should be in [0, ")
                     + to_str(this->n_-1)
                     + "], but is equal to " + to_str(j) + ".");
#endif
    PetscInt idxm[1] = {i};
    PetscInt idxn[1] = {j};
    PetscScalar v[1];
    MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
    return v[0];
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, PETScSeqDense, Allocator>
  ::Copy(const Matrix<T, Prop, PETScSeqDense, Allocator>& A)
  {
    this->mpi_communicator_ = A.mpi_communicator_;
    PetscMatrix<T, Prop, RowMajor, Allocator>::Copy(A.GetPetscMatrix());
    this->petsc_matrix_deallocated_ = false;
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScSeqDense, Allocator>&
  Matrix<T, Prop, PETScSeqDense, Allocator>
  ::operator= (const Matrix<T, Prop, PETScSeqDense, Allocator>& A)
  {
    this->Copy(A);
    return *this;
  }



  template <class T, class Prop, class Allocator>
  void Matrix<T, Prop, PETScSeqDense, Allocator>::Print() const
  {
    int ierr;
    ierr = MatView(this->petsc_matrix_, PETSC_VIEWER_STDOUT_SELF);
    CHKERRABORT(this->mpi_communicator_, ierr);
  }


  // Matrix<PETScMPIDense> //


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



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScMPIDense, Allocator>::Matrix():
    PetscMatrix<T, Prop, RowMajor, Allocator>()
  {
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScMPIDense, Allocator>::Matrix(int i, int j):
    PetscMatrix<T, Prop, RowMajor, Allocator>(i, j)
  {
    MatSetType(this->petsc_matrix_, MATMPIDENSE);
    MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScMPIDense, Allocator>::Matrix(Mat& A):
  PetscMatrix<T, Prop, RowMajor, Allocator>(A)
  {
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScMPIDense, Allocator>
  ::Matrix(const Matrix<T, Prop, PETScMPIDense, Allocator>& A)
  {
    this->petsc_matrix_deallocated_ = true;
    Copy(A);
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, PETScMPIDense, Allocator>
  ::Reallocate(int i, int j, int local_i, int local_j)
  {
    this->Clear();
    int ierr;
    MatCreate(this->mpi_communicator_, &this->petsc_matrix_);
    MatSetType(this->petsc_matrix_, MATMPIDENSE);
    MatSetSizes(this->petsc_matrix_, local_i, local_j, i, j);
    this->m_ = i;
    this->n_ = j;
    this->petsc_matrix_deallocated_ = false;
    this->Flush();
  }



  template <class T, class Prop, class Allocator>
  inline typename Matrix<T, Prop, PETScMPIDense, Allocator>::value_type
  Matrix<T, Prop, PETScMPIDense, Allocator>::operator() (int i, int j)
  {
#ifdef SELDON_CHECK_BOUNDS
    int start, end;
    this->GetProcessorRowRange(start, end);
    if (i < start || i >= end)
      throw WrongRow("PetscMatrix<PETScMPIDense>::operator()",
                     string("Index should be in [")
                     + to_str(start)
                     + ", " + to_str(end - 1)
                     + "], but is equal to "
                     + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("PetscMatrix<PETScMPIDense>::operator()",
                     string("Index should be in [0, ") + to_str(this->n_-1)
                     + "], but is equal to " + to_str(j) + ".");
#endif
    PetscInt idxm[1] = {i};
    PetscInt idxn[1] = {j};
    PetscScalar v[1];
    MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
    return v[0];
  }



  template <class T, class Prop, class Allocator>
  inline typename Matrix<T, Prop, PETScMPIDense, Allocator>::value_type
  Matrix<T, Prop, PETScMPIDense, Allocator>::operator() (int i, int j) const
  {
#ifdef SELDON_CHECK_BOUNDS
    int start, end;
    this->GetProcessorRowRange(start, end);
    if (i < start || i >= end)
      throw WrongRow("PetscMatrix<PETScMPIDense>::operator()",
                     string("Index should be in [")
                     + to_str(start)
                     + ", " + to_str(end - 1) + "], but is equal to "
                     + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("PetscMatrix<PETScMPIDense>::operator()",
                     string("Index should be in [0, ")
                     + to_str(this->n_-1)
                     + "], but is equal to " + to_str(j) + ".");
#endif
    PetscInt idxm[1] = {i};
    PetscInt idxn[1] = {j};
    PetscScalar v[1];
    MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
    return v[0];
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, PETScMPIDense, Allocator>
  ::Copy(const Matrix<T, Prop, PETScMPIDense, Allocator>& A)
  {
    this->mpi_communicator_ = A.mpi_communicator_;
    PetscMatrix<T, Prop, RowMajor, Allocator>::Copy(A.GetPetscMatrix());
    this->petsc_matrix_deallocated_ = false;
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScMPIDense, Allocator>&
  Matrix<T, Prop, PETScMPIDense, Allocator>
  ::operator= (const Matrix<T, Prop, PETScMPIDense, Allocator>& A)
  {
    this->Copy(A);
    return *this;
  }



  template <class T, class Prop, class Allocator>
  void Matrix<T, Prop, PETScMPIDense, Allocator>::Print() const
  {
    int ierr;
    ierr = MatView(this->petsc_matrix_, PETSC_VIEWER_STDOUT_WORLD);
    CHKERRABORT(this->mpi_communicator_, ierr);
  }


  // Matrix<PETScMPIAIJ> //


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



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScMPIAIJ, Allocator>::Matrix():
    PetscMatrix<T, Prop, RowMajor, Allocator>()
  {
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScMPIAIJ, Allocator>::Matrix(int i, int j):
    PetscMatrix<T, Prop, RowMajor, Allocator>(i, j)
  {
    MatSetType(this->petsc_matrix_, MATMPIAIJ);
    MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE, i, j);
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScMPIAIJ, Allocator>::Matrix(Mat& A):
  PetscMatrix<T, Prop, RowMajor, Allocator>(A)
  {
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScMPIAIJ, Allocator>
  ::Matrix(const Matrix<T, Prop, PETScMPIAIJ, Allocator>& A)
  {
    this->petsc_matrix_deallocated_ = true;
    Copy(A);
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, PETScMPIAIJ, Allocator>
  ::Reallocate(int i, int j, int local_i, int local_j)
  {
    this->Clear();
    int ierr;
    MatCreate(this->mpi_communicator_, &this->petsc_matrix_);
    MatSetType(this->petsc_matrix_, MATMPIAIJ);
    MatSetSizes(this->petsc_matrix_, local_i, local_j, i, j);
    this->m_ = i;
    this->n_ = j;
    this->petsc_matrix_deallocated_ = false;
    this->Flush();
  }



  template <class T, class Prop, class Allocator>
  inline typename Matrix<T, Prop, PETScMPIAIJ, Allocator>::value_type
  Matrix<T, Prop, PETScMPIAIJ, Allocator>::operator() (int i, int j)
  {
#ifdef SELDON_CHECK_BOUNDS
    int start, end;
    this->GetProcessorRowRange(start, end);
    if (i < start || i >= end)
      throw WrongRow("PetscMatrix<PETScMPIAIJ>::operator()",
                     string("Index should be in [")
                     + to_str(start)
                     + ", " + to_str(end - 1)
                     + "], but is equal to "
                     + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("PetscMatrix<PETScMPIAIJ>::operator()",
                     string("Index should be in [0, ") + to_str(this->n_-1)
                     + "], but is equal to " + to_str(j) + ".");
#endif
    PetscInt idxm[1] = {i};
    PetscInt idxn[1] = {j};
    PetscScalar v[1];
    MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
    return v[0];
  }



  template <class T, class Prop, class Allocator>
  inline typename Matrix<T, Prop, PETScMPIAIJ, Allocator>::value_type
  Matrix<T, Prop, PETScMPIAIJ, Allocator>::operator() (int i, int j) const
  {
#ifdef SELDON_CHECK_BOUNDS
    int start, end;
    this->GetProcessorRowRange(start, end);
    if (i < start || i >= end)
      throw WrongRow("PetscMatrix<PETScMPIAIJ>::operator()",
                     string("Index should be in [")
                     + to_str(start)
                     + ", " + to_str(end - 1) + "], but is equal to "
                     + to_str(i) + ".");
    if (j < 0 || j >= this->n_)
      throw WrongCol("PetscMatrix<PETScMPIAIJ>::operator()",
                     string("Index should be in [0, ")
                     + to_str(this->n_-1)
                     + "], but is equal to " + to_str(j) + ".");
#endif
    PetscInt idxm[1] = {i};
    PetscInt idxn[1] = {j};
    PetscScalar v[1];
    MatGetValues(this->petsc_matrix_, 1, idxm, 1, idxn, v);
    return v[0];
  }



  template <class T, class Prop, class Allocator>
  inline void Matrix<T, Prop, PETScMPIAIJ, Allocator>
  ::Copy(const Matrix<T, Prop, PETScMPIAIJ, Allocator>& A)
  {
    this->mpi_communicator_ = A.mpi_communicator_;
    PetscMatrix<T, Prop, RowMajor, Allocator>::Copy(A.GetPetscMatrix());
    this->petsc_matrix_deallocated_ = false;
  }



  template <class T, class Prop, class Allocator>
  inline Matrix<T, Prop, PETScMPIAIJ, Allocator>&
  Matrix<T, Prop, PETScMPIAIJ, Allocator>
  ::operator= (const Matrix<T, Prop, PETScMPIAIJ, Allocator>& A)
  {
    this->Copy(A);
    return *this;
  }



  template <class T, class Prop, class Allocator>
  template <class T0,  class Allocator0>
  inline void Matrix<T, Prop, PETScMPIAIJ, Allocator>
  ::Copy(const Matrix<T0, General, RowSparse, Allocator0>& A)
  {
    this->Clear();

    int ierr;
    ierr = MatCreate(this->mpi_communicator_, &this->petsc_matrix_);
    CHKERRABORT(this->mpi_communicator_, ierr);

    int ma = A.GetM();
    int na = A.GetN();
    int nnz = A.GetDataSize();
    double *value = A.GetData();
    int *column = A.GetInd();
    int *ptr = A.GetPtr();

    this->m_ = ma;
    this->n_ = na;

    ierr = MatSetType(this->petsc_matrix_, MATMPIAIJ);
    CHKERRABORT(this->mpi_communicator_, ierr);
    ierr = MatSetSizes(this->petsc_matrix_, PETSC_DECIDE, PETSC_DECIDE,
                       this->m_, this->n_);
    CHKERRABORT(this->mpi_communicator_, ierr);
    int loc_start, loc_end;
    ierr = MatGetOwnershipRange(this->petsc_matrix_,
                                &loc_start, &loc_end);
    CHKERRABORT(this->mpi_communicator_, ierr);
    for (int i = loc_start; i < loc_end; i++)
      for (int j = ptr[i]; j < ptr[i + 1]; j++)
        ierr = MatSetValues(this->petsc_matrix_, 1, &i, 1, &column[j],
                            &value[j], INSERT_VALUES);
    ierr = MatAssemblyBegin(this->petsc_matrix_,MAT_FINAL_ASSEMBLY);
    CHKERRABORT(this->mpi_communicator_, ierr);
    ierr = MatAssemblyEnd(this->petsc_matrix_,MAT_FINAL_ASSEMBLY);
    CHKERRABORT(this->mpi_communicator_, ierr);
  }



  template <class T, class Prop, class Allocator>
  int Matrix<T, Prop, PETScMPIAIJ, Allocator>::GetLocalM() const
  {
    int ierr;
    int m;
    ierr = MatGetLocalSize(this->petsc_matrix_, &m, PETSC_NULL);
    CHKERRABORT(this->mpi_communicator_, ierr);
    return m;
  }



  template <class T, class Prop, class Allocator>
  int Matrix<T, Prop, PETScMPIAIJ, Allocator>::GetLocalN() const
  {
    int ierr;
    int n;
    ierr = MatGetLocalSize(this->petsc_matrix_, PETSC_NULL, &n);
    CHKERRABORT(this->mpi_communicator_, ierr);
    return n;
  }



  template <class T, class Prop, class Allocator>
  void Matrix<T, Prop, PETScMPIAIJ, Allocator>::Print() const
  {
    int ierr;
    ierr = MatView(this->petsc_matrix_, PETSC_VIEWER_STDOUT_WORLD);
    CHKERRABORT(this->mpi_communicator_, ierr);
  }


}


#define SELDON_FILE_MATRIX_PETSCMATRIX_CXX
#endif