From 06f4882c7c51ce876464fa4eaff82c363e05b36c Mon Sep 17 00:00:00 2001
From: Timo Heister <timo.heister@gmail.com>
Date: Thu, 18 Apr 2013 22:33:08 +0000
Subject: [PATCH] better reinit for PETSc matrix, now identical to Trilinos

git-svn-id: https://svn.dealii.org/branches/branch_unify_linear_algebra@29336 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../lac/petsc_parallel_sparse_matrix.h        |  20 ++
 .../lac/petsc_parallel_sparse_matrix.cc       | 186 ++++++++++++++++++
 2 files changed, 206 insertions(+)

diff --git a/deal.II/include/deal.II/lac/petsc_parallel_sparse_matrix.h b/deal.II/include/deal.II/lac/petsc_parallel_sparse_matrix.h
index 29dded8591..7eb8b9d087 100644
--- a/deal.II/include/deal.II/lac/petsc_parallel_sparse_matrix.h
+++ b/deal.II/include/deal.II/lac/petsc_parallel_sparse_matrix.h
@@ -368,6 +368,18 @@ namespace PETScWrappers
                    const unsigned int               this_process,
                    const bool                       preset_nonzero_locations = true);
 
+      /**
+       * Create a matrix where the size() of the IndexSets determine the global
+       * number of rows and columns and the entries of the IndexSet give
+       * the rows and columns for the calling processor.
+       * Note that only contiguous IndexSets are supported.
+       */
+      template <typename SparsityType>
+      void reinit (const IndexSet & local_rows,
+                  const IndexSet & local_columns,
+                   const SparsityType         &sparsity_pattern,
+                   const MPI_Comm                  &communicator);
+
       /**
        * Return a reference to the MPI
        * communicator object in use with
@@ -472,6 +484,14 @@ namespace PETScWrappers
                       const unsigned int               this_process,
                       const bool                       preset_nonzero_locations);
 
+      /**
+       * Same as previous functions.
+       */
+      template <typename SparsityType>
+      void do_reinit (const IndexSet & local_rows,
+                      const IndexSet & local_columns,
+                       const SparsityType         &sparsity_pattern);
+
       /**
        *  To allow calling protected
        *  prepare_add() and
diff --git a/deal.II/source/lac/petsc_parallel_sparse_matrix.cc b/deal.II/source/lac/petsc_parallel_sparse_matrix.cc
index b90774a75f..405bc0c2bf 100644
--- a/deal.II/source/lac/petsc_parallel_sparse_matrix.cc
+++ b/deal.II/source/lac/petsc_parallel_sparse_matrix.cc
@@ -178,7 +178,27 @@ namespace PETScWrappers
                  preset_nonzero_locations);
     }
 
+    template <typename SparsityType>
+    void
+    SparseMatrix::
+    reinit (const IndexSet & local_rows,
+        const IndexSet & local_columns,
+        const SparsityType &sparsity_pattern,
+        const MPI_Comm &communicator)
+    {
+      this->communicator = communicator;
 
+      // get rid of old matrix and generate a
+      // new one
+#if DEAL_II_PETSC_VERSION_LT(3,2,0)
+      const int ierr = MatDestroy (matrix);
+#else
+      const int ierr = MatDestroy (&matrix);
+#endif
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+      do_reinit (local_rows, local_columns, sparsity_pattern);
+    }
 
     void
     SparseMatrix::do_reinit (const unsigned int m,
@@ -308,6 +328,159 @@ namespace PETScWrappers
     }
 
 
+    template <typename SparsityType>
+    void
+    SparseMatrix::
+    do_reinit (const IndexSet & local_rows,
+        const IndexSet & local_columns,
+        const SparsityType         &sparsity_pattern)
+    {
+      Assert(sparsity_pattern.n_rows()==local_rows.size(),
+          ExcMessage("SparsityPattern and IndexSet have different number of rows"));
+      Assert(sparsity_pattern.n_cols()==local_columns.size(),
+          ExcMessage("SparsityPattern and IndexSet have different number of columns"));
+      Assert(local_rows.is_contiguous() && local_columns.is_contiguous(),
+          ExcMessage("PETSc only supports contiguous row/column ranges"));
+
+
+
+
+            // create the matrix. We do not set row length but set the
+            // correct SparsityPattern later.
+            int ierr;
+
+            ierr = MatCreate(communicator,&matrix);
+            AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+            ierr = MatSetSizes(matrix,
+                                local_rows.n_elements(),
+                                local_columns.n_elements(),
+                               sparsity_pattern.n_rows(),
+                               sparsity_pattern.n_cols());
+            AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+            ierr = MatSetType(matrix,MATMPIAIJ);
+            AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+
+            // next preset the exact given matrix
+            // entries with zeros. this doesn't avoid any
+            // memory allocations, but it at least
+            // avoids some searches later on. the
+            // key here is that we can use the
+            // matrix set routines that set an
+            // entire row at once, not a single
+            // entry at a time
+            //
+            // for the usefulness of this option
+            // read the documentation of this
+            // class.
+            //if (preset_nonzero_locations == true)
+              {
+                // MatMPIAIJSetPreallocationCSR
+                // can be used to allocate the sparsity
+                // pattern of a matrix
+
+              const PetscInt local_row_start = local_rows.nth_index_in_set(0);
+              const PetscInt local_col_start = local_columns.nth_index_in_set(0);
+              const PetscInt
+                local_row_end = local_row_start + local_rows.n_elements();
+
+
+                // first set up the column number
+                // array for the rows to be stored
+                // on the local processor. have one
+                // dummy entry at the end to make
+                // sure petsc doesn't read past the
+                // end
+                std::vector<PetscInt>
+
+                rowstart_in_window (local_row_end - local_row_start + 1, 0),
+                                   colnums_in_window;
+                {
+                  unsigned int n_cols = 0;
+                  for (PetscInt i=local_row_start; i<local_row_end; ++i)
+                    {
+                      const PetscInt row_length = sparsity_pattern.row_length(i);
+                      rowstart_in_window[i+1-local_row_start]
+                        = rowstart_in_window[i-local_row_start] + row_length;
+                      n_cols += row_length;
+                    }
+                  colnums_in_window.resize (n_cols+1, -1);
+                }
+
+                // now copy over the information
+                // from the sparsity pattern.
+                {
+                  PetscInt* ptr = & colnums_in_window[0];
+
+                  for (PetscInt i=local_row_start; i<local_row_end; ++i)
+                    {
+                      typename SparsityType::row_iterator
+                      row_start = sparsity_pattern.row_begin(i),
+                      row_end = sparsity_pattern.row_end(i);
+
+                      std::copy(row_start, row_end, ptr);
+                      ptr += row_end - row_start;
+                    }
+                }
+
+
+                // then call the petsc function
+                // that summarily allocates these
+                // entries:
+                MatMPIAIJSetPreallocationCSR (matrix,
+                                              &rowstart_in_window[0],
+                                              &colnums_in_window[0],
+                                              0);
+                compress ();
+
+
+                // Tell PETSc that we are not
+                // planning on adding new entries
+                // to the matrix. Generate errors
+                // in debug mode.
+                int ierr;
+      #if DEAL_II_PETSC_VERSION_LT(3,0,0)
+      #ifdef DEBUG
+                ierr = MatSetOption (matrix, MAT_NEW_NONZERO_LOCATION_ERR);
+                AssertThrow (ierr == 0, ExcPETScError(ierr));
+      #else
+                ierr = MatSetOption (matrix, MAT_NO_NEW_NONZERO_LOCATIONS);
+                AssertThrow (ierr == 0, ExcPETScError(ierr));
+      #endif
+      #else
+      #ifdef DEBUG
+                ierr = MatSetOption (matrix, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
+                AssertThrow (ierr == 0, ExcPETScError(ierr));
+      #else
+                ierr = MatSetOption (matrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
+                AssertThrow (ierr == 0, ExcPETScError(ierr));
+      #endif
+      #endif
+
+                // Tell PETSc to keep the
+                // SparsityPattern entries even if
+                // we delete a row with
+                // clear_rows() which calls
+                // MatZeroRows(). Otherwise one can
+                // not write into that row
+                // afterwards.
+      #if DEAL_II_PETSC_VERSION_LT(3,0,0)
+                ierr = MatSetOption (matrix, MAT_KEEP_ZEROED_ROWS);
+                AssertThrow (ierr == 0, ExcPETScError(ierr));
+      #elif DEAL_II_PETSC_VERSION_LT(3,1,0)
+                ierr = MatSetOption (matrix, MAT_KEEP_ZEROED_ROWS, PETSC_TRUE);
+                AssertThrow (ierr == 0, ExcPETScError(ierr));
+      #else
+                ierr = MatSetOption (matrix, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
+                AssertThrow (ierr == 0, ExcPETScError(ierr));
+      #endif
+
+              }
+
+    }
+
 
     template <typename SparsityType>
     void
@@ -608,6 +781,13 @@ namespace PETScWrappers
                           const unsigned int,
                           const bool);
 
+    template void
+    SparseMatrix::
+        reinit (const IndexSet &,
+            const IndexSet &,
+            const CompressedSimpleSparsityPattern &,
+            const MPI_Comm                  &);
+
     template void
     SparseMatrix::do_reinit (const SparsityPattern &,
                              const std::vector<unsigned int> &,
@@ -627,6 +807,12 @@ namespace PETScWrappers
                              const unsigned int ,
                              const bool);
 
+    template void
+        SparseMatrix::
+        do_reinit (const IndexSet &,
+            const IndexSet &,
+            const CompressedSimpleSparsityPattern &);
+
 
     PetscScalar
     SparseMatrix::matrix_norm_square (const Vector &v) const
-- 
2.39.5