Merge branch 'deterministic-indicesSyncer' into 'master'

Add flag useFixedOrder to the coarsen method of AMGs ParallelIndicesCoarsener.

See merge request !553

CHANGELOG.md

@@ -15,6 +15,14 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
   performance. Hence, `MatrixIndexSet` can no longer be used for very large matrices with more
   than 2^32 columns.
 
+- Added flag 'useFixedOrder' to the coarsen method of AMGs ParallelIndicesCoarsener.
+  If set to true, during the creation of the coarser matrix (by accumulation and restriction
+  to fewer processes) the coarse matrix rows are ordered in a fixed manner, making parallel runs
+  reproducible; the runtime is possibly not ideal though.
+  If set to false (which is the default), the row order depends on the order of messages received from the
+  processes responsible for the respective parts of the finer grid. Then the indices on the coarser grid
+  may differ from run to run.
+
 ## Deprecations and removals
 
 - The deprecated CMake variables `SUPERLU_INCLUDE_DIRS`, `SUPERLU_LIBRARIES`,

dune/istl/paamg/indicescoarsener.hh

@@ -83,6 +83,10 @@ namespace Dune
        * @param aggregates The mapping of unknowns onto aggregates.
        * @param coarseInfo The information about the parallel data decomposition
        * on the coarse level.
+       * @param noAggregates
+       * @param useFixedOrder Flag indicating if creating indices for the coarser level
+       * should be done in a fixed order, i.e., the order in which the rows were sent.
+       * If set to true, this makes the runs reproducible but it might slow down performance.
        * @return The number of unknowns on the coarse level.
        */
       template<typename Graph, typename VM>
@@ -92,7 +96,8 @@ namespace Dune
               VM& visitedMap,
               AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
               ParallelInformation& coarseInfo,
-              typename Graph::VertexDescriptor noAggregates);
+              typename Graph::VertexDescriptor noAggregates,
+              bool useFixedOrder = false);
 
     private:
       template<typename G, typename I>
@@ -187,7 +192,8 @@ namespace Dune
                                            const AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
                                            ParallelIndexSet& coarseIndices,
                                            RemoteIndices& coarseRemote,
-                                           ParallelAggregateRenumberer<Graph,I>& renumberer);
+                                           ParallelAggregateRenumberer<Graph,I>& renumberer,
+                                           bool useFixedOrder);
 
     };
 
@@ -219,7 +225,8 @@ namespace Dune
               VM& visitedMap,
               AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
               SequentialInformation& coarseInfo,
-              typename Graph::VertexDescriptor noAggregates);
+              typename Graph::VertexDescriptor noAggregates,
+              bool useFixedOrder = false);
     };
 
 #if HAVE_MPI
@@ -231,13 +238,14 @@ namespace Dune
                                            VM& visitedMap,
                                            AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
                                            ParallelInformation& coarseInfo,
-                                           [[maybe_unused]] typename Graph::VertexDescriptor noAggregates)
+                                           [[maybe_unused]] typename Graph::VertexDescriptor noAggregates,
+                                           bool useFixedOrder)
     {
       ParallelAggregateRenumberer<Graph,typename ParallelInformation::GlobalLookupIndexSet> renumberer(aggregates, fineInfo.globalLookup());
       buildCoarseIndexSet(fineInfo, fineGraph, visitedMap, aggregates,
                           coarseInfo.indexSet(), renumberer);
       buildCoarseRemoteIndices(fineInfo.remoteIndices(), aggregates, coarseInfo.indexSet(),
-                               coarseInfo.remoteIndices(), renumberer);
+                               coarseInfo.remoteIndices(), renumberer, useFixedOrder);
 
       return renumberer;
     }
@@ -318,7 +326,8 @@ namespace Dune
                                                                  const AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
                                                                  ParallelIndexSet& coarseIndices,
                                                                  RemoteIndices& coarseRemote,
-                                                                 ParallelAggregateRenumberer<Graph,I>& renumberer)
+                                                                 ParallelAggregateRenumberer<Graph,I>& renumberer,
+                                                                 bool useFixedOrder)
     {
       std::vector<char> attributes(static_cast<std::size_t>(renumberer));
 
@@ -375,7 +384,7 @@ namespace Dune
       // sync the index set and the remote indices to recompute missing
       // indices
       IndicesSyncer<ParallelIndexSet> syncer(coarseIndices, coarseRemote);
-      syncer.sync(renumberer);
+      syncer.sync(renumberer, useFixedOrder);
 
     }
 
@@ -390,7 +399,8 @@ namespace Dune
       [[maybe_unused]] VM& visitedMap,
       [[maybe_unused]] AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
       [[maybe_unused]] SequentialInformation& coarseInfo,
-      [[maybe_unused]] typename Graph::VertexDescriptor noAggregates)
+      [[maybe_unused]] typename Graph::VertexDescriptor noAggregates,
+      [[maybe_unused]] bool useFixedOrder)
     {
       return noAggregates;
     }

dune/istl/paamg/matrixhierarchy.hh

@@ -297,10 +297,13 @@ namespace Dune
        * coarsening will stop (default: 1.2)
        * @param prolongDamp The damping factor to apply to the prolongated update (default: 1.6)
        * @param accumulate Whether to accumulate the data onto fewer processors on coarser levels.
+       * @param useFixedOrder Flag indicating if creating indices for the coarser level
+       * should be done in a fixed order, i.e., the order in which the rows were sent.
+       * If set to true, this makes the runs reproducible but it might slow down performance.
        */
       CoarsenCriterion(int maxLevel=100, int coarsenTarget=1000, double minCoarsenRate=1.2,
-                       double prolongDamp=1.6, AccumulationMode accumulate=successiveAccu)
-        : AggregationCriterion(Dune::Amg::Parameters(maxLevel, coarsenTarget, minCoarsenRate, prolongDamp, accumulate))
+                       double prolongDamp=1.6, AccumulationMode accumulate=successiveAccu, bool useFixedOrder = false)
+        : AggregationCriterion(Dune::Amg::Parameters(maxLevel, coarsenTarget, minCoarsenRate, prolongDamp, accumulate, useFixedOrder))
       {}
 
       CoarsenCriterion(const Dune::Amg::Parameters& parms)
@@ -610,7 +613,8 @@ namespace Dune
                                    visitedMap,
                                    *aggregatesMap,
                                    *infoLevel,
-                                   noAggregates);
+                                   noAggregates,
+                                   criterion.useFixedOrder());
         GraphCreator::free(graphs);
 
         if(criterion.debugLevel()>2) {

dune/istl/paamg/parameters.hh

@@ -314,7 +314,7 @@ namespace Dune
         return accumulate_;
       }
       /**
-       * @brief Set whether he data should be accumulated on fewer processes on coarser levels.
+       * @brief Set whether the data should be accumulated on fewer processes on coarser levels.
        */
       void setAccumulate(AccumulationMode accu)
       {
@@ -324,6 +324,20 @@ namespace Dune
       void setAccumulate(bool accu){
         accumulate_=accu ? successiveAccu : noAccu;
       }
+
+      /**
+       * @brief Check if the indices for the coarser levels should be created in a fixed order.
+       */
+      bool useFixedOrder() const
+      {
+        return useFixedOrder_;
+      }
+
+      void setUseFixedOrder(bool useFixedOrder)
+      {
+        useFixedOrder_ = useFixedOrder;
+      }
+
       /**
        * @brief Set the damping factor for the prolongation.
        *
@@ -352,11 +366,15 @@ namespace Dune
        * coarsening will stop (default: 1.2)
        * @param prolongDamp The damping factor to apply to the prolongated update (default: 1.6)
        * @param accumulate Whether to accumulate the data onto fewer processors on coarser levels.
+       * @param useFixedOrder Flag indicating if creating indices for the coarser level
+       * should be done in a fixed order, i.e., the order in which the rows were sent.
+       * If set to true, this makes the runs reproducible but it might slow down performance.
        */
       CoarseningParameters(int maxLevel=100, int coarsenTarget=1000, double minCoarsenRate=1.2,
-                           double prolongDamp=1.6, AccumulationMode accumulate=successiveAccu)
+                           double prolongDamp=1.6, AccumulationMode accumulate=successiveAccu,
+                           bool useFixedOrder = false)
         : maxLevel_(maxLevel), coarsenTarget_(coarsenTarget), minCoarsenRate_(minCoarsenRate),
-          dampingFactor_(prolongDamp), accumulate_( accumulate)
+          dampingFactor_(prolongDamp), accumulate_( accumulate), useFixedOrder_(useFixedOrder)
       {}
 
     private:
@@ -381,6 +399,12 @@ namespace Dune
        * coarser levels.
        */
       AccumulationMode accumulate_;
+      /**
+       * @brief useFixedOrder Flag indicating if creating indices for the coarser
+       * level should be done in a fixed order. If set to true, this makes the runs reproducible
+       * but it might slow down performance.
+       */
+      bool useFixedOrder_;
     };
 
     /**
@@ -491,8 +515,8 @@ namespace Dune
        * @param accumulate Whether to accumulate the data onto fewer processors on coarser levels.
        */
       Parameters(int maxLevel=100, int coarsenTarget=1000, double minCoarsenRate=1.2,
-                 double prolongDamp=1.6, AccumulationMode accumulate=successiveAccu)
-        : CoarseningParameters(maxLevel, coarsenTarget, minCoarsenRate, prolongDamp, accumulate)
+                 double prolongDamp=1.6, AccumulationMode accumulate=successiveAccu, bool useFixedOrder = false)
+        : CoarseningParameters(maxLevel, coarsenTarget, minCoarsenRate, prolongDamp, accumulate, useFixedOrder)
           , debugLevel_(2), preSmoothSteps_(2), postSmoothSteps_(2), gamma_(1),
           additive_(false)
       {}

dune/istl/paamg/test/CMakeLists.txt

@@ -92,6 +92,11 @@ dune_add_test(SOURCES hierarchytest.cc
               TIMEOUT 600
               CMAKE_GUARD MPI_FOUND)
 
+dune_add_test(SOURCES coarsentest.cc
+              MPI_RANKS 1 2 4
+              TIMEOUT 600
+              CMAKE_GUARD MPI_FOUND)
+
 dune_add_test(NAME pamgtest
               SOURCES parallelamgtest.cc
               MPI_RANKS 1 2 4

dune/istl/paamg/test/coarsentest.cc

+// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
+// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include <config.h>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/istl/paamg/matrixhierarchy.hh>
+#include <dune/istl/owneroverlapcopy.hh>
+#include <dune/istl/schwarz.hh>
+#include "anisotropic.hh"
+
+constexpr int blockSize = 10;
+
+using Communication = Dune::OwnerOverlapCopyCommunication<int,int>;
+using ParallelIndexSet = Communication::ParallelIndexSet;
+using MatrixBlock = Dune::FieldMatrix<double,blockSize,blockSize>;
+using BCRSMat = Dune::BCRSMatrix<MatrixBlock>;
+using VectorBlock = Dune::FieldVector<double,blockSize>;
+using Vector = Dune::BlockVector<VectorBlock>;
+
+/**
+ * Check if the matrices A and B are the same.
+ *
+ * @param BCRSMat A
+ * @param BCRSMat B
+ * @return bool - True if A and B are the same, false if they are not.
+ */
+bool areSame(BCRSMat A, BCRSMat B)
+{
+  for (auto rowIt = A.begin(); rowIt != A.end(); ++rowIt)
+    for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
+      for (size_t k = 0; k < blockSize; k++)
+        for (size_t l = 0; l < blockSize; l++) {
+          auto entryA = (*colIt)[k][l];
+          try {
+            auto entryB = B[rowIt.index()][colIt.index()][k][l];
+          } catch (Dune::ISTLError e) {
+            return false; // If the entry B[rowIt.index()][colIt.index()][k][l] does not exist, then the matrices are not the same.
+          }
+          if (std::abs((*colIt)[k][l]) > 1e-7 and std::abs(B[rowIt.index()][colIt.index()][k][l]) > 1e-7
+              and std::abs((*colIt)[k][l] - B[rowIt.index()][colIt.index()][k][l])/(*colIt)[k][l] > 1e-7) {
+            return false; // Check if the entries are approximately the same.
+          }
+        }
+  return true;
+}
+
+/**
+ * Builds two AMG matrix hierarchies and checks if they are equal.
+ * When setting up such matrix hierarchies on multiple processes, they get created in a non-deterministic way and the two matrix hierarchies might differ.
+ *
+ * @param int N Total number of unknowns
+ * @param int rank Rank of current process
+ * @param int noProcesses Number of all processes
+ * @return int - 0 if the hierarchies are the same, 1 if they are not.
+ */
+int testHierarchy(int N, int rank, int noProcesses)
+{
+  int n;
+  Communication pinfo(Dune::MPIHelper::getCommunicator());
+  ParallelIndexSet& indices = pinfo.indexSet();
+
+  using RemoteIndices = Dune::RemoteIndices<ParallelIndexSet>;
+  RemoteIndices& remoteIndices = pinfo.remoteIndices();
+
+  BCRSMat mat = setupAnisotropic2d<MatrixBlock>(N, indices, Dune::MPIHelper::getCommunication(), &n);
+
+  remoteIndices.rebuild<false>();
+
+  using OverlapFlags = Dune::NegateSet<Communication::OwnerSet>;
+  using Operator = Dune::OverlappingSchwarzOperator<BCRSMat,Vector,Vector,Communication>;
+  using Hierarchy = Dune::Amg::MatrixHierarchy<Operator,Communication>;
+
+  Operator op(mat, pinfo);
+  Hierarchy hierarchyA(
+    Dune::stackobject_to_shared_ptr(op),
+    Dune::stackobject_to_shared_ptr(pinfo));
+
+  Hierarchy hierarchyB(
+    Dune::stackobject_to_shared_ptr(op),
+    Dune::stackobject_to_shared_ptr(pinfo));
+
+  using Criterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >;
+
+  Criterion criterionA(1000,10,1.2,1.6,Dune::Amg::AccumulationMode::successiveAccu,true);
+  Criterion criterionB(1000,10,1.2,1.6,Dune::Amg::AccumulationMode::successiveAccu,true);
+
+  hierarchyA.template build<OverlapFlags>(criterionA);
+  hierarchyB.template build<OverlapFlags>(criterionB);
+
+  MPI_Barrier(MPI_COMM_WORLD);
+  //Check if the second finest matrices are the same
+  auto finestMatrixA = hierarchyA.matrices().finest();
+  finestMatrixA++;
+  auto finestMatrixB = hierarchyB.matrices().finest();
+  finestMatrixB++;
+  if (!areSame(finestMatrixA->getmat(), finestMatrixB->getmat())) {
+      std::cerr << "During run with " << noProcesses << " processes: The matrix hierarchy on rank " << rank << " was created in a non-deterministic way." << std::endl;
+      return 1;
+  }
+  return 0;
+}
+
+
+int main(int argc, char** argv)
+{
+  Dune::MPIHelper& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+  int N=10;
+
+  if(argc>1)
+    N = atoi(argv[1]);
+
+  return testHierarchy(N, mpiHelper.rank(), mpiHelper.size());
+}