single pass algorithm for both Image and Graph connectivity

vtkm/worklet/connectivities/GraphConnectivity.h

@@ -36,6 +36,69 @@ public:
 
   using InputDomain = _1;
 
+  template <typename Parents>
+  VTKM_EXEC vtkm::Id findRoot(const Parents& parents, vtkm::Id index) const
+  {
+    while (parents.Get(index) != index)
+      index = parents.Get(index);
+    return index;
+  }
+
+  template <typename Parents>
+  VTKM_EXEC void Unite(Parents& parents, vtkm::Id u, vtkm::Id v) const
+  {
+    // Data Race Resolutions
+    // Since this function modifies the Union-Find data structure, concurrent
+    // invocation of it by 2 or more threads cause potential data race. Here
+    // is a case analysis why the potential data race does no harm in the
+    // context of the iterative connected component algorithm.
+
+    // Case 1, Two threads calling Unite(u, v) concurrently.
+    // Problem: One thread might attach u to v while the other thread attach
+    // v to u, causing a cycle in the Union-Find data structure.
+    // Resolution: This is not necessary a data race issue. This is resolved by
+    // "linking by index" as in SV Jayanti et.al. with less than as the total order.
+    // The two threads will make the same decision on how to Unite the two tree
+    // (e.g. from root with larger id to root with smaller id.) This avoids cycles in
+    // the resulting graph and maintains the rooted forest structure of Union-Find.
+
+    // Case 2, T0 calling Unite(u, v), T1 calling Unite(u, w) and T2 calling
+    // Unite(v, s) concurrently.
+    // Problem I: There is a potential write after read data race. After T0
+    // calls findRoot for u and v, T1 might have called parents.Set(root_u, root_w)
+    // thus changed root_u to root_w thus making root_u "obsolete" before T0 calls
+    // parents.Set() on root_u/root_v.
+    // When the root of the tree to be attached to (e.g. root_u, when root_u <  root_v)
+    // is changed, there is no hazard, since we are just attaching a tree to a
+    // now a non-root node, root_u, (thus, root_w <- root_u <- root_v) and three
+    // components merged.
+    // However, when the root of the attaching tree (root_v) is change, it
+    // means that the root_u has been attached to yet some other root_s and became
+    // a non-root node. If we are now attaching this non-root node to root_w we
+    // would leave root_s behind and undoing previous work.
+    // Resolution:
+    auto root_u = findRoot(parents, u);
+    auto root_v = findRoot(parents, v);
+
+    // There is a potential concurrent write data race as it is possible for
+    // two threads to try to change the same old root to different new roots,
+    // e.g. threadA calls parents.Set(root, rootB) while threadB calls
+    // parents(root, rootB) where rootB < root and rootC < root (but the order
+    // of rootA and rootB is unspecified.) Each thread assumes success while
+    // the outcome is actually unspecified. An atomic Compare and Swap is
+    // suggested in SV Janati et. al. to "resolve" data race. However, I don't
+    // see any need to use CAS, it looks like the data race will always correct
+    // itself by the algorithm in later iterations as long as atomic Store of
+    // memory_order_release and Load of memory_order_acquire is used (as provided
+    // by AtomicArrayInOut.) This memory consistency model is the default mode
+    // for x86, thus having zero extra cost but might be required for CUDA and/or ARM.
+    if (root_u < root_v)
+      parents.Set(root_v, root_u);
+    else if (root_u > root_v)
+      parents.Set(root_u, root_v);
+    // else, no need to do anything when they are the same set.
+  }
+
   // TODO: Use Scatter?
   template <typename InPortalType, typename InOutPortalType>
   VTKM_EXEC void operator()(vtkm::Id index,
@@ -47,13 +110,9 @@ public:
     for (vtkm::Id offset = start; offset < start + degree; offset++)
     {
       vtkm::Id neighbor = conn.Get(offset);
-      // Note: comp.Get(index) == comp.Get(comp.Get(index)) applies for both the
-      // root of the tree and the first level vertices. To check for root, use
-      // index == comp.Get(index).
-      if ((comp.Get(index) == comp.Get(comp.Get(index))) && (comp.Get(neighbor) < comp.Get(index)))
-      {
-        comp.Set(comp.Get(index), comp.Get(neighbor));
-      }
+      auto thisComp = comp.Get(index);
+      auto thatComp = comp.Get(neighbor);
+      Unite(comp, thisComp, thatComp);
     }
   }
 };
@@ -70,29 +129,15 @@ public:
            const InputPortalType& connectivityArray,
            OutputPortalType& componentsOut) const
   {
-    bool everythingIsAStar = false;
     vtkm::cont::ArrayHandle<vtkm::Id> components;
     Algorithm::Copy(
       vtkm::cont::ArrayHandleCounting<vtkm::Id>(0, 1, numIndicesArray.GetNumberOfValues()),
       components);
 
-    //used as an atomic bool, so we use Int32 as it the
-    //smallest type that VTK-m supports as atomics
-    vtkm::cont::ArrayHandle<vtkm::Int32> allStars;
-    allStars.Allocate(1);
-
     vtkm::cont::Invoker invoke;
 
-    do
-    {
-      allStars.WritePortal().Set(0, 1); //reset the atomic state
-      invoke(detail::Graft{}, indexOffsetsArray, numIndicesArray, connectivityArray, components);
-
-      // Detection of allStars has to come before pointer jumping. Don't try to rearrange it.
-      invoke(IsStar{}, components, allStars);
-      everythingIsAStar = (allStars.WritePortal().Get(0) == 1);
-      invoke(PointerJumping{}, components);
-    } while (!everythingIsAStar);
+    invoke(detail::Graft{}, indexOffsetsArray, numIndicesArray, connectivityArray, components);
+    invoke(PointerJumping{}, components);
 
     // renumber connected component to the range of [0, number of components).
     vtkm::cont::ArrayHandle<vtkm::Id> uniqueComponents;

vtkm/worklet/connectivities/ImageConnectivity.h

@@ -34,10 +34,9 @@ public:
   using ControlSignature = void(CellSetIn,
                                 FieldInNeighborhood compIn,
                                 FieldInNeighborhood color,
-                                AtomicArrayInOut compOut,
-                                AtomicArrayInOut changed);
+                                AtomicArrayInOut compOut);
 
-  using ExecutionSignature = void(Boundary, _2, _3, _4, _5);
+  using ExecutionSignature = void(Boundary, _2, _3, _4);
 
   // This is the naive find() without path compaction in SV Jayanti et. al.
   // Since the parents array is read-only there is no data race.

vtkm/worklet/connectivities/UnionFind.h

@@ -68,29 +68,6 @@ public:
   }
 };
 
-class IsStar : public vtkm::worklet::WorkletMapField
-{
-public:
-  using ControlSignature = void(WholeArrayIn comp, AtomicArrayInOut);
-  using ExecutionSignature = void(WorkIndex, _1, _2);
-  using InputDomain = _1;
-
-  template <typename InPortalType, typename AtomicInOut>
-  VTKM_EXEC void operator()(vtkm::Id index, const InPortalType& comp, AtomicInOut& hasStar) const
-  {
-    //hasStar emulates a LogicalAnd across all the values
-    //where we start with a value of 'true'|1.
-    // Note: comp.Get(index) == comp.Get(comp.Get(index)) applies for both the
-    // root of the tree and the first level vertices. If all vertices
-    // is either a root or first level vertices, it is a rooted star.
-    const bool isAStar = (comp.Get(index) == comp.Get(comp.Get(index)));
-    if (!isAStar && hasStar.Get(0) == 1)
-    {
-      hasStar.Set(0, 0);
-    }
-  }
-};
-
 } // connectivity
 } // worklet
 } // vtkm