Fix TrilinosWrappers::PreconditionBlock* with no local rows

work around the ifpack error by pretending to use a point smoother on
processors without any local rows.

diff --git a/doc/news/changes/minor/20170704Heister b/doc/news/changes/minor/20170704Heister
new file mode 100644 (file)
index 0000000..7d12d2a
--- /dev/null
+++ b/doc/news/changes/minor/20170704Heister
@@ -0,0 +1,3 @@
+Fixed: TrilinosWrappers::PreconditionBlock* classes now work correctly if one of the processors does not own any rows.
+<br>
+(Timo Heister, 2017/07/04)

diff --git a/source/lac/trilinos_precondition.cc b/source/lac/trilinos_precondition.cc
index f0dfb4f38e425936c06ea098594f117997ae415c..5172593a14f539060269025bbd5e0303a2caa459 100644 (file)
--- a/source/lac/trilinos_precondition.cc
+++ b/source/lac/trilinos_precondition.cc
@@ -283,10 +283,15 @@ namespace TrilinosWrappers
   {
     // release memory before reallocation
     preconditioner.reset ();
+
+    // Block relaxation setup fails if we have no locally owned rows. As a work-around
+    // we just pretend to use point relaxation on those processors:
     preconditioner.reset (Ifpack().Create
-                          ("block relaxation",
-                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
-                           0));
+                          (
+                            (matrix.trilinos_matrix().NumMyRows()==0) ?
+                            "point relaxation" : "block relaxation",
+                            const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
+                            0));
 
     Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                     (preconditioner.get());
@@ -343,10 +348,15 @@ namespace TrilinosWrappers
                                      const AdditionalData &additional_data)
   {
     preconditioner.reset ();
+
+    // Block relaxation setup fails if we have no locally owned rows. As a work-around
+    // we just pretend to use point relaxation on those processors:
     preconditioner.reset (Ifpack().Create
-                          ("block relaxation",
-                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
-                           additional_data.overlap));
+                          (
+                            (matrix.trilinos_matrix().NumMyRows()==0) ?
+                            "point relaxation" : "block relaxation",
+                            const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
+                            additional_data.overlap));
 
     Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                     (preconditioner.get());
@@ -404,10 +414,15 @@ namespace TrilinosWrappers
                                     const AdditionalData &additional_data)
   {
     preconditioner.reset ();
+
+    // Block relaxation setup fails if we have no locally owned rows. As a work-around
+    // we just pretend to use point relaxation on those processors:
     preconditioner.reset (Ifpack().Create
-                          ("block relaxation",
-                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
-                           additional_data.overlap));
+                          (
+                            (matrix.trilinos_matrix().NumMyRows()==0) ?
+                            "point relaxation" : "block relaxation",
+                            const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
+                            additional_data.overlap));
 
     Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                     (preconditioner.get());

diff --git a/tests/trilinos/block_relax_no_rows.cc b/tests/trilinos/block_relax_no_rows.cc
new file mode 100644 (file)
index 0000000..91d4f2c
--- /dev/null
+++ b/tests/trilinos/block_relax_no_rows.cc
@@ -0,0 +1,77 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2013 - 2016 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Trilinos PreconditionBlock* used to fail if one processor has no locally
+// owned rows
+
+#include "../tests.h"
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_precondition.h>
+
+
+template <class Prec>
+void test()
+{
+  const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  const unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  // only proc 0 and 1 own rows
+  IndexSet rows(4*n_procs);
+  if (myid==0)
+    rows.add_range(0, 2*n_procs);
+  else if (myid==1)
+    rows.add_range(2*n_procs, 4*n_procs);
+  rows.compress();
+
+  TrilinosWrappers::MPI::Vector src (rows, MPI_COMM_WORLD),
+                   dst (rows, MPI_COMM_WORLD);
+
+  TrilinosWrappers::SparseMatrix mat(rows, rows, MPI_COMM_WORLD);
+  for (const auto &row : rows)
+    {
+      const unsigned int i = row;
+      mat.set(i, i, 100.);
+      for (unsigned int j=0; j<mat.n(); ++j)
+        if (i!=j)
+          mat.set (i,j, i*j*.5+.5);
+    }
+
+  for (unsigned int i=0; i<src.size(); ++i)
+    src(i) = i;
+
+  mat.compress (VectorOperation::insert);
+  src.compress (VectorOperation::insert);
+
+  Prec preconditioner;
+  typename Prec::AdditionalData data;
+  data.block_size = 4;
+
+  preconditioner.initialize(mat, data);
+  preconditioner.vmult(dst, src);
+
+  deallog << "dst: " << dst.l2_norm() << std::endl;
+}
+
+
+int main (int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+  MPILogInitAll log;
+
+  test<TrilinosWrappers::PreconditionBlockJacobi>();
+  test<TrilinosWrappers::PreconditionBlockSOR>();
+  test<TrilinosWrappers::PreconditionBlockSSOR>();
+}

diff --git a/tests/trilinos/block_relax_no_rows.mpirun=3.output b/tests/trilinos/block_relax_no_rows.mpirun=3.output
new file mode 100644 (file)
index 0000000..beaca9e
--- /dev/null
+++ b/tests/trilinos/block_relax_no_rows.mpirun=3.output
@@ -0,0 +1,14 @@
+
+DEAL:0::dst: 0.137153
+DEAL:0::dst: 0.112530
+DEAL:0::dst: 0.0911212
+
+DEAL:1::dst: 0.137153
+DEAL:1::dst: 0.112530
+DEAL:1::dst: 0.0911212
+
+
+DEAL:2::dst: 0.137153
+DEAL:2::dst: 0.112530
+DEAL:2::dst: 0.0911212
+