PETSc ReinitHelper for BlockVector

diff --git a/include/deal.II/lac/block_vector.h b/include/deal.II/lac/block_vector.h
index ac8a467d424dec4c8dd969078c01e77b07880e7d..78342c74afb41f9b283a79f649e2f6206c4a2441 100644 (file)
--- a/include/deal.II/lac/block_vector.h
+++ b/include/deal.II/lac/block_vector.h
@@ -488,7 +488,7 @@ namespace internal
     template <typename> class ReinitHelper;
 
     /**
-     * A helper class internally used in linear_operator.h. Specialization for
+     * A helper class used internally in linear_operator.h. Specialization for
      * BlockVector<number>.
      */
     template<typename number>

diff --git a/include/deal.II/lac/petsc_block_sparse_matrix.h b/include/deal.II/lac/petsc_block_sparse_matrix.h
index c4d7233448d212e37a17d04be8d44a38c65cbf64..41115858ca3fffdf612a163ac70bf833d262cd18 100644 (file)
--- a/include/deal.II/lac/petsc_block_sparse_matrix.h
+++ b/include/deal.II/lac/petsc_block_sparse_matrix.h
@@ -206,6 +206,19 @@ namespace PETScWrappers
     void Tvmult (Vector       &dst,
                  const Vector &src) const;
 
+    /**
+     * Return the partitioning of the domain space of this matrix, i.e., the
+     * partitioning of the vectors this matrix has to be multiplied with.
+     */
+    std::vector< size_type > locally_domain_sizes() const;
+
+    /**
+     * Return the partitioning of the range space of this matrix, i.e., the
+     * partitioning of the vectors that are result from matrix-vector
+     * products.
+     */
+    std::vector< size_type > locally_range_sizes() const;
+
     /**
      * Make the clear() function in the base class visible, though it is
      * protected.

diff --git a/include/deal.II/lac/petsc_block_vector.h b/include/deal.II/lac/petsc_block_vector.h
index 9dc85b424ce6d520b950766e18dc011cdd3462cf..e912a3367c4cc4bb8862229cbc8da9fd41a6a625 100644 (file)
--- a/include/deal.II/lac/petsc_block_vector.h
+++ b/include/deal.II/lac/petsc_block_vector.h
@@ -459,6 +459,42 @@ namespace PETScWrappers
 }
 
 
+namespace internal
+{
+  namespace LinearOperator
+  {
+    template <typename> class ReinitHelper;
+
+    /**
+     * A helper class used internally in linear_operator.h. Specialization for
+     * PETScWrappers::BlockVector.
+     */
+    template<>
+    class ReinitHelper<PETScWrappers::BlockVector>
+    {
+    public:
+      template <typename Matrix>
+      static
+      void reinit_range_vector (const Matrix &matrix,
+                                PETScWrappers::BlockVector &v,
+                                bool omit_zeroing_entries)
+      {
+        v.reinit(matrix.locally_range_sizes(), omit_zeroing_entries);
+      }
+
+      template <typename Matrix>
+      static
+      void reinit_domain_vector(const Matrix &matrix,
+                                PETScWrappers::BlockVector &v,
+                                bool omit_zeroing_entries)
+      {
+        v.reinit(matrix.locally_domain_sizes(), omit_zeroing_entries);
+      }
+    };
+
+  } /* namespace LinearOperator */
+} /* namespace internal */
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif  // DEAL_II_WITH_PETSC

diff --git a/include/deal.II/lac/petsc_parallel_block_sparse_matrix.h b/include/deal.II/lac/petsc_parallel_block_sparse_matrix.h
index 1f79fa60edc41316f37ea061c380aff00634de2a..fb5fb17786bf8c56004d536cb8edbf174f0de5ce 100644 (file)
--- a/include/deal.II/lac/petsc_parallel_block_sparse_matrix.h
+++ b/include/deal.II/lac/petsc_parallel_block_sparse_matrix.h
@@ -234,6 +234,19 @@ namespace PETScWrappers
        */
       void collect_sizes ();
 
