Enable PETScWrappers::MPI::SparseMatrix interface to account for off-diagonal element...

diff --git a/include/deal.II/lac/petsc_parallel_sparse_matrix.h b/include/deal.II/lac/petsc_parallel_sparse_matrix.h
index a28dc035ff3ed5c5332cb6f1e1ab28f974dde360..1bf0bc01cd409a08dcbb81fd8df3e1b6105173c6 100644 (file)
--- a/include/deal.II/lac/petsc_parallel_sparse_matrix.h
+++ b/include/deal.II/lac/petsc_parallel_sparse_matrix.h
@@ -166,11 +166,15 @@ namespace PETScWrappers
        * Create a sparse matrix of
        * dimensions @p m times @p n, with
        * an initial guess of @p
-       * n_nonzero_per_row nonzero elements
-       * per row. PETSc is able to cope
-       * with the situation that more than
-       * this number of elements are later
-       * allocated for a row, but this
+       * n_nonzero_per_row and @p
+       * n_offdiag_nonzero_per_row nonzero 
+       * elements per row (see documentation 
+       * of the MatCreateAIJ PETSc function 
+       * for more information about these 
+       * parameters). PETSc is able to 
+       * cope with the situation that more 
+       * than this number of elements are 
+       * later allocated for a row, but this
        * involves copying data, and is thus
        * expensive.
        *
@@ -198,14 +202,17 @@ namespace PETScWrappers
                     const size_type  local_rows,
                     const size_type  local_columns,
                     const size_type  n_nonzero_per_row,
