]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Interface to MatShell matrix-type in PETSc
authorsteigemann <steigemann@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 24 Sep 2012 21:57:25 +0000 (21:57 +0000)
committersteigemann <steigemann@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 24 Sep 2012 21:57:25 +0000 (21:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@26691 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/lac/petsc_matrix_free.h [new file with mode: 0644]
deal.II/source/lac/petsc_matrix_free.cc [new file with mode: 0755]

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 (file)
index 0000000..d6a6678
--- /dev/null
@@ -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 <deal.II/base/config.h>
+
+#ifdef DEAL_II_USE_PETSC
+
+#  include <deal.II/lac/exceptions.h>
+#  include <deal.II/lac/petsc_matrix_base.h>
+#  include <deal.II/lac/petsc_vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+namespace PETScWrappers
+{
+/**
+ * Implementation of a parallel matrix class based on PETSc <tt>MatShell</tt> 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
+ *   <tt>vmult(VectorBase &dst, const VectorBase &src)</tt>
+ * which is pure virtual and must be reimplemented in derived classes.
+ * Besides the usual interface, this class has a matrix-vector multiplication
+ *   <tt>vmult(Vec  &dst, const Vec  &src)</tt>
+ * taking PETSc Vec objects, which will be called by
+ *   <tt>matrix_free_mult(Mat A, Vec src, Vec dst)</tt>
+ * 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 <tt>Vec*</tt> 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 <tt>is_symmetric</tt>.
+                                          */
+      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 <tt>is_symmetric</tt>.
+                                          */
+      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);
+
+                                         /**
+                                          * Constructor for the serial case:
+                                          * Same function as
+                                          * <tt>MatrixFree()</tt>, see above,
+                                          * with <tt>communicator = MPI_COMM_WORLD</tt>.
+                                          */
+      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
+                                          * <tt>MatrixFree()</tt>, see above,
+                                          * with <tt>communicator = MPI_COMM_WORLD</tt>.
+                                          */
+      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);
+
+                                         /**
+                                          * 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<unsigned int> &local_rows_per_process,
+                   const std::vector<unsigned int> &local_columns_per_process,
+                   const unsigned int  this_process);
+
+                                         /**
+                                          * Calls the @p reinit() function
+                                          * above with <tt>communicator = MPI_COMM_WORLD</tt>.
+                                          */
+      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 <tt>communicator = MPI_COMM_WORLD</tt>.
+                                          */
+      void 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);
+
+                                       /**
+                                        * 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 <i>dst = M*src</i> with
+                                          * <i>M</i> 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
+                                          * <i>dst = M<sup>T</sup>*src</i> with
+                                          * <i>M</i> 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
+                                          * <i>M*src</i> on <i>dst</i>
+                                          * with <i>M</i> 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
+                                          * <i>M<sup>T</sup>*src</i> to
+                                          * <i>dst</i> with <i>M</i> 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 <tt>vmult(VectorBase &dst, const VectorBase &src)</tt>
+                                          * 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 <i>dst = A*src</i>
+                                          *
+                                          * Source and destination must
+                                          * not be the same vector.
+                                          *
+                                          * This function calls
+                                          * <tt>vmult(Vec &dst, const Vec &src)</tt>
+                                          * 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 (executable)
index 0000000..d9a6ae2
--- /dev/null
@@ -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 <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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.