From: steigemann Date: Mon, 24 Sep 2012 21:57:25 +0000 (+0000) Subject: Interface to MatShell matrix-type in PETSc X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c0bd0dae2e5f71e7bdb6bb2ead18b8b48cfb6d67;p=dealii-svn.git Interface to MatShell matrix-type in PETSc git-svn-id: https://svn.dealii.org/trunk@26691 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/lac/petsc_matrix_free.h b/deal.II/include/deal.II/lac/petsc_matrix_free.h new file mode 100644 index 0000000000..d6a6678503 --- /dev/null +++ b/deal.II/include/deal.II/lac/petsc_matrix_free.h @@ -0,0 +1,377 @@ +//--------------------------------------------------------------------------- +// $Id: petsc_matrix_free.h 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. +// +//--------------------------------------------------------------------------- + + +#ifndef __deal2__petsc_matrix_free_h +#define __deal2__petsc_matrix_free_h + + +#include + +#ifdef DEAL_II_USE_PETSC + +# include +# include +# include + +DEAL_II_NAMESPACE_OPEN + + + +namespace PETScWrappers +{ +/** + * Implementation of a parallel matrix class based on PETSc MatShell matrix-type. + * This base class implements only the interface to the PETSc matrix object, + * while all the functionality is contained in the matrix-vector + * multiplication which must be reimplmented in derived classes. + * + * This interface is an addition to the dealii::MatrixFree class to realize + * user-defined matrix-classes together with PETSc solvers and functionalities. + * See also the documentation of dealii::MatrixFree class and step-37 and step-48. + * + * Similar to other @p PETScWrappers::*::*Matrix classes, the MatrxiFree class + * provides the usual matrix-vector multiplication + * vmult(VectorBase &dst, const VectorBase &src) + * which is pure virtual and must be reimplemented in derived classes. + * Besides the usual interface, this class has a matrix-vector multiplication + * vmult(Vec &dst, const Vec &src) + * taking PETSc Vec objects, which will be called by + * matrix_free_mult(Mat A, Vec src, Vec dst) + * registered as matrix-vector multiplication of this PETSc matrix object. + * The default implementation of the vmult function in the base class translates + * the given PETSc Vec* vectors into a deal.II vector, calls + * the usual vmult function with the usual interface and converts + * the result back to PETSc Vec*. This could be made much more efficient + * in derived classes without allocating new memory. + * + * @ingroup PETScWrappers + * @ingroup Matrix1 + * @author Wolfgang Bangerth, Martin Steigemann, 2012 + */ + class MatrixFree : public MatrixBase + { + public: + + /** + * Default constructor. Create an + * empty matrix object. + */ + MatrixFree (); + + /** + * Create a matrix object of + * dimensions @p m times @p n + * with communication happening + * over the provided @p communicator. + * + * For the meaning of the @p local_rows + * and @p local_columns parameters, + * see the PETScWrappers::MPI::SparseMatrix + * class documentation. + * + * As other PETSc matrices, also the + * the matrix-free object needs to + * have a size and to perform matrix + * vector multiplications efficiently + * in parallel also @p local_rows + * and @p local_columns. But in contrast + * to PETSc::SparseMatrix classes a + * PETSc matrix-free object does not need + * any estimation of non_zero entries + * and has no option is_symmetric. + */ + MatrixFree (const MPI_Comm &communicator, + const unsigned int m, + const unsigned int n, + const unsigned int local_rows, + const unsigned int local_columns); + + /** + * Create a matrix object of + * dimensions @p m times @p n + * with communication happening + * over the provided @p communicator. + * + * As other PETSc matrices, also the + * the matrix-free object needs to + * have a size and to perform matrix + * vector multiplications efficiently + * in parallel also @p local_rows + * and @p local_columns. But in contrast + * to PETSc::SparseMatrix classes a + * PETSc matrix-free object does not need + * any estimation of non_zero entries + * and has no option is_symmetric. + */ + MatrixFree (const MPI_Comm &communicator, + const unsigned int m, + const unsigned int n, + const std::vector &local_rows_per_process, + const std::vector &local_columns_per_process, + const unsigned int this_process); + + /** + * Constructor for the serial case: + * Same function as + * MatrixFree(), see above, + * with communicator = MPI_COMM_WORLD. + */ + MatrixFree (const unsigned int m, + const unsigned int n, + const unsigned int local_rows, + const unsigned int local_columns); + + /** + * Constructor for the serial case: + * Same function as + * MatrixFree(), see above, + * with communicator = MPI_COMM_WORLD. + */ + MatrixFree (const unsigned int m, + const unsigned int n, + const std::vector &local_rows_per_process, + const std::vector &local_columns_per_process, + const unsigned int this_process); + + /** + * Throw away the present matrix and + * generate one that has the same + * properties as if it were created by + * the constructor of this class with + * the same argument list as the + * present function. + */ + void reinit (const MPI_Comm &communicator, + const unsigned int m, + const unsigned int n, + const unsigned int local_rows, + const unsigned int local_columns); + + /** + * Throw away the present matrix and + * generate one that has the same + * properties as if it were created by + * the constructor of this class with + * the same argument list as the + * present function. + */ + void reinit (const MPI_Comm &communicator, + const unsigned int m, + const unsigned int n, + const std::vector &local_rows_per_process, + const std::vector &local_columns_per_process, + const unsigned int this_process); + + /** + * Calls the @p reinit() function + * above with communicator = MPI_COMM_WORLD. + */ + void reinit (const unsigned int m, + const unsigned int n, + const unsigned int local_rows, + const unsigned int local_columns); + + /** + * Calls the @p reinit() function + * above with communicator = MPI_COMM_WORLD. + */ + void reinit (const unsigned int m, + const unsigned int n, + const std::vector &local_rows_per_process, + const std::vector &local_columns_per_process, + const unsigned int this_process); + + /** + * Release all memory and return + * to a state just like after + * having called the default + * constructor. + */ + void clear (); + + /** + * Return a reference to the MPI + * communicator object in use with + * this matrix. + */ + const MPI_Comm & get_mpi_communicator () const; + + /** + * Matrix-vector multiplication: + * let dst = M*src with + * M being this matrix. + * + * Source and destination must + * not be the same vector. + * + * Note that if the current object + * represents a parallel distributed + * matrix (of type + * PETScWrappers::MPI::SparseMatrix), + * then both vectors have to be + * distributed vectors as + * well. Conversely, if the matrix is + * not distributed, then neither of the + * vectors may be. + */ + virtual + void vmult (VectorBase &dst, + const VectorBase &src) const = 0; + + /** + * Matrix-vector multiplication: let + * dst = MT*src with + * M being this matrix. This + * function does the same as @p vmult() + * but takes the transposed matrix. + * + * Source and destination must + * not be the same vector. + * + * Note that if the current object + * represents a parallel distributed + * matrix then both vectors have to be + * distributed vectors as + * well. Conversely, if the matrix is + * not distributed, then neither of the + * vectors may be. + */ + virtual + void Tvmult (VectorBase &dst, + const VectorBase &src) const = 0; + + /** + * Adding Matrix-vector + * multiplication. Add + * M*src on dst + * with M being this + * matrix. + * + * Source and destination must + * not be the same vector. + * + * Note that if the current object + * represents a parallel distributed + * matrix then both vectors have to be + * distributed vectors as + * well. Conversely, if the matrix is + * not distributed, then neither of the + * vectors may be. + */ + virtual + void vmult_add (VectorBase &dst, + const VectorBase &src) const = 0; + + /** + * Adding Matrix-vector + * multiplication. Add + * MT*src to + * dst with M being + * this matrix. This function + * does the same as @p vmult_add() + * but takes the transposed + * matrix. + * + * Source and destination must + * not be the same vector. + * + * Note that if the current object + * represents a parallel distributed + * matrix then both vectors have to be + * distributed vectors as + * well. Conversely, if the matrix is + * not distributed, then neither of the + * vectors may be. + */ + virtual + void Tvmult_add (VectorBase &dst, + const VectorBase &src) const = 0; + + /** + * The matrix-vector multiplication + * called by @p matrix_free_mult(). + * This function can be reimplemented + * in derived classes for efficiency. The default + * implementation copies the given vectors + * into PETScWrappers::*::Vector + * and calls vmult(VectorBase &dst, const VectorBase &src) + * which is purely virtual and must be reimplemented + * in derived classes. + */ + virtual + void vmult (Vec &dst, const Vec &src) const; + + private: + + /** + * Copy of the communicator object to + * be used for this parallel matrix-free object. + */ + MPI_Comm communicator; + + /** + * Callback-function registered + * as the matrix-vector multiplication + * of this matrix-free object + * called by PETSc routines. + * This function must be static and + * takes a PETSc matrix @p A, + * and vectors @p src and @p dst, + * where dst = A*src + * + * Source and destination must + * not be the same vector. + * + * This function calls + * vmult(Vec &dst, const Vec &src) + * which should be reimplemented in + * derived classes. + */ + static int matrix_free_mult (Mat A, Vec src, Vec dst); + + /** + * Do the actual work for the + * respective @p reinit() function and + * the matching constructor, + * i.e. create a matrix object. Getting rid + * of the previous matrix is left to + * the caller. + */ + void do_reinit (const unsigned int m, + const unsigned int n, + const unsigned int local_rows, + const unsigned int local_columns); + }; + + + +// -------- template and inline functions ---------- + + inline + const MPI_Comm & + MatrixFree::get_mpi_communicator () const + { + return communicator; + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_USE_PETSC + + +/*---------------------------- petsc_matrix_free.h ---------------------------*/ + +#endif +/*---------------------------- petsc_matrix_free.h ---------------------------*/ diff --git a/deal.II/source/lac/petsc_matrix_free.cc b/deal.II/source/lac/petsc_matrix_free.cc new file mode 100755 index 0000000000..d9a6ae207b --- /dev/null +++ b/deal.II/source/lac/petsc_matrix_free.cc @@ -0,0 +1,289 @@ +//--------------------------------------------------------------------------- +// $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 + +#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 &local_rows_per_process, + const std::vector &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 &local_rows_per_process, + const std::vector &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 &local_rows_per_process, + const std::vector &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 &local_rows_per_process, + const std::vector &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