//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
//    $Id: petsc_matrix_free.cc 26043 2012-09-24 19:25:57Z steigemann $
//    Version: $Name$
//
//    Copyright (C) 2012 by the deal.II authors
//
//    This file is subject to QPL and may not be  distributed
//    without copyright and license information. Please refer
//    to the file deal.II/doc/license.html for the  text  and
//    further information on this license.
//
//---------------------------------------------------------------------------

#include <deal.II/lac/petsc_matrix_free.h>

#ifdef DEAL_II_USE_PETSC

DEAL_II_NAMESPACE_OPEN

namespace PETScWrappers
{
  MatrixFree::MatrixFree ()
            : communicator (PETSC_COMM_SELF)
  {
    const int m=0;
    do_reinit (m, m, m, m);
  }

  MatrixFree::MatrixFree (const MPI_Comm     &communicator,
                          const unsigned int  m,
                          const unsigned int  n,
                          const unsigned int  local_rows,
                          const unsigned int  local_columns)
            : communicator (communicator)
  {
    do_reinit (m, n, local_rows, local_columns);
  }

  MatrixFree::MatrixFree (const MPI_Comm     &communicator,
                          const unsigned int  m,
                          const unsigned int  n,
                          const std::vector<unsigned int> &local_rows_per_process,
                          const std::vector<unsigned int> &local_columns_per_process,
                          const unsigned int  this_process)
            : communicator (communicator)
  {
    Assert (local_rows_per_process.size() == local_columns_per_process.size(),
            ExcDimensionMismatch (local_rows_per_process.size(),
                                  local_columns_per_process.size()));
    Assert (this_process < local_rows_per_process.size(),
            ExcInternalError());

    do_reinit (m, n,
               local_rows_per_process[this_process],
               local_columns_per_process[this_process]);
  }

  MatrixFree::MatrixFree (const unsigned int  m,
                          const unsigned int  n,
                          const unsigned int  local_rows,
                          const unsigned int  local_columns)
            : communicator (MPI_COMM_WORLD)
  {
    do_reinit (m, n, local_rows, local_columns);
  }

  MatrixFree::MatrixFree (const unsigned int  m,
                          const unsigned int  n,
                          const std::vector<unsigned int> &local_rows_per_process,
                          const std::vector<unsigned int> &local_columns_per_process,
                          const unsigned int  this_process)
            : communicator (MPI_COMM_WORLD)
  {
    Assert (local_rows_per_process.size() == local_columns_per_process.size(),
            ExcDimensionMismatch (local_rows_per_process.size(),
                                  local_columns_per_process.size()));
    Assert (this_process < local_rows_per_process.size(),
            ExcInternalError());

    do_reinit (m, n,
               local_rows_per_process[this_process],
               local_columns_per_process[this_process]);
  }

  void MatrixFree::reinit (const MPI_Comm     &communicator,
                           const unsigned int  m,
                           const unsigned int  n,
                           const unsigned int  local_rows,
                           const unsigned int  local_columns)
  {
    this->communicator = communicator;

                                     // destroy the matrix and
                                     // generate a new one
#if DEAL_II_PETSC_VERSION_LT(3,2,0)
    int ierr = MatDestroy (matrix);
#else
    int ierr = MatDestroy (&matrix);
#endif
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    do_reinit (m, n, local_rows, local_columns);
  }

  void MatrixFree::reinit (const MPI_Comm     &communicator,
                           const unsigned int  m,
                           const unsigned int  n,
                           const std::vector<unsigned int> &local_rows_per_process,
                           const std::vector<unsigned int> &local_columns_per_process,
                           const unsigned int  this_process)
  {
    Assert (local_rows_per_process.size() == local_columns_per_process.size(),
            ExcDimensionMismatch (local_rows_per_process.size(),
                                  local_columns_per_process.size()));
    Assert (this_process < local_rows_per_process.size(),
            ExcInternalError());

    this->communicator = communicator;

#if DEAL_II_PETSC_VERSION_LT(3,2,0)
    int ierr = MatDestroy (matrix);
#else
    int ierr = MatDestroy (&matrix);
#endif
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    do_reinit (m, n,
               local_rows_per_process[this_process],
               local_columns_per_process[this_process]);
  }

  void MatrixFree::reinit (const unsigned int  m,
                           const unsigned int  n,
                           const unsigned int  local_rows,
                           const unsigned int  local_columns)
  {
    reinit (MPI_COMM_WORLD, m, n, local_rows, local_columns);
  }

