Add ScanExtended device algorithm.

This behaves just like `ScanExclusive`, but rather than returning the
total sum, it is appended to the end of the output array.

This is in preparation for the CellSetExplicit refactoring described in
issue #408.

docs/changelog/scan-extended.md

@@ -0,0 +1,29 @@
+# A `ScanExtended` device algorithm has been added.
+
+This new scan algorithm produces an array that contains both an inclusive scan
+and an exclusive scan in the same array:
+
+```
+#include <vtkm/cont/Algorithm.h>
+#include <vtkm/cont/ArrayGetValue.h>
+#include <vtkm/cont/ArrayHandle.h>
+#include <vtkm/cont/ArrayHandleView.h>
+
+vtkm::cont::ArrayHandle<T> inputData = ...;
+const vtkm::Id size = inputData.GetNumberOfValues();
+
+vtkm::cont::ArrayHandle<T> extendedScan;
+vtkm::cont::Algorithm::ScanExtended(inputData, extendedScan);
+
+// The exclusive scan is the first `inputSize` values starting at index 0:
+auto exclusiveScan = vtkm::cont::make_ArrayHandleView(extendedScan, 0, size);
+
+// The inclusive scan is the first `inputSize` values starting at index 1:
+auto inclusiveScan = vtkm::cont::make_ArrayHandleView(extendedScan, 1, size);
+
+// The total sum of the input data is the last value in the extended scan.
+const T totalSum = vtkm::cont::ArrayGetValue(size, extendedScan);
+```
+
+This can also be thought of as an exclusive scan that appends the total sum,
+rather than returning it.

vtkm/cont/Algorithm.h

@@ -317,6 +317,19 @@ struct ScanExclusiveByKeyFunctor
   }
 };
 
+template <typename T>
+struct ScanExtendedFunctor
+{
+  template <typename Device, typename... Args>
+  VTKM_CONT bool operator()(Device, Args&&... args)
+  {
+    VTKM_IS_DEVICE_ADAPTER_TAG(Device);
+    vtkm::cont::DeviceAdapterAlgorithm<Device>::ScanExtended(
+      PrepareArgForExec<Device>(std::forward<Args>(args))...);
+    return true;
+  }
+};
+
 struct ScheduleFunctor
 {
   template <typename Device, typename... Args>
@@ -976,6 +989,42 @@ struct Algorithm
   }
 
 
+  template <typename T, class CIn, class COut>
+  VTKM_CONT static void ScanExtended(vtkm::cont::DeviceAdapterId devId,
+                                     const vtkm::cont::ArrayHandle<T, CIn>& input,
+                                     vtkm::cont::ArrayHandle<T, COut>& output)
+  {
+    detail::ScanExtendedFunctor<T> functor;
+    vtkm::cont::TryExecuteOnDevice(devId, functor, input, output);
+  }
+  template <typename T, class CIn, class COut>
+  VTKM_CONT static void ScanExtended(const vtkm::cont::ArrayHandle<T, CIn>& input,
+                                     vtkm::cont::ArrayHandle<T, COut>& output)
+  {
+    ScanExtended(vtkm::cont::DeviceAdapterTagAny(), input, output);
+  }
+
+
+  template <typename T, class CIn, class COut, class BinaryFunctor>
+  VTKM_CONT static void ScanExtended(vtkm::cont::DeviceAdapterId devId,
+                                     const vtkm::cont::ArrayHandle<T, CIn>& input,
+                                     vtkm::cont::ArrayHandle<T, COut>& output,
+                                     BinaryFunctor binaryFunctor,
+                                     const T& initialValue)
+  {
+    detail::ScanExtendedFunctor<T> functor;
+    vtkm::cont::TryExecuteOnDevice(devId, functor, input, output, binaryFunctor, initialValue);
+  }
+  template <typename T, class CIn, class COut, class BinaryFunctor>
+  VTKM_CONT static void ScanExtended(const vtkm::cont::ArrayHandle<T, CIn>& input,
+                                     vtkm::cont::ArrayHandle<T, COut>& output,
+                                     BinaryFunctor binaryFunctor,
+                                     const T& initialValue)
+  {
+    ScanExtended(vtkm::cont::DeviceAdapterTagAny(), input, output, binaryFunctor, initialValue);
+  }
+
+
   template <class Functor>
   VTKM_CONT static void Schedule(vtkm::cont::DeviceAdapterId devId,
                                  Functor functor,

vtkm/cont/DeviceAdapterAlgorithm.h

@@ -378,6 +378,55 @@ struct DeviceAdapterAlgorithm
                                            const vtkm::cont::ArrayHandle<U, VIn>& values,
                                            vtkm::cont::ArrayHandle<U, VOut>& output);
 
