]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduce min and max to VectorOperation for parallel vectors
authorSvenja Schoeder <schoeder@lnm.mw.tum.de>
Mon, 18 Jun 2018 14:24:22 +0000 (16:24 +0200)
committerSvenja Schoeder <schoeder@lnm.mw.tum.de>
Fri, 29 Jun 2018 11:19:45 +0000 (13:19 +0200)
Add tests for min max VectorOperation

add changelog

disallow complex numbers for VectorOperation::min and VectorOperation::max

testing of VectorOperation::max if several processors have same ghost element but with different values

fixed typo on cases

12 files changed:
doc/news/changes/minor/20180618SvenjaSchoeder [new file with mode: 0644]
include/deal.II/base/partitioner.templates.h
include/deal.II/lac/la_parallel_vector.h
include/deal.II/lac/read_write_vector.templates.h
include/deal.II/lac/vector_operation.h
source/lac/cuda_vector.cu
source/lac/trilinos_epetra_vector.cc
source/lac/trilinos_sparse_matrix.cc
source/lac/trilinos_vector.cc
tests/mpi/parallel_vector_22.cc [new file with mode: 0644]
tests/mpi/parallel_vector_22.mpirun=10.output [new file with mode: 0644]
tests/mpi/parallel_vector_22.mpirun=4.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20180618SvenjaSchoeder b/doc/news/changes/minor/20180618SvenjaSchoeder
new file mode 100644 (file)
index 0000000..7c5f98b
--- /dev/null
@@ -0,0 +1,6 @@
+New: The LinearAlgebra::distributed::Vector::compress() function has
+gained two new combine operations VectorOperation::min and
+VectorOperation::max to define the entry at the owning processor as
+the minimum/maximum of its own value and the one coming from a ghost.
+<br>
+(Svenja Schoeder, 2018/06/18)
index 2b5c8bfb1618e09f19d14a3d955e09e112809c37..87b5c601ff04611091ae095d9422f30c52784234 100644 (file)
@@ -329,6 +329,43 @@ namespace Utilities
         return a;
       }
 
+      // In the import_from_ghosted_array_finish we need to calculate maximal
+      // and minimal value on number types, which is not straight forward for
+      // complex numbers. Therfore, comparison of complex numbers is
+      // prohibited and throw an assert.
+      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 max 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::min max not "
+                               "implemented for complex numbers"));
+        return a;
+      }
     } // namespace internal
 
 
@@ -403,11 +440,29 @@ namespace Utilities
           // local values. For insert, nothing is done here (but in debug mode
           // we assert that the specified value is either zero or matches with
           // the ones already present
-          if (vector_operation != dealii::VectorOperation::insert)
+          if (vector_operation == dealii::VectorOperation::add)
             for (; my_imports != import_indices_data.end(); ++my_imports)
               for (unsigned int j = my_imports->first; j < my_imports->second;
                    j++)
                 locally_owned_array[j] += *read_position++;
+          else if (vector_operation == dealii::VectorOperation::min)
+            for (; my_imports != import_indices_data.end(); ++my_imports)
+              for (unsigned int j = my_imports->first; j < my_imports->second;
+                   j++)
+                {
+                  locally_owned_array[j] =
+                    internal::get_min(*read_position, locally_owned_array[j]);
+                  read_position++;
+                }
+          else if (vector_operation == dealii::VectorOperation::max)
+            for (; my_imports != import_indices_data.end(); ++my_imports)
+              for (unsigned int j = my_imports->first; j < my_imports->second;
+                   j++)
+                {
+                  locally_owned_array[j] =
+                    internal::get_max(*read_position, locally_owned_array[j]);
+                  read_position++;
+                }
           else
             for (; my_imports != import_indices_data.end(); ++my_imports)
               for (unsigned int j = my_imports->first; j < my_imports->second;
index 819d88cb5ea4302f4e6cd465f7ca72aa513f4590..0a7fbc9dfd256644a1e114c800876e8872c53120 100644 (file)
@@ -402,7 +402,7 @@ namespace LinearAlgebra
        * @ref GlossCompress "Compressing distributed vectors and matrices"
        * in the glossary.
        *
-       * There are two variants for this function. If called with argument @p
+       * There are four variants for this function. If called with argument @p
        * VectorOperation::add adds all the data accumulated in ghost elements
        * to the respective elements on the owning processor and clears the
        * ghost array afterwards. If called with argument @p
@@ -416,6 +416,13 @@ namespace LinearAlgebra
        * actually consistent between processors, i.e., whenever a non-zero
        * ghost element is found, it is compared to the value on the owning
        * processor and an exception is thrown if these elements do not agree.
+       * If called with VectorOperation::min or VectorOperation::max the
+       * minimum or maximum on all elements across the processors is set. If
+       * a processor does not set values for the ghosts, the value 0.0 is
+       * assumed, which can be the minimal or maximal value across all
+       * processors. The operations min and max are reasonable if all
+       * processors set relevant values or set values that do not alter the
+       * results, e.g. small values for the max operation.
        */
       virtual void
       compress(::dealii::VectorOperation::values operation) override;
index d577e4c7c73e15af5ad38050fa78aaeac1f9cf53..62608caf1089816c7e918ad4ce7574c57c0254b8 100644 (file)
@@ -250,7 +250,46 @@ namespace LinearAlgebra
     return *this;
   }
 
+  namespace internal
+  {
+    // In the import_from_ghosted_array_finish we need to calculate maximal
+    // and minimal value on number types, which is not straight forward for
+    // complex numbers. Therfore, comparison of complex numbers is
+    // prohibited and throw an assert.
+    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 max 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::min max not "
+                             "implemented for complex numbers"));
+      return a;
+    }
+  } // namespace internal
 
   template <typename Number>
   void
@@ -288,6 +327,16 @@ namespace LinearAlgebra
     if (operation == VectorOperation::add)
       for (size_type i = 0; i < stored.n_elements(); ++i)
         local_element(i) += tmp_vector(stored.nth_index_in_set(i));
+    else if (operation == VectorOperation::min)
+      for (size_type i = 0; i < stored.n_elements(); ++i)
+        local_element(i) =
+          internal::get_min(tmp_vector(stored.nth_index_in_set(i)),
+                            local_element(i));
+    else if (operation == VectorOperation::max)
+      for (size_type i = 0; i < stored.n_elements(); ++i)
+        local_element(i) =
+          internal::get_max(tmp_vector(stored.nth_index_in_set(i)),
+                            local_element(i));
     else
       for (size_type i = 0; i < stored.n_elements(); ++i)
         local_element(i) = tmp_vector(stored.nth_index_in_set(i));
@@ -438,6 +487,42 @@ namespace LinearAlgebra
         for (int i = 0; i < size; ++i)
           values[i] += new_values[i];
       }
+    else if (operation == VectorOperation::min)
+      {
+        const int err = target_vector.Import(multivector, import, Add);
+        AssertThrow(err == 0,
+                    ExcMessage("Epetra Import() failed with error code: " +
+                               Utilities::to_string(err)));
+
+        const double *new_values = target_vector.Values();
+        const int     size       = target_vector.MyLength();
+        Assert(size == 0 || values != nullptr,
+               ExcInternalError("Import failed."));
+
+        // To ensure that this code also compiles with complex
+        // numbers, we only compare the real part of the
+        // variable. Note that min/max do not make sense with complex
+        // numbers.
+        for (int i = 0; i < size; ++i)
+          if (std::real(new_values[i]) - std::real(values[i]) < 0.0)
+            values[i] = new_values[i];
+      }
+    else if (operation == VectorOperation::max)
+      {
+        const int err = target_vector.Import(multivector, import, Add);
+        AssertThrow(err == 0,
+                    ExcMessage("Epetra Import() failed with error code: " +
+                               Utilities::to_string(err)));
+
+        const double *new_values = target_vector.Values();
+        const int     size       = target_vector.MyLength();
+        Assert(size == 0 || values != nullptr,
+               ExcInternalError("Import failed."));
+
+        for (int i = 0; i < size; ++i)
+          if (std::real(new_values[i]) - std::real(values[i]) > 0.0)
+            values[i] = new_values[i];
+      }
     else
       AssertThrow(false, ExcNotImplemented());
   }
@@ -500,7 +585,7 @@ namespace LinearAlgebra
                                             cudaMemcpyDeviceToHost);
         AssertCuda(error_code);
       }
-    else
+    else if (operation == VectorOperation::add)
       {
         // Copy the vector from the device to a temporary vector on the host
         std::vector<Number> tmp(n_elements);
@@ -514,6 +599,40 @@ namespace LinearAlgebra
         for (unsigned int i = 0; i < n_elements; ++i)
           values[i] += tmp[i];
       }
+    else if (operation == VectorOperation::min)
+      {
+        // Copy the vector from the device to a temporary vector on the host
+        std::vector<Number> tmp(n_elements);
+        cudaError_t         error_code = cudaMemcpy(tmp.data(),
+                                            cuda_vec.get_values(),
+                                            n_elements * sizeof(Number),
+                                            cudaMemcpyDeviceToHost);
+        AssertCuda(error_code);
+
+        // To ensure that this code also compiles with complex
+        // numbers, we only compare the real part of the
+        // variable. Note that min/max do not make sense with complex
+        // numbers.
+        for (unsigned int i = 0; i < n_elements; ++i)
+          if (std::real(tmp[i]) - std::real(values[i]) < 0.0)
+            values[i] = tmp[i];
+      }
+    else if (operation == VectorOperation::max)
+      {
+        // Copy the vector from the device to a temporary vector on the host
+        std::vector<Number> tmp(n_elements);
+        cudaError_t         error_code = cudaMemcpy(tmp.data(),
+                                            cuda_vec.get_values(),
+                                            n_elements * sizeof(Number),
+                                            cudaMemcpyDeviceToHost);
+        AssertCuda(error_code);
+
+        for (unsigned int i = 0; i < n_elements; ++i)
+          if (std::real(tmp[i]) - std::real(values[i]) > 0.0)
+            values[i] = tmp[i];
+      }
+    else
+      AssertThrow(false, ExcNotImplemented());
   }
 #endif
 
index 6e31772cfb56ff98c409cfc4b057319e6ced47b5..4b3fc163de11999c0955d3b86ea4809631ebc959 100644 (file)
@@ -50,7 +50,15 @@ struct VectorOperation
     /**
      * The current operation is an addition.
      */
-    add
+    add,
+    /**
+     * The current operation is a minimization.
+     */
+    min,
+    /**
+     * The current operation is a maximization.
+     */
+    max
   };
 };
 
index 7f6b1c40e292e564c8a82c0668316d0ac11d1759..725dc483c3edd66019208f6330186f252f2e426d 100644 (file)
@@ -615,7 +615,7 @@ namespace LinearAlgebra
                                               cudaMemcpyHostToDevice);
           AssertCuda(error_code);
         }
-      else
+      else if (operation == VectorOperation::add)
         {
           // Create a temporary vector on the device
           Number *    tmp;
@@ -644,6 +644,8 @@ namespace LinearAlgebra
           error_code = cudaFree(tmp);
           AssertCuda(error_code);
         }
+      else
+        AssertThrow(false, ExcNotImplemented());
     }
 
 
index 654a3668edcee1107e0ee2466a69ec94afadc7ee..36d38cf164c07904d73f1928b87a849d84d6ad51 100644 (file)
@@ -187,8 +187,14 @@ namespace LinearAlgebra
 
       if (operation == VectorOperation::insert)
         vector->Export(source_vector, import, Insert);
-      else
+      else if (operation == VectorOperation::add)
         vector->Export(source_vector, import, Add);
+      else if (operation == VectorOperation::max)
+        vector->Export(source_vector, import, Epetra_Max);
+      else if (operation == VectorOperation::min)
+        vector->Export(source_vector, import, Epetra_Min);
+      else
+        AssertThrow(false, ExcNotImplemented());
     }
 
 
index fe657772ef6d0225f4072c85cdbef2ae7ffb9856..1eee3a76cd7590dcd91b83db15b37455db65b671 100644 (file)
@@ -1155,6 +1155,10 @@ namespace TrilinosWrappers
           mode = Add;
         else if (operation == ::dealii::VectorOperation::insert)
           mode = Insert;
+        else
+          Assert(
+            false,
+            "compress() can only be called with VectorOperation add, insert, or unknown");
       }
     else
       {
index f2acfcec4bfce13e0ba33761d74c89fe2dfbc07a..f143f26209a25236b0a8d7afb194f52e0fd3c0d1 100644 (file)
@@ -560,6 +560,10 @@ namespace TrilinosWrappers
             mode = Add;
           else if (given_last_action == ::dealii::VectorOperation::insert)
             mode = Insert;
+          else
+            Assert(
+              false,
+              "compress() can only be called with VectorOperation add, insert, or unknown");
         }
       else
         {
diff --git a/tests/mpi/parallel_vector_22.cc b/tests/mpi/parallel_vector_22.cc
new file mode 100644 (file)
index 0000000..ecae110
--- /dev/null
@@ -0,0 +1,114 @@
+// ---------------------------------------------------------------------
+//
+// 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// check LA::Vector::compress(VectorOperation::min/max) from ghosts
+
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+
+void
+test()
+{
+  unsigned int myid    = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  unsigned int numproc = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+
+  if (myid == 0)
+    deallog << "numproc=" << numproc << std::endl;
+
+
+  // each processor owns 2 indices and all
+  // are ghosting element 1 (the second)
+  IndexSet local_owned(numproc * 2);
+  local_owned.add_range(myid * 2, myid * 2 + 2);
+  IndexSet local_relevant(numproc * 2);
+  local_relevant = local_owned;
+  local_relevant.add_range(1, 2);
+
+  // create vector
+  LinearAlgebra::distributed::Vector<double> v(local_owned,
+                                               local_relevant,
+                                               MPI_COMM_WORLD);
+
+  // set local values
+  v(myid * 2)     = myid * 2.0;
+  v(myid * 2 + 1) = myid * 2.0 + 1.0;
+  v.compress(VectorOperation::add);
+  v *= 2.0;
+
+  // check setup of vectors
+  deallog << myid << ":"
+          << "first owned entry: " << v(myid * 2) << std::endl;
+  deallog << myid << ":"
+          << "second owned entry: " << v(myid * 2 + 1) << std::endl;
+
+  // set ghost dof on not owning processor and maximize
+  if (myid)
+    v(1) = 7. * myid;
+  v.compress(VectorOperation::max);
+
+  // import ghosts onto all procs
+  v.update_ghost_values();
+
+  // check
+  deallog << myid << ":"
+          << "ghost entry: " << v(1) << std::endl;
+
+  // ghosts are set to zero
+  v.zero_out_ghosts();
+
+  // minimize
+  v.compress(VectorOperation::min);
+  v.update_ghost_values();
+
+  // check
+  deallog << myid << ":"
+          << "ghost entry: " << v(1) << std::endl;
+
+  // update of ghost value from owner and minimize
+  v.zero_out_ghosts();
+  if (!myid)
+    v(1) = -1.;
+  v.compress(VectorOperation::min);
+  v.update_ghost_values();
+
+  // check
+  deallog << myid << ":"
+          << "ghost entry: " << v(1) << std::endl;
+
+  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());
+
+  MPILogInitAll log;
+
+  test();
+}
diff --git a/tests/mpi/parallel_vector_22.mpirun=10.output b/tests/mpi/parallel_vector_22.mpirun=10.output
new file mode 100644 (file)
index 0000000..8a13c7f
--- /dev/null
@@ -0,0 +1,71 @@
+
+DEAL:0::numproc=10
+DEAL:0::0:first owned entry: 0.00000
+DEAL:0::0:second owned entry: 2.00000
+DEAL:0::0:ghost entry: 63.0000
+DEAL:0::0:ghost entry: 0.00000
+DEAL:0::0:ghost entry: -1.00000
+DEAL:0::OK
+
+DEAL:1::1:first owned entry: 4.00000
+DEAL:1::1:second owned entry: 6.00000
+DEAL:1::1:ghost entry: 63.0000
+DEAL:1::1:ghost entry: 0.00000
+DEAL:1::1:ghost entry: -1.00000
+
+
+DEAL:2::2:first owned entry: 8.00000
+DEAL:2::2:second owned entry: 10.0000
+DEAL:2::2:ghost entry: 63.0000
+DEAL:2::2:ghost entry: 0.00000
+DEAL:2::2:ghost entry: -1.00000
+
+
+DEAL:3::3:first owned entry: 12.0000
+DEAL:3::3:second owned entry: 14.0000
+DEAL:3::3:ghost entry: 63.0000
+DEAL:3::3:ghost entry: 0.00000
+DEAL:3::3:ghost entry: -1.00000
+
+
+DEAL:4::4:first owned entry: 16.0000
+DEAL:4::4:second owned entry: 18.0000
+DEAL:4::4:ghost entry: 63.0000
+DEAL:4::4:ghost entry: 0.00000
+DEAL:4::4:ghost entry: -1.00000
+
+
+DEAL:5::5:first owned entry: 20.0000
+DEAL:5::5:second owned entry: 22.0000
+DEAL:5::5:ghost entry: 63.0000
+DEAL:5::5:ghost entry: 0.00000
+DEAL:5::5:ghost entry: -1.00000
+
+
+DEAL:6::6:first owned entry: 24.0000
+DEAL:6::6:second owned entry: 26.0000
+DEAL:6::6:ghost entry: 63.0000
+DEAL:6::6:ghost entry: 0.00000
+DEAL:6::6:ghost entry: -1.00000
+
+
+DEAL:7::7:first owned entry: 28.0000
+DEAL:7::7:second owned entry: 30.0000
+DEAL:7::7:ghost entry: 63.0000
+DEAL:7::7:ghost entry: 0.00000
+DEAL:7::7:ghost entry: -1.00000
+
+
+DEAL:8::8:first owned entry: 32.0000
+DEAL:8::8:second owned entry: 34.0000
+DEAL:8::8:ghost entry: 63.0000
+DEAL:8::8:ghost entry: 0.00000
+DEAL:8::8:ghost entry: -1.00000
+
+
+DEAL:9::9:first owned entry: 36.0000
+DEAL:9::9:second owned entry: 38.0000
+DEAL:9::9:ghost entry: 63.0000
+DEAL:9::9:ghost entry: 0.00000
+DEAL:9::9:ghost entry: -1.00000
+
diff --git a/tests/mpi/parallel_vector_22.mpirun=4.output b/tests/mpi/parallel_vector_22.mpirun=4.output
new file mode 100644 (file)
index 0000000..4c36987
--- /dev/null
@@ -0,0 +1,29 @@
+
+DEAL:0::numproc=4
+DEAL:0::0:first owned entry: 0.00000
+DEAL:0::0:second owned entry: 2.00000
+DEAL:0::0:ghost entry: 21.0000
+DEAL:0::0:ghost entry: 0.00000
+DEAL:0::0:ghost entry: -1.00000
+DEAL:0::OK
+
+DEAL:1::1:first owned entry: 4.00000
+DEAL:1::1:second owned entry: 6.00000
+DEAL:1::1:ghost entry: 21.0000
+DEAL:1::1:ghost entry: 0.00000
+DEAL:1::1:ghost entry: -1.00000
+
+
+DEAL:2::2:first owned entry: 8.00000
+DEAL:2::2:second owned entry: 10.0000
+DEAL:2::2:ghost entry: 21.0000
+DEAL:2::2:ghost entry: 0.00000
+DEAL:2::2:ghost entry: -1.00000
+
+
+DEAL:3::3:first owned entry: 12.0000
+DEAL:3::3:second owned entry: 14.0000
+DEAL:3::3:ghost entry: 21.0000
+DEAL:3::3:ghost entry: 0.00000
+DEAL:3::3:ghost entry: -1.00000
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.