+      /**
+       * Return the partitioning of the domain space of this matrix, i.e., the
+       * partitioning of the vectors this matrix has to be multiplied with.
+       */
+      std::vector< IndexSet > locally_owned_domain_indices() const;
+
+      /**
+       * Return the partitioning of the range space of this matrix, i.e., the
+       * partitioning of the vectors that are result from matrix-vector
+       * products.
+       */
+      std::vector< IndexSet > locally_owned_range_indices() const;
+
       /**
        * Return a reference to the MPI communicator object in use with this
        * matrix.

diff --git a/include/deal.II/lac/petsc_parallel_block_vector.h b/include/deal.II/lac/petsc_parallel_block_vector.h
index ac8fa4e1b4d4e36a19c82f59b8f7bfd224cb9b4b..9263128aab50e86e60fd0f410708b317734bcf2a 100644 (file)
--- a/include/deal.II/lac/petsc_parallel_block_vector.h
+++ b/include/deal.II/lac/petsc_parallel_block_vector.h
@@ -538,6 +538,42 @@ namespace PETScWrappers
 
 }
 
+namespace internal
+{
+  namespace LinearOperator
+  {
+    template <typename> class ReinitHelper;
+
+    /**
+     * A helper class used internally in linear_operator.h. Specialization for
+     * PETScWrappers::MPI::BlockVector.
+     */
+    template<>
+    class ReinitHelper<PETScWrappers::MPI::BlockVector>
+    {
+    public:
+      template <typename Matrix>
+      static
+      void reinit_range_vector (const Matrix &matrix,
+                                PETScWrappers::MPI::BlockVector &v,
+                                bool omit_zeroing_entries)
+      {
+        v.reinit(matrix.locally_owned_range_indices(), matrix.get_mpi_communicator());
+      }
+
+      template <typename Matrix>
+      static
+      void reinit_domain_vector(const Matrix &matrix,
+                                PETScWrappers::MPI::BlockVector &v,
+                                bool omit_zeroing_entries)
+      {
+        v.reinit(matrix.locally_owned_domain_indices(), matrix.get_mpi_communicator());
+      }
+    };
+
+  } /* namespace LinearOperator */
+} /* namespace internal */
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif  // DEAL_II_WITH_PETSC

diff --git a/include/deal.II/lac/petsc_parallel_sparse_matrix.h b/include/deal.II/lac/petsc_parallel_sparse_matrix.h
index 6d17876da623328b7ee22e005ffe237d9e085995..f8ced1a8b7cf22d5e749ffc07d281a5e018bab88 100644 (file)
--- a/include/deal.II/lac/petsc_parallel_sparse_matrix.h
+++ b/include/deal.II/lac/petsc_parallel_sparse_matrix.h
@@ -387,7 +387,7 @@ namespace PETScWrappers
 
       /**
        * Return the partitioning of the range space of this matrix, i.e., the
-       * partitioning of the vectors that are result from matrix-vector
+       * partitioning of the vectors that result from matrix-vector
        * products.
        */
       IndexSet locally_owned_range_indices() const;

diff --git a/include/deal.II/lac/petsc_parallel_vector.h b/include/deal.II/lac/petsc_parallel_vector.h
index 522a45f937427b566d45ca6570d616baa44fa586..3f005fea9c7d7544f049ab145fd901524b39295e 100644 (file)
--- a/include/deal.II/lac/petsc_parallel_vector.h
+++ b/include/deal.II/lac/petsc_parallel_vector.h
@@ -589,8 +589,8 @@ namespace internal
     template <typename> class ReinitHelper;
 
     /**
-     * A helper class internally used in linear_operator.h. Specialization for
-     * TrilinosWrappers::MPI::Vector.
+     * A helper class used internally in linear_operator.h. Specialization for
+     * PETScWrappers::MPI::Vector.
      */
     template<>
     class ReinitHelper<PETScWrappers::MPI::Vector>

diff --git a/include/deal.II/lac/petsc_vector.h b/include/deal.II/lac/petsc_vector.h
index ba4e3c4d5d13c688e18c31f3b675e5bfddb355c1..5e113f761857b77f633cd1f92246a9fd5ace82ba 100644 (file)
--- a/include/deal.II/lac/petsc_vector.h
+++ b/include/deal.II/lac/petsc_vector.h
@@ -398,7 +398,7 @@ namespace internal
     template <typename> class ReinitHelper;
 
     /**
-     * A helper class internally used in linear_operator.h. Specialization for
+     * A helper class used internally in linear_operator.h. Specialization for
      * PETScWrappers::Vector.
      */
     template<>

