Fix a bug in distributed::Vector::import

The import function would fail if the ReadWriteVector would only contains
ghosted values

diff --git a/include/deal.II/lac/la_parallel_vector.templates.h b/include/deal.II/lac/la_parallel_vector.templates.h
index d4d233d820e8686c6e6643b7720d81df27f9a0e4..28984635529bccd350d8b4a00cd51e3e6e79827c 100644 (file)
--- a/include/deal.II/lac/la_parallel_vector.templates.h
+++ b/include/deal.II/lac/la_parallel_vector.templates.h
@@ -36,6 +36,49 @@ namespace LinearAlgebra
 {
   namespace distributed
   {
+    namespace internal
+    {
+      // In the import_from_ghosted_array_finish we might need to calculate the
+      // maximal and minimal value for the given number type, which is not
+      // straight forward for complex numbers. Therefore, comparison of complex
+      // numbers is prohibited and throws an exception.
+      template <typename Number>
+      Number
+      get_min(const Number a, const Number b)
+      {
+        return std::min(a, b);
+      }
+
+      template <typename Number>
+      std::complex<Number>
+      get_min(const std::complex<Number> a, const std::complex<Number>)
+      {
+        AssertThrow(false,
+                    ExcMessage("VectorOperation::min not "
+                               "implemented for complex numbers"));
+        return a;
+      }
+
+      template <typename Number>
+      Number
+      get_max(const Number a, const Number b)
+      {
+        return std::max(a, b);
+      }
+
+      template <typename Number>
+      std::complex<Number>
+      get_max(const std::complex<Number> a, const std::complex<Number>)
+      {
+        AssertThrow(false,
+                    ExcMessage("VectorOperation::max not "
+                               "implemented for complex numbers"));
+        return a;
+      }
+    } // namespace internal
+
+
+
     template <typename Number>
     void
     Vector<Number>::clear_mpi_requests()
@@ -669,8 +712,7 @@ namespace LinearAlgebra
           // then create the communication pattern
           IndexSet ghost_indices(V.get_stored_elements());
           ghost_indices.subtract_set(locally_owned_elem);
-          IndexSet local_indices(V.get_stored_elements());
-          local_indices.subtract_set(ghost_indices);
+          IndexSet local_indices(locally_owned_elem);
           comm_pattern = std::make_shared<Utilities::MPI::Partitioner>(
             local_indices, ghost_indices, get_mpi_communicator());
         }
@@ -683,28 +725,61 @@ namespace LinearAlgebra
                       ExcMessage("The communication pattern is not of type "
                                  "Utilities::MPI::Partitioner."));
         }
+      IndexSet ghost_indices(V.get_stored_elements());
+      ghost_indices.subtract_set(locally_owned_elem);
+      IndexSet       local_indices(locally_owned_elem);
       Vector<Number> tmp_vector(comm_pattern);
+      std::copy(begin(), end(), tmp_vector.begin());
 
       // fill entries from ReadWriteVector into the distributed vector,
       // including ghost entries. this is not really efficient right now
       // because indices are translated twice, once by nth_index_in_set(i) and
       // once for operator() of tmp_vector
-      const IndexSet &v_stored = V.get_stored_elements();
-      for (size_type i = 0; i < v_stored.n_elements(); ++i)
-        tmp_vector(v_stored.nth_index_in_set(i)) = V.local_element(i);
+      const IndexSet &v_stored     = V.get_stored_elements();
+      const size_type v_n_elements = v_stored.n_elements();
+      switch (operation)
+        {
+          case VectorOperation::insert:
+            {
+              for (size_type i = 0; i < v_n_elements; ++i)
+                tmp_vector(v_stored.nth_index_in_set(i)) = V.local_element(i);
 
-      tmp_vector.compress(operation);
+              break;
+            }
+          case VectorOperation::add:
+            {
+              for (size_type i = 0; i < v_n_elements; ++i)
+                tmp_vector(v_stored.nth_index_in_set(i)) += V.local_element(i);
 
-      // Copy the local elements of tmp_vector to the right place in val
-      IndexSet tmp_index_set = tmp_vector.locally_owned_elements();
-      for (size_type i = 0; i < tmp_index_set.n_elements(); ++i)
-        {
-          values[locally_owned_elem.index_within_set(
-            tmp_index_set.nth_index_in_set(i))] = tmp_vector.local_element(i);
-        }
-    }
+              break;
+            }
+          case VectorOperation::min:
+            {
+              for (size_type i = 0; i < v_n_elements; ++i)
+                tmp_vector(v_stored.nth_index_in_set(i)) =
+                  internal::get_min(tmp_vector(v_stored.nth_index_in_set(i)),
+                                    V.local_element(i));
+
+              break;
+            }
+          case VectorOperation::max:
+            {
+              for (size_type i = 0; i < v_n_elements; ++i)
+                tmp_vector(v_stored.nth_index_in_set(i)) =
+                  internal::get_max(tmp_vector(v_stored.nth_index_in_set(i)),
+                                    V.local_element(i));
 
+              break;
+            }
+          default:
+            {
+              Assert(false, ExcMessage("This operation is not supported."));
+            }
+        }
+      tmp_vector.compress(operation);
 