-                    const bool       is_symmetric = false);
+                    const bool       is_symmetric = false,
+                    const size_type  n_offdiag_nonzero_per_row = 0);
 
       /**
        * Initialize a rectangular matrix
        * with @p m rows and @p n columns.
        * The maximal number of nonzero
-       * entries for each row separately is
-       * given by the @p row_lengths array.
+       * entries for diagonal and off-
+       * diagonal blocks of each row is
+       * given by the @p row_lengths and 
+       * @p offdiag_row_lengths arrays.
        *
        * For the meaning of the @p
        * local_row and @p local_columns
@@ -239,7 +246,8 @@ namespace PETScWrappers
                     const size_type               local_rows,
                     const size_type               local_columns,
                     const std::vector<size_type> &row_lengths,
-                    const bool                    is_symmetric = false);
+                    const bool                    is_symmetric = false,
+                    const std::vector<size_type> &offdiag_row_lengths = std::vector<size_type>());
 
       /**
        * Initialize using the given
@@ -328,7 +336,8 @@ namespace PETScWrappers
                    const size_type local_rows,
                    const size_type local_columns,
                    const size_type n_nonzero_per_row,
-                   const bool      is_symmetric = false);
+                   const bool      is_symmetric = false,
+                   const size_type n_offdiag_nonzero_per_row = 0);
 
       /**
        * Throw away the present matrix and
@@ -344,7 +353,8 @@ namespace PETScWrappers
                    const size_type               local_rows,
                    const size_type               local_columns,
                    const std::vector<size_type> &row_lengths,
-                   const bool                    is_symmetric = false);
+                   const bool                    is_symmetric = false,
+                  const std::vector<size_type> &offdiag_row_lengths = std::vector<size_type>());
 
       /**
        * Initialize using the given
@@ -482,7 +492,8 @@ namespace PETScWrappers
                       const size_type local_rows,
                       const size_type local_columns,
                       const size_type n_nonzero_per_row,
-                      const bool      is_symmetric = false);
+                      const bool      is_symmetric = false,
+                      const size_type n_offdiag_nonzero_per_row = 0);
 
       /**
        * Same as previous function.
@@ -492,7 +503,8 @@ namespace PETScWrappers
                       const size_type               local_rows,
                       const size_type               local_columns,
                       const std::vector<size_type> &row_lengths,
-                      const bool                    is_symmetric = false);
+                      const bool                    is_symmetric = false,
+                     const std::vector<size_type> &offdiag_row_lengths = std::vector<size_type>());
 
       /**
        * Same as previous functions.

diff --git a/source/lac/petsc_parallel_sparse_matrix.cc b/source/lac/petsc_parallel_sparse_matrix.cc
index 7a842188a048dfbc4d1b718245367e3cd10c0053..3e59ffe75ef240753e2f7f6d8b56c25915e36e70 100644 (file)
--- a/source/lac/petsc_parallel_sparse_matrix.cc
+++ b/source/lac/petsc_parallel_sparse_matrix.cc
@@ -60,12 +60,14 @@ namespace PETScWrappers
                                 const size_type  local_rows,
                                 const size_type  local_columns,
                                 const size_type  n_nonzero_per_row,
-                                const bool       is_symmetric)
+                                const bool       is_symmetric,
+                                const size_type  n_offdiag_nonzero_per_row)
       :
       communicator (communicator)
     {
       do_reinit (m, n, local_rows, local_columns,
-                 n_nonzero_per_row, is_symmetric);
+                 n_nonzero_per_row, is_symmetric,
+                 n_offdiag_nonzero_per_row);
     }
 
 
@@ -76,12 +78,13 @@ namespace PETScWrappers
                                 const size_type               local_rows,
                                 const size_type               local_columns,
                                 const std::vector<size_type> &row_lengths,
-                                const bool                    is_symmetric)
+                                const bool                    is_symmetric,
+                                const std::vector<size_type> &offdiag_row_lengths)
       :
       communicator (communicator)
     {
       do_reinit (m, n, local_rows, local_columns,
-                 row_lengths, is_symmetric);
+                 row_lengths, is_symmetric, offdiag_row_lengths);
     }
 
 
@@ -130,7 +133,8 @@ namespace PETScWrappers
                           const size_type  local_rows,
                           const size_type  local_columns,
                           const size_type  n_nonzero_per_row,
-                          const bool       is_symmetric)
+                          const bool       is_symmetric,
+                          const size_type  n_offdiag_nonzero_per_row)
     {
       this->communicator = communicator;
 
@@ -144,7 +148,8 @@ namespace PETScWrappers
       AssertThrow (ierr == 0, ExcPETScError(ierr));
 
       do_reinit (m, n, local_rows, local_columns,
-                 n_nonzero_per_row, is_symmetric);
+                 n_nonzero_per_row, is_symmetric,
+                 n_offdiag_nonzero_per_row);
     }
 
 
@@ -156,7 +161,8 @@ namespace PETScWrappers
                           const size_type               local_rows,
                           const size_type               local_columns,
                           const std::vector<size_type> &row_lengths,
-                          const bool                    is_symmetric)
+                          const bool                    is_symmetric,
+                          const std::vector<size_type> &offdiag_row_lengths)
     {
       this->communicator = communicator;
 
@@ -169,7 +175,8 @@ namespace PETScWrappers
 #endif
       AssertThrow (ierr == 0, ExcPETScError(ierr));
 
-      do_reinit (m, n, local_rows, local_columns, row_lengths, is_symmetric);
+      do_reinit (m, n, local_rows, local_columns, 
+                 row_lengths, is_symmetric, offdiag_row_lengths);
     }
 
 
@@ -228,7 +235,8 @@ namespace PETScWrappers
                              const size_type local_rows,
                              const size_type local_columns,
                              const size_type n_nonzero_per_row,
-                             const bool      is_symmetric)
+                             const bool      is_symmetric,
+                             const size_type n_offdiag_nonzero_per_row)
     {
       Assert (local_rows <= m, ExcLocalRowsTooLarge (local_rows, m));
 
@@ -242,14 +250,16 @@ namespace PETScWrappers
         = MatCreateMPIAIJ (communicator,
                            local_rows, local_columns,
                            m, n,
-                           n_nonzero_per_row, 0, 0, 0,
+                           n_nonzero_per_row, 0, 
+                           n_offdiag_nonzero_per_row, 0,
                            &matrix);
 #else
       ierr
         = MatCreateAIJ (communicator,
                         local_rows, local_columns,
                         m, n,
-                        n_nonzero_per_row, 0, 0, 0,
+                        n_nonzero_per_row, 0, 
+                        n_offdiag_nonzero_per_row, 0,
                         &matrix);
       AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -280,7 +290,8 @@ namespace PETScWrappers
                              const size_type local_rows,
                              const size_type local_columns,
                              const std::vector<size_type> &row_lengths,
-                             const bool      is_symmetric)
+                             const bool      is_symmetric,
+                             const std::vector<size_type> &offdiag_row_lengths)
     {
       Assert (local_rows <= m, ExcLocalRowsTooLarge (local_rows, m));
 
@@ -307,6 +318,8 @@ namespace PETScWrappers
       // tricks with conversions of pointers
       const std::vector<PetscInt> int_row_lengths (row_lengths.begin(),
                                                    row_lengths.end());
+      const std::vector<PetscInt> int_offdiag_row_lengths (offdiag_row_lengths.begin(),
+                                                           offdiag_row_lengths.end());
 
 //TODO: There must be a significantly better way to provide information about the off-diagonal blocks of the matrix. this way, petsc keeps allocating tiny chunks of memory, and gets completely hung up over this
 
@@ -317,14 +330,16 @@ namespace PETScWrappers
         = MatCreateMPIAIJ (communicator,
                            local_rows, local_columns,
                            m, n,
-                           0, &int_row_lengths[0], 0, 0,
+                           0, &int_row_lengths[0], 
+                           0, offdiag_row_lengths.size() ? &int_offdiag_row_lengths[0] : 0,
                            &matrix);
 #else
       ierr
         = MatCreateAIJ (communicator,
                         local_rows, local_columns,
                         m, n,
-                        0, &int_row_lengths[0], 0, 0,
+                        0, &int_row_lengths[0], 
+                        0, offdiag_row_lengths.size() ? &int_offdiag_row_lengths[0] : 0,
                         &matrix);
       AssertThrow (ierr == 0, ExcPETScError(ierr));