  void MatrixFree::reinit (const unsigned int  m,
                           const unsigned int  n,
                           const std::vector<unsigned int> &local_rows_per_process,
                           const std::vector<unsigned int> &local_columns_per_process,
                           const unsigned int  this_process)
  {
    reinit (MPI_COMM_WORLD, m, n, local_rows_per_process, local_columns_per_process, this_process);
  }

  void MatrixFree::clear ()
  {
#if DEAL_II_PETSC_VERSION_LT(3,2,0)
    int ierr = MatDestroy (matrix);
#else
    int ierr = MatDestroy (&matrix);
#endif
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    const int m=0;
    do_reinit (m, m, m, m);
  }

  void MatrixFree::vmult (Vec  &dst, const Vec  &src) const
  {

//TODO: Translate the given PETSc Vec* vector into a deal.II
// vector so we can call the vmult function with the usual
// interface; then convert back. This could be much more
// efficient, if the PETScWrappers::*::Vector classes
// had a way to simply generate such a vector object from
// a given PETSc Vec* object without allocating new memory
// and without taking ownership of the Vec*

    VectorBase  *x = 0;
    VectorBase  *y = 0;
                                     // because we do not know,
                                     // if dst and src are sequential
                                     // or distributed vectors,
                                     // we ask for the vector-type
                                     // and reinit x and y with
                                     // dealii::PETScWrappers::*::Vector:
    const char  *vec_type;
    int ierr = VecGetType (src, &vec_type);

    PetscInt  local_size;
    ierr = VecGetLocalSize (src, &local_size);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    if (strcmp(vec_type,"mpi") == 0)
      {
        PetscInt  size;
        ierr = VecGetSize (src, &size);
        AssertThrow (ierr == 0, ExcPETScError(ierr));

        x = new PETScWrappers::MPI::Vector (this->get_mpi_communicator (), size, local_size);
        y = new PETScWrappers::MPI::Vector (this->get_mpi_communicator (), size, local_size);
      }
    else if (strcmp(vec_type,"seq") == 0)
      {
        x = new PETScWrappers::Vector (local_size);
        y = new PETScWrappers::Vector (local_size);
      }
    else
      AssertThrow (false, ExcMessage("PETScWrappers::MPI::MatrixFree::do_matrix_vector_action: "
                                     "This only works for Petsc Vec Type = VECMPI | VECSEQ"));

                                     // copy src to x
    x->equ(1., PETScWrappers::VectorBase(src));
                                     // and call vmult(x,y) which must
                                     // be reimplemented in derived classes
    vmult (*y, *x);

    y->compress();
                                     // copy the result back to dst
    ierr = VecCopy (&(*(*y)), dst);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    delete (x);
    delete (y);
  }

  int MatrixFree::matrix_free_mult (Mat  A, Vec  src, Vec  dst)
  {
                                     // create a pointer to this MatrixFree
                                     // object and link the given matrix A
                                     // to the matrix-vector multiplication
                                     // of this MatrixFree object,
    MatrixFree  *this_object;
    int ierr = MatShellGetContext (A, &this_object);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

                                     // call vmult of this object:
    this_object->vmult (dst, src);

    return (0);
  }

  void MatrixFree::do_reinit (const unsigned int  m,
                              const unsigned int  n,
                              const unsigned int  local_rows,
                              const unsigned int  local_columns)
  {
    Assert (local_rows <= m, ExcDimensionMismatch (local_rows, m));
    Assert (local_columns <= n, ExcDimensionMismatch (local_columns, n));

    int ierr;
                                     // create a PETSc MatShell matrix-type
                                     // object of dimension m x n and local size
                                     // local_rows x local_columns
    ierr = MatCreateShell(communicator, local_rows, local_columns, m, n, (void*)this, &matrix);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
                                     // register the MatrixFree::matrix_free_mult function
                                     // as the matrix multiplication used by this matrix
    ierr = MatShellSetOperation (matrix, MATOP_MULT,
               (void(*)(void))&dealii::PETScWrappers::MatrixFree::matrix_free_mult);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = MatSetFromOptions (matrix);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }
}

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_USE_PETSC

git-svn-id: https://svn.dealii.org/trunk@26688 0785d39b-7218-0410-832d-ea1e28bc413d