diff --git a/include/deal.II/lac/trilinos_block_vector.h b/include/deal.II/lac/trilinos_block_vector.h
index 2792388a3cdc6bea798f67c5783dbe2f93796e39..965bb147515fb777d37663dbe86c92d470533c6e 100644 (file)
--- a/include/deal.II/lac/trilinos_block_vector.h
+++ b/include/deal.II/lac/trilinos_block_vector.h
@@ -464,7 +464,7 @@ namespace internal
     template <typename> class ReinitHelper;
 
     /**
-     * A helper class internally used in linear_operator.h. Specialization for
+     * A helper class used internally in linear_operator.h. Specialization for
      * TrilinosWrappers::BlockVector.
      */
     template<>

diff --git a/include/deal.II/lac/trilinos_parallel_block_vector.h b/include/deal.II/lac/trilinos_parallel_block_vector.h
index 9213f88cce7f2d3c2d8d98ef32dc8b23fb368b86..f7a9cee6230fcaa4b29042a58a3a22878ecd3481 100644 (file)
--- a/include/deal.II/lac/trilinos_parallel_block_vector.h
+++ b/include/deal.II/lac/trilinos_parallel_block_vector.h
@@ -488,7 +488,7 @@ namespace internal
     template <typename> class ReinitHelper;
 
     /**
-     * A helper class internally used in linear_operator.h. Specialization for
+     * A helper class used internally in linear_operator.h. Specialization for
      * TrilinosWrappers::MPI::BlockVector.
      */
     template<>

diff --git a/include/deal.II/lac/trilinos_vector.h b/include/deal.II/lac/trilinos_vector.h
index 990311d95854affd3098c95c500235f581055fa5..2954a48f44ea6280703cf55eced38b2dad63d365 100644 (file)
--- a/include/deal.II/lac/trilinos_vector.h
+++ b/include/deal.II/lac/trilinos_vector.h
@@ -975,7 +975,7 @@ namespace internal
     template <typename> class ReinitHelper;
 
     /**
-     * A helper class internally used in linear_operator.h. Specialization for
+     * A helper class used internally in linear_operator.h. Specialization for
      * TrilinosWrappers::MPI::Vector.
      */
     template<>
@@ -1002,7 +1002,7 @@ namespace internal
     };
 
     /**
-     * A helper class internally used in linear_operator.h. Specialization for
+     * A helper class used internally in linear_operator.h. Specialization for
      * TrilinosWrappers::Vector.
      */
     template<>

diff --git a/source/lac/petsc_block_sparse_matrix.cc b/source/lac/petsc_block_sparse_matrix.cc
index bbe6817d027a471be287c56d6420a95dc553e249..23b5be69d918d9838e7a05de7b2dbb419e5b113c 100644 (file)
--- a/source/lac/petsc_block_sparse_matrix.cc
+++ b/source/lac/petsc_block_sparse_matrix.cc
@@ -66,7 +66,29 @@ namespace PETScWrappers
         }
   }
 
+  std::vector<BlockSparseMatrix::size_type >
+  BlockSparseMatrix::
+  locally_domain_sizes() const
+  {
+    std::vector< size_type > index_sets;
+
+    for ( unsigned int i=0; i<this->n_block_cols(); ++i)
+      index_sets.push_back(this->block(0,i).n());
+
+    return index_sets;
+  }
 
+  std::vector<BlockSparseMatrix::size_type >
+  BlockSparseMatrix::
+  locally_range_sizes() const
+  {
+    std::vector< size_type > index_sets;
+
+    for ( unsigned int i=0; i<this->n_block_rows(); ++i)
+      index_sets.push_back(this->block(i,0).m());
+
+    return index_sets;
+  }
 
   void
   BlockSparseMatrix::collect_sizes ()

diff --git a/source/lac/petsc_parallel_block_sparse_matrix.cc b/source/lac/petsc_parallel_block_sparse_matrix.cc
index 3c3555f2e096106af96c1d405ee24479725cdd53..751f9f896c86d6030a61302df930882e766b5345 100644 (file)
--- a/source/lac/petsc_parallel_block_sparse_matrix.cc
+++ b/source/lac/petsc_parallel_block_sparse_matrix.cc
@@ -126,7 +126,27 @@ namespace PETScWrappers
       BaseClass::collect_sizes ();
     }
 
+    std::vector< IndexSet >
+    BlockSparseMatrix::locally_owned_domain_indices() const
+    {
+      std::vector< IndexSet > index_sets;
+
+      for ( unsigned int i=0; i<this->n_block_cols(); ++i)
+        index_sets.push_back(this->block(0,i).locally_owned_domain_indices());
+
+      return index_sets;
+    }
 
+    std::vector< IndexSet >
+    BlockSparseMatrix::locally_owned_range_indices() const
+    {
+      std::vector< IndexSet > index_sets;
+
+      for ( unsigned int i=0; i<this->n_block_rows(); ++i)
+        index_sets.push_back(this->block(i,0).locally_owned_range_indices());
+
+      return index_sets;
+    }
 
     const MPI_Comm &
     BlockSparseMatrix::get_mpi_communicator () const

diff --git a/source/lac/petsc_parallel_sparse_matrix.cc b/source/lac/petsc_parallel_sparse_matrix.cc
index 043a8f5e0073d39416146a095b4476151688c808..0888d3582a9045d6e739c8210dd92500ce00481a 100644 (file)
--- a/source/lac/petsc_parallel_sparse_matrix.cc
+++ b/source/lac/petsc_parallel_sparse_matrix.cc
@@ -899,8 +899,11 @@ namespace PETScWrappers
       IS *cols = nullptr;
 
       ierr = MatGetSize (matrix, &n_rows, &n_cols);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
       ierr = MatGetOwnershipIS(matrix, rows, cols);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
       ierr = ISGetMinMax(*rows, &min, &max);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
 
       IndexSet locally_owned_domain_indices(n_rows);
       locally_owned_domain_indices.add_range(min, max);
@@ -917,8 +920,11 @@ namespace PETScWrappers
       IS *cols = nullptr;
 
       ierr = MatGetSize (matrix, &n_rows, &n_cols);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
       ierr = MatGetOwnershipIS(matrix, rows, cols);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
       ierr = ISGetMinMax(*cols, &min, &max);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
 
       IndexSet locally_owned_range_indices(n_cols);
       locally_owned_range_indices.add_range(min, max);

diff --git a/source/lac/petsc_sparse_matrix.cc b/source/lac/petsc_sparse_matrix.cc
index ec0eaa8471fc028b2a03de9204836806989927f8..750f54db0b6237b710b20cf2ddb4165bfc51f8db 100644 (file)
--- a/source/lac/petsc_sparse_matrix.cc
+++ b/source/lac/petsc_sparse_matrix.cc
@@ -313,6 +313,7 @@ namespace PETScWrappers
   {
     PetscInt m,n;
     PetscErrorCode ierr = MatGetSize(matrix, &m, &n);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     return m;
   }
@@ -322,6 +323,7 @@ namespace PETScWrappers
   {
     PetscInt m,n;
     PetscErrorCode ierr = MatGetSize(matrix, &m, &n);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     return n;
   }

diff --git a/tests/lac/linear_operator_09.cc b/tests/lac/linear_operator_09.cc
index 41c48cc1a7f4781eef614bf088cfd6d2bf157af8..43a3279f1be032312242d4ad933fd69d2c8d6118 100644 (file)
--- a/tests/lac/linear_operator_09.cc
+++ b/tests/lac/linear_operator_09.cc
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2015 by the deal.II authors
+// Copyright (C) 2016 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -14,7 +14,7 @@
 // ---------------------------------------------------------------------
 
 // Test that it is possible to instantiate a LinearOperator object for all
-// different kinds of Trilinos matrices and vectors
+// different kinds of PETSc matrices and vectors
 // TODO: A bit more tests...
 
 #include "../tests.h"
@@ -22,12 +22,21 @@
 
 #include <deal.II/lac/linear_operator.h>
 
+// Vectors:
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/petsc_sparse_matrix.h>
 
 #include <deal.II/lac/petsc_parallel_vector.h>
 #include <deal.II/lac/petsc_parallel_sparse_matrix.h>
 
+// Block Matrix and Vectors:
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_block_sparse_matrix.h>
+
+#include <deal.II/lac/petsc_parallel_block_vector.h>
+#include <deal.II/lac/petsc_parallel_block_sparse_matrix.h>
+
+
 using namespace dealii;
 
 int main(int argc, char *argv[])
@@ -47,6 +56,15 @@ int main(int argc, char *argv[])
     auto op_a  = linear_operator<PETScWrappers::MPI::Vector>(a);
   }
 
+  {
+    PETScWrappers::BlockSparseMatrix a;
+    auto op_a  = linear_operator<PETScWrappers::BlockVector>(a);
+  }
+
+  {
+    PETScWrappers::MPI::BlockSparseMatrix a;
+    auto op_a  = linear_operator<PETScWrappers::MPI::BlockVector>(a);
+  }
 
   deallog << "OK" << std::endl;