+  /// \brief Compute an extended prefix sum operation on the input ArrayHandle.
+  ///
+  /// Computes an extended prefix sum operation on the \c input ArrayHandle,
+  /// storing the results in the \c output ArrayHandle. This produces an output
+  /// array that contains both an inclusive scan (in elements [1, size)) and an
+  /// exclusive scan (in elements [0, size-1)). By using ArrayHandleView,
+  /// arrays containing both inclusive and exclusive scans can be generated
+  /// from an extended scan with minimal memory usage.
+  ///
+  /// This algorithm may also be more efficient than ScanInclusive and
+  /// ScanExclusive on some devices, since it may be able to avoid copying the
+  /// total sum to the control environment to return.
+  ///
+  /// ScanExtended is similar to the stl partial sum function, exception that
+  /// ScanExtended doesn't do a serial summation. This means that if you have
+  /// defined a custom plus operator for T it must be associative, or you will
+  /// get inconsistent results.
+  ///
+  /// This overload of ScanExtended uses vtkm::Add for the binary functor, and
+  /// uses zero for the initial value of the scan operation.
+  ///
+  template <typename T, class CIn, class COut>
+  VTKM_CONT static void ScanExtended(const vtkm::cont::ArrayHandle<T, CIn>& input,
+                                     vtkm::cont::ArrayHandle<T, COut>& output);
+
+  /// \brief Compute an extended prefix sum operation on the input ArrayHandle.
+  ///
+  /// Computes an extended prefix sum operation on the \c input ArrayHandle,
+  /// storing the results in the \c output ArrayHandle. This produces an output
+  /// array that contains both an inclusive scan (in elements [1, size)) and an
+  /// exclusive scan (in elements [0, size-1)). By using ArrayHandleView,
+  /// arrays containing both inclusive and exclusive scans can be generated
+  /// from an extended scan with minimal memory usage.
+  ///
+  /// This algorithm may also be more efficient than ScanInclusive and
+  /// ScanExclusive on some devices, since it may be able to avoid copying the
+  /// total sum to the control environment to return.
+  ///
+  /// ScanExtended is similar to the stl partial sum function, exception that
+  /// ScanExtended doesn't do a serial summation. This means that if you have
+  /// defined a custom plus operator for T it must be associative, or you will
+  /// get inconsistent results.
+  ///
+  template <typename T, class CIn, class COut, class BinaryFunctor>
+  VTKM_CONT static void ScanExtended(const vtkm::cont::ArrayHandle<T, CIn>& input,
+                                     vtkm::cont::ArrayHandle<T, COut>& output,
+                                     BinaryFunctor binaryFunctor,
+                                     const T& initialValue);
+
   /// \brief Schedule many instances of a function to run on concurrent threads.
   ///
   /// Calls the \c functor on several threads. This is the function used in the

vtkm/cont/internal/DeviceAdapterAlgorithmGeneral.h

@@ -16,6 +16,7 @@
 #include <vtkm/cont/ArrayHandleImplicit.h>
 #include <vtkm/cont/ArrayHandleIndex.h>
 #include <vtkm/cont/ArrayHandleStreaming.h>
+#include <vtkm/cont/ArrayHandleView.h>
 #include <vtkm/cont/ArrayHandleZip.h>
 #include <vtkm/cont/BitField.h>
 #include <vtkm/cont/Logging.h>
@@ -660,6 +661,7 @@ public:
     vtkm::Id numValues = input.GetNumberOfValues();
     if (numValues <= 0)
     {
+      output.Shrink(0);
       return initialValue;
     }
 
@@ -687,6 +689,50 @@ public:
       input, output, vtkm::Sum(), vtkm::TypeTraits<T>::ZeroInitialization());
   }
 
+  //--------------------------------------------------------------------------
+  // Scan Exclusive Extend
+  template <typename T, class CIn, class COut, class BinaryFunctor>
+  VTKM_CONT static void ScanExtended(const vtkm::cont::ArrayHandle<T, CIn>& input,
+                                     vtkm::cont::ArrayHandle<T, COut>& output,
+                                     BinaryFunctor binaryFunctor,
+                                     const T& initialValue)
+  {
+    VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf);
+
+    vtkm::Id numValues = input.GetNumberOfValues();
+    if (numValues <= 0)
+    {
+      output.Allocate(1);
+      output.GetPortalControl().Set(0, initialValue);
+      return;
+    }
+
+    vtkm::cont::ArrayHandle<T, vtkm::cont::StorageTagBasic> inclusiveScan;
+    T result = DerivedAlgorithm::ScanInclusive(input, inclusiveScan, binaryFunctor);
+
+    auto inputPortal = inclusiveScan.PrepareForInput(DeviceAdapterTag());
+    auto outputPortal = output.PrepareForOutput(numValues + 1, DeviceAdapterTag());
+
+    InclusiveToExtendedKernel<decltype(inputPortal), decltype(outputPortal), BinaryFunctor>
+      inclusiveToExtended(inputPortal,
+                          outputPortal,
+                          binaryFunctor,
+                          initialValue,
+                          binaryFunctor(initialValue, result));
+
+    DerivedAlgorithm::Schedule(inclusiveToExtended, numValues + 1);
+  }
+
+  template <typename T, class CIn, class COut>
+  VTKM_CONT static void ScanExtended(const vtkm::cont::ArrayHandle<T, CIn>& input,
+                                     vtkm::cont::ArrayHandle<T, COut>& output)
+  {
+    VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf);
+
+    DerivedAlgorithm::ScanExtended(
+      input, output, vtkm::Sum(), vtkm::TypeTraits<T>::ZeroInitialization());
+  }
+
   //--------------------------------------------------------------------------
   // Scan Exclusive By Key
   template <typename KeyT,

