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>
  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> //
 
 
 
  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> //
 
  
 
  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> //
 
 
 
  template <class T, class Prop, class Allocator>
  template <class T0,  class Allocator0>
  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