Allow iterating over TrilinosSparsityPattern in parallel and after compress

diff --git a/include/deal.II/lac/trilinos_sparsity_pattern.h b/include/deal.II/lac/trilinos_sparsity_pattern.h
index a61b358c21d063afc6fd49928ba8b325e69fe665..e2bf25e04d04b36447f77e4b32ee97b4b3e26a03 100644 (file)
--- a/include/deal.II/lac/trilinos_sparsity_pattern.h
+++ b/include/deal.II/lac/trilinos_sparsity_pattern.h
@@ -1288,6 +1288,7 @@ namespace TrilinosWrappers
     {}
 
 
+
     inline Iterator::Iterator(const Iterator &) = default;
 
 
@@ -1300,18 +1301,23 @@ namespace TrilinosWrappers
 
       ++accessor.a_index;
 
-      // If at end of line: do one
-      // step, then cycle until we
-      // find a row with a nonzero
-      // number of entries.
+      // If at end of line: do one step, then cycle until we find a row with a
+      // nonzero number of entries that is stored locally.
       if (accessor.a_index >= accessor.colnum_cache->size())
         {
           accessor.a_index = 0;
           ++accessor.a_row;
 
-          while ((accessor.a_row < accessor.sparsity_pattern->n_rows()) &&
-                 (accessor.sparsity_pattern->row_length(accessor.a_row) == 0))
-            ++accessor.a_row;
+          while (accessor.a_row < accessor.sparsity_pattern->n_rows())
+            {
+              const auto row_length =
+                accessor.sparsity_pattern->row_length(accessor.a_row);
+              if (row_length == 0 ||
+                  row_length == static_cast<SparsityPattern::size_type>(-1))
+                ++accessor.a_row;
+              else
+                break;
+            }
 
           accessor.visit_present_row();
         }
@@ -1376,7 +1382,8 @@ namespace TrilinosWrappers
   inline SparsityPattern::const_iterator
   SparsityPattern::begin() const
   {
-    return const_iterator(this, 0, 0);
+    const size_type first_valid_row = this->local_range().first;
+    return const_iterator(this, first_valid_row, 0);
   }
 
 

diff --git a/source/lac/trilinos_sparsity_pattern.cc b/source/lac/trilinos_sparsity_pattern.cc
index 6960fc94d7aa3d7a923fe5ea5688b0fd07e6d2b1..5a57541144e5698ab24a475ea55579e3aa52b302 100644 (file)
--- a/source/lac/trilinos_sparsity_pattern.cc
+++ b/source/lac/trilinos_sparsity_pattern.cc
@@ -41,12 +41,12 @@ namespace TrilinosWrappers
       if (this->a_row == sparsity_pattern->n_rows())
         {
           colnum_cache.reset();
-
           return;
         }
 
-      // otherwise first flush Trilinos caches
-      sparsity_pattern->compress();
+      // otherwise first flush Trilinos caches if necessary
+      if (!sparsity_pattern->is_compressed())
+        sparsity_pattern->compress();
 
       colnum_cache = std::make_shared<std::vector<size_type>>(
         sparsity_pattern->row_length(this->a_row));
@@ -815,14 +815,16 @@ namespace TrilinosWrappers
         ierr = graph->FillComplete(*column_space_map,
                                    static_cast<const Epetra_Map &>(
                                      graph->RangeMap()));
+        AssertThrow(ierr == 0, ExcTrilinosError(ierr));
       }
     else
-      ierr = graph->GlobalAssemble(*column_space_map,
-                                   static_cast<const Epetra_Map &>(
-                                     graph->RangeMap()),
-                                   true);
-
-    AssertThrow(ierr == 0, ExcTrilinosError(ierr));
+      {
+        ierr = graph->GlobalAssemble(*column_space_map,
+                                     static_cast<const Epetra_Map &>(
+                                       graph->RangeMap()),
+                                     true);
+        AssertThrow(ierr == 0, ExcTrilinosError(ierr));
+      }
 
     ierr = graph->OptimizeStorage();
     AssertThrow(ierr == 0, ExcTrilinosError(ierr));

diff --git a/tests/trilinos/sparsity_pattern_06.cc b/tests/trilinos/sparsity_pattern_06.cc
new file mode 100644 (file)
index 0000000..d4200b2
--- /dev/null
+++ b/tests/trilinos/sparsity_pattern_06.cc
@@ -0,0 +1,94 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test accessing the entries of TrilinosWrappers::SparsityPattern
+// after compress has been called.
+// The sparsity pattern used is the same as in sparse_matrix_add_03.
+
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+
+#include "../tests.h"
+
+
+void
+test()
+{
+  const unsigned int MyPID   = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  const unsigned int NumProc = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+
+  if (MyPID == 0)
+    deallog << "NumProc=" << NumProc << std::endl;
+
+  // create non-contiguous index set for NumProc > 3
+  dealii::IndexSet parallel_partitioning(NumProc * 2);
+
+  // non-contiguous
+  parallel_partitioning.add_index(MyPID);
+  parallel_partitioning.add_index(NumProc + MyPID);
+
+  // create sparsity pattern from parallel_partitioning
+
+  // The sparsity pattern corresponds to a [FE_DGQ<1>(p=0)]^2 FESystem,
+  // on a periodic triangulation in which each MPI process owns 2 cells,
+  // with reordered dofs by its components, such that the rows in the
+  // final matrix are locally not in a contiguous set.
+
+  dealii::TrilinosWrappers::SparsityPattern sp(parallel_partitioning,
+                                               MPI_COMM_WORLD,
+                                               2);
+
+  sp.add(MyPID, (NumProc + MyPID - 1) % NumProc);
+  sp.add(MyPID, MyPID);
+  sp.add(MyPID, (MyPID + 1) % NumProc);
+  sp.add(MyPID + NumProc, (NumProc + MyPID - 1) % NumProc + NumProc);
+  sp.add(MyPID + NumProc, MyPID + NumProc);
+  sp.add(MyPID + NumProc, (MyPID + 1) % NumProc + NumProc);
+
+  deallog << "before compress:" << std::endl;
+
+  for (const auto &el : sp)
+    {
+      deallog << "index: " << el.index() << " = "
+              << "(" << el.row() << " , " << el.column() << ")" << std::endl;
+    }
+
+  sp.compress();
+
+  deallog << "after compress:" << std::endl;
+  for (const auto &el : sp)
+    {
+      deallog << "index: " << el.index() << " = "
+              << "(" << el.row() << " , " << el.column() << ")" << std::endl;
+    }
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+  MPILogInitAll mpi_log;
+
+  test();
+}

diff --git a/tests/trilinos/sparsity_pattern_06.with_trilinos=true.with_mpi=true.mpirun=3.output b/tests/trilinos/sparsity_pattern_06.with_trilinos=true.with_mpi=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..6c15a08
--- /dev/null
+++ b/tests/trilinos/sparsity_pattern_06.with_trilinos=true.with_mpi=true.mpirun=3.output
@@ -0,0 +1,48 @@
+
+DEAL:0::NumProc=3
+DEAL:0::before compress:
+DEAL:0::index: 0 = (0 , 0)
+DEAL:0::index: 1 = (0 , 1)
+DEAL:0::index: 2 = (0 , 2)
+DEAL:0::index: 0 = (3 , 3)
+DEAL:0::index: 1 = (3 , 4)
+DEAL:0::index: 2 = (3 , 5)
+DEAL:0::after compress:
+DEAL:0::index: 0 = (0 , 0)
+DEAL:0::index: 1 = (0 , 1)
+DEAL:0::index: 2 = (0 , 2)
+DEAL:0::index: 0 = (3 , 3)
+DEAL:0::index: 1 = (3 , 4)
+DEAL:0::index: 2 = (3 , 5)
+
+DEAL:1::before compress:
+DEAL:1::index: 0 = (1 , 1)
+DEAL:1::index: 1 = (1 , 0)
+DEAL:1::index: 2 = (1 , 2)
+DEAL:1::index: 0 = (4 , 4)
+DEAL:1::index: 1 = (4 , 3)
+DEAL:1::index: 2 = (4 , 5)
+DEAL:1::after compress:
+DEAL:1::index: 0 = (1 , 1)
+DEAL:1::index: 1 = (1 , 0)
+DEAL:1::index: 2 = (1 , 2)
+DEAL:1::index: 0 = (4 , 4)
+DEAL:1::index: 1 = (4 , 3)
+DEAL:1::index: 2 = (4 , 5)
+
+
+DEAL:2::before compress:
+DEAL:2::index: 0 = (2 , 2)
+DEAL:2::index: 1 = (2 , 0)
+DEAL:2::index: 2 = (2 , 1)
+DEAL:2::index: 0 = (5 , 5)
+DEAL:2::index: 1 = (5 , 3)
+DEAL:2::index: 2 = (5 , 4)
+DEAL:2::after compress:
+DEAL:2::index: 0 = (2 , 2)
+DEAL:2::index: 1 = (2 , 0)
+DEAL:2::index: 2 = (2 , 1)
+DEAL:2::index: 0 = (5 , 5)
+DEAL:2::index: 1 = (5 , 3)
+DEAL:2::index: 2 = (5 , 4)
+