+      std::copy(tmp_vector.begin(), tmp_vector.end(), begin());
+    }
 
     template <typename Number>
     void

diff --git a/tests/mpi/parallel_vector_24.cc b/tests/mpi/parallel_vector_24.cc
new file mode 100644 (file)
index 0000000..9d510bc
--- /dev/null
+++ b/tests/mpi/parallel_vector_24.cc
@@ -0,0 +1,143 @@
+// ---------------------------------------------------------------------
+//
+// 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 based on parallel_vector_03a but using import function. This test was
+// not working because we would only set ghost elements and no local elements
+// which resulted in an error during the import
+
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/read_write_vector.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+void
+check(const unsigned int                                myid,
+      const LinearAlgebra::distributed::Vector<double> &v)
+{
+  v.print(std::cout);
+  if (myid == 0)
+    {
+      AssertThrow(v(10) == 10.0, ExcInternalError());
+      AssertThrow(v(11) == 0., ExcInternalError());
+      AssertThrow(v(12) == 0., ExcInternalError());
+      AssertThrow(v(14) == 14., ExcInternalError());
+
+      AssertThrow(v(5) == 55., ExcInternalError());
+    }
+  else
+    {
+      AssertThrow(v(4) == 0., ExcInternalError());
+      AssertThrow(v(5) == 55., ExcInternalError());
+      AssertThrow(v(6) == 66., ExcInternalError());
+    }
+}
+
+
+void
+test()
+{
+  unsigned int myid    = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  unsigned int numproc = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+
+  Assert(numproc == 2, ExcNotImplemented());
+
+  const unsigned int size = 20;
+  IndexSet           local_owned(size);
+  IndexSet           local_nonzero(size);
+  IndexSet           local_relevant(size);
+  if (myid == 0)
+    {
+      local_owned.add_range(0, 10);
+      local_nonzero.add_range(5, 10);
+      local_relevant = local_owned;
+      local_relevant.add_range(10, 13);
+      local_relevant.add_range(14, 15);
+    }
+  else
+    {
+      local_owned.add_range(10, size);
+      local_nonzero.add_range(10, 11);
+      local_nonzero.add_range(13, 15);
+      local_relevant = local_owned;
+      local_relevant.add_range(4, 7);
+    }
+
+  LinearAlgebra::distributed::Vector<double> v(local_owned,
+                                               local_relevant,
+                                               MPI_COMM_WORLD);
+  v = 0.;
+
+  // set local values
+  for (unsigned int i = 0; i < local_nonzero.n_elements(); i++)
+    v(local_nonzero.nth_index_in_set(i)) = local_nonzero.nth_index_in_set(i);
+
+  // set value from processor which does not own it:
+  v(5) = 55.;
+  v.compress(VectorOperation::insert);
+
+  // add to value from processor which has it as a ghost
+  IndexSet                               ghost_set(size);
+  LinearAlgebra::ReadWriteVector<double> rw_vector;
+  if (myid == 1)
+    {
+      ghost_set.add_index(6);
+      rw_vector.reinit(ghost_set);
+      rw_vector(6) = 60;
+    }
+  else
+    {
+      rw_vector.reinit(ghost_set);
+    }
+  v.import(rw_vector, VectorOperation::add); // 60 + 6
+  // compress(insert) used to leave ghosts un-touched which resulted in
+  // the wrong 55+55 for this compress(add) operation.
+
+  v.update_ghost_values();
+
+  check(myid, v);
+
+  if (myid == 0)
+    deallog << "OK" << std::endl;
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+  unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  deallog.push(Utilities::int_to_string(myid));
+
+  if (myid == 0)
+    {
+      std::ofstream logfile("output");
+      deallog.attach(logfile);
+      deallog << std::setprecision(4);
+
+      test();
+    }
+  else
+    test();
+}

diff --git a/tests/mpi/parallel_vector_24.mpirun=2.output b/tests/mpi/parallel_vector_24.mpirun=2.output
new file mode 100644 (file)
index 0000000..be8d055
--- /dev/null
+++ b/tests/mpi/parallel_vector_24.mpirun=2.output
@@ -0,0 +1,2 @@
+
+DEAL:0::OK