vtkm/cont/internal/FunctorsGeneral.h

@@ -1106,9 +1106,50 @@ struct InclusiveToExclusiveKernel : vtkm::exec::FunctorBase
   VTKM_EXEC
   void operator()(vtkm::Id index) const
   {
-    ValueType result = (index == 0)
+    const ValueType result = (index == 0)
       ? this->InitialValue
       : this->BinaryOperator(this->InitialValue, this->InPortal.Get(index - 1));
+
+    this->OutPortal.Set(index, result);
+  }
+};
+
+template <typename InPortalType, typename OutPortalType, typename BinaryFunctor>
+struct InclusiveToExtendedKernel : vtkm::exec::FunctorBase
+{
+  using ValueType = typename InPortalType::ValueType;
+
+  InPortalType InPortal;
+  OutPortalType OutPortal;
+  BinaryFunctor BinaryOperator;
+  ValueType InitialValue;
+  ValueType FinalValue;
+
+  VTKM_CONT
+  InclusiveToExtendedKernel(const InPortalType& inPortal,
+                            const OutPortalType& outPortal,
+                            BinaryFunctor& binaryOperator,
+                            ValueType initialValue,
+                            ValueType finalValue)
+    : InPortal(inPortal)
+    , OutPortal(outPortal)
+    , BinaryOperator(binaryOperator)
+    , InitialValue(initialValue)
+    , FinalValue(finalValue)
+  {
+  }
+
+  VTKM_SUPPRESS_EXEC_WARNINGS
+  VTKM_EXEC
+  void operator()(vtkm::Id index) const
+  {
+    // The output array has one more value than the input, which holds the
+    // total sum.
+    const ValueType result =
+      (index == 0) ? this->InitialValue : (index == this->InPortal.GetNumberOfValues())
+        ? this->FinalValue
+        : this->BinaryOperator(this->InitialValue, this->InPortal.Get(index - 1));
+
     this->OutPortal.Set(index, result);
   }
 };

vtkm/cont/testing/TestingDeviceAdapter.h

@@ -14,6 +14,7 @@
 #include <vtkm/BinaryPredicates.h>
 #include <vtkm/TypeTraits.h>
 
+#include <vtkm/cont/ArrayGetValues.h>
 #include <vtkm/cont/ArrayHandle.h>
 #include <vtkm/cont/ArrayHandleConstant.h>
 #include <vtkm/cont/ArrayHandleIndex.h>
@@ -1985,6 +1986,108 @@ private:
     }
   }
 
+  static VTKM_CONT void TestScanExtended()
+  {
+    std::cout << "-------------------------------------------" << std::endl;
+    std::cout << "Testing Extended Scan" << std::endl;
+
+    {
+      std::cout << "  size " << ARRAY_SIZE << std::endl;
+
+      //construct the index array
+      IdArrayHandle array;
+      Algorithm::Schedule(ClearArrayKernel(array.PrepareForOutput(ARRAY_SIZE, DeviceAdapterTag())),
+                          ARRAY_SIZE);
+
+      // we now have an array whose sum = (OFFSET * ARRAY_SIZE),
+      // let's validate that
+      Algorithm::ScanExtended(array, array);
+      VTKM_TEST_ASSERT(array.GetNumberOfValues() == ARRAY_SIZE + 1, "Output size incorrect.");
+      auto portal = array.GetPortalConstControl();
+      for (vtkm::Id i = 0; i < ARRAY_SIZE + 1; ++i)
+      {
+        const vtkm::Id value = portal.Get(i);
+        VTKM_TEST_ASSERT(value == i * OFFSET, "Incorrect partial sum");
+      }
+
+      std::cout << "  size 1" << std::endl;
+      array.Shrink(1);
+      array.GetPortalControl().Set(0, OFFSET);
+      Algorithm::ScanExtended(array, array);
+      VTKM_TEST_ASSERT(array.GetNumberOfValues() == 2);
+      portal = array.GetPortalConstControl();
+      VTKM_TEST_ASSERT(portal.Get(0) == 0, "Incorrect initial value");
+      VTKM_TEST_ASSERT(portal.Get(1) == OFFSET, "Incorrect total sum");
+
+      std::cout << "  size 0" << std::endl;
+      array.Shrink(0);
+      Algorithm::ScanExtended(array, array);
+      VTKM_TEST_ASSERT(array.GetNumberOfValues() == 1);
+      portal = array.GetPortalConstControl();
+      VTKM_TEST_ASSERT(portal.Get(0) == 0, "Incorrect initial value");
+    }
+
+    std::cout << "-------------------------------------------" << std::endl;
+    std::cout << "Testing Extended Scan with multiplication operator" << std::endl;
+    {
+      std::vector<vtkm::Float64> inputValues(ARRAY_SIZE);
+      for (std::size_t i = 0; i < ARRAY_SIZE; ++i)
+      {
+        inputValues[i] = 1.01;
+      }
+
+      std::size_t mid = ARRAY_SIZE / 2;
+      inputValues[mid] = 0.0;
+
+      vtkm::cont::ArrayHandle<vtkm::Float64> array =
+        vtkm::cont::make_ArrayHandle(inputValues, vtkm::CopyFlag::On);
+
+      vtkm::Float64 initialValue = 2.00;
+      Algorithm::ScanExtended(array, array, vtkm::Multiply(), initialValue);
+
+      VTKM_TEST_ASSERT(array.GetNumberOfValues() == ARRAY_SIZE + 1,
+                       "ScanExtended output size incorrect.");
+
+      auto portal = array.GetPortalConstControl();
+      VTKM_TEST_ASSERT(portal.Get(0) == initialValue,
+                       "ScanExtended result's first value != initialValue");
+
+      for (std::size_t i = 1; i <= mid; ++i)
+      {
+        vtkm::Id index = static_cast<vtkm::Id>(i);
+        vtkm::Float64 expected = pow(1.01, static_cast<vtkm::Float64>(i)) * initialValue;
+        vtkm::Float64 got = portal.Get(index);
+        VTKM_TEST_ASSERT(test_equal(got, expected), "Incorrect results for ScanExtended");
+      }
+      for (std::size_t i = mid + 1; i < ARRAY_SIZE + 1; ++i)
+      {
+        vtkm::Id index = static_cast<vtkm::Id>(i);
+        VTKM_TEST_ASSERT(portal.Get(index) == 0.0f, "Incorrect results for ScanExtended");
+      }
+    }
+
+    std::cout << "-------------------------------------------" << std::endl;
+    std::cout << "Testing Extended Scan with a vtkm::Vec" << std::endl;
+
+    {
+      using Vec3 = vtkm::Vec3f_64;
+      using Vec3ArrayHandle = vtkm::cont::ArrayHandle<Vec3, StorageTag>;
+
+      std::vector<Vec3> testValues(ARRAY_SIZE);
+
+      for (std::size_t i = 0; i < ARRAY_SIZE; ++i)
+      {
+        testValues[i] = TestValue(1, Vec3());
+      }
+      Vec3ArrayHandle values = vtkm::cont::make_ArrayHandle(testValues, vtkm::CopyFlag::On);
+
+      Algorithm::ScanExtended(values, values);
+      VTKM_TEST_ASSERT(test_equal(vtkm::cont::ArrayGetValue(ARRAY_SIZE, values),
+                                  (TestValue(1, Vec3()) * ARRAY_SIZE)),
+                       "Got bad sum from ScanExtended");
+    }
+  }
+
   static VTKM_CONT void TestErrorExecution()
   {
     std::cout << "-------------------------------------------" << std::endl;
@@ -2763,6 +2866,7 @@ private:
       TestReduceByKeyWithFancyArrays();
 
       TestScanExclusive();
+      TestScanExtended();
 
       TestScanInclusive();
       TestScanInclusiveWithComparisonObject();

vtkm/cont/testing/UnitTestAlgorithm.cxx

@@ -99,6 +99,8 @@ void ScanTest()
   out = vtkm::cont::Algorithm::ScanExclusive(input, output, vtkm::Maximum(), vtkm::Id(0));
   vtkm::cont::Algorithm::ScanExclusiveByKey(keys, input, output, vtkm::Id(0), vtkm::Maximum());
   vtkm::cont::Algorithm::ScanExclusiveByKey(keys, input, output);
+  vtkm::cont::Algorithm::ScanExtended(input, output);
+  vtkm::cont::Algorithm::ScanExtended(input, output, vtkm::Maximum(), vtkm::Id(0));
   (void)out;
 }