Make summations autovectorizable

Signed-off-by: Szymon Gizler <sgizler@antmicro.com>

Author: sgizler

src/gpl2/src/densityOp.cpp

@@ -35,6 +35,7 @@
 #include "densityOp.h"
 
 #include <Kokkos_Core.hpp>
+#include "kokkosUtil.h"
 
 #include "placerBase.h"
 #include "placerObjects.h"
@@ -292,10 +293,10 @@ void DensityOp::updateDensityForceBin()
 
   sumOverflow_ = 0.0;
   Kokkos::DefaultHostExecutionSpace hostSpace;
-  auto hBinOverflowArea = Kokkos::create_mirror_view_and_copy(hostSpace, dBinOverflowArea);
-  for(int i = 0; i<numBins; ++i) {
-    sumOverflow_ += hBinOverflowArea[i];
-  }
+
+  Kokkos::View<float*> hBinOverflowArea("hBinOverflowArea", dBinOverflowArea.size());
+  Kokkos::deep_copy(hBinOverflowArea, dBinOverflowArea);
+  sumOverflow_ = sumFloats(hBinOverflowArea, numBins);
 
   Kokkos::fence();
 

src/gpl2/src/kokkosUtil.h

@@ -81,3 +81,51 @@ KOKKOS_INLINE_FUNCTION float consistentCosf(float x) {
 KOKKOS_INLINE_FUNCTION float consistentExpf(float x) {
   return exp((double) x);
 }
+
+#ifdef KOKKOS_ENABLE_CUDA
+  #define HOST_FUNCTION __host__
+#else
+  #define HOST_FUNCTION KOKKOS_FUNCTION
+#endif
+
+#ifdef KOKKOS_ENABLE_CUDA
+  #define HOST_INLINE_FUNCTION inline __host__
+#else
+  #define HOST_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
+#endif
+
+// We can't use parallel_reduce as we would lose consisiency between platforms
+// In order to ensure consistency with as low performance penalty as possible, we do it with host-only functions
+// that are autovectorizable by compiler.
+HOST_INLINE_FUNCTION float sumFloats(const Kokkos::View<const float*> arr, size_t size) {
+  float partialSums[4] = {0.0, 0.0, 0.0, 0.0};
+  auto hArr = Kokkos::create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), arr);
+  for(int i = 0; i<size/4*4; i+=4) {
+    partialSums[0] += hArr[i+0];
+    partialSums[1] += hArr[i+1];
+    partialSums[2] += hArr[i+2];
+    partialSums[3] += hArr[i+3];
+  }
+  float leftover = 0.0;
+  for(int i = size/4*4; i<size; ++i) {
+    leftover += hArr[i];
+  }
+  return partialSums[0] + partialSums[1] + partialSums[2] + partialSums[3] + leftover;
+}
+
+// More accurate version of sumFloats() that use double as accumulator. TODO: Consider using Kahan summation algorithm
+HOST_INLINE_FUNCTION float sumFloatsAccurate(const Kokkos::View<const float*> arr, size_t size) {
+  auto hArr = Kokkos::create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), arr);
+  double partialSums[4] = {0.0, 0.0, 0.0, 0.0};
+  for(int i = 0; i<size/4*4; i+=4) {
+    partialSums[0] += hArr[i+0];
+    partialSums[1] += hArr[i+1];
+    partialSums[2] += hArr[i+2];
+    partialSums[3] += hArr[i+3];
+  }
+  double leftover = 0.0;
+  for(int i = size/4*4; i<size; ++i) {
+    leftover += hArr[i];
+  }
+  return partialSums[0] + partialSums[1] + partialSums[2] + partialSums[3] + leftover;
+}

src/gpl2/src/placerBase.cpp

@@ -48,6 +48,7 @@
 #include <unordered_set>
 
 #include "db_sta/dbNetwork.hh"
+#include "kokkosUtil.h"
 #include "placerObjects.h"
 #include "sta/Liberty.hh"
 #include "utl/Logger.h"
@@ -58,13 +59,6 @@ using utl::GPL2;
 
 #define REPLACE_SQRT2 1.414213562373095048801L
 
-#ifdef KOKKOS_ENABLE_CUDA
-  #define BACKEND_DEPENDENT_FUNCTION __host__
-#else
-  #define BACKEND_DEPENDENT_FUNCTION KOKKOS_FUNCTION
-#endif
-
-
 ///////////////////////////////////////////////////////////////////////////////////
 // PlacerBaseVars
 ///////////////////////////////////////////////////////////////////////////////////
@@ -844,17 +838,14 @@ void PlacerBase::initDensity1()
 
 // (a)  // (a) define the get distance method
 // getDistance is only defined on the host side
-BACKEND_DEPENDENT_FUNCTION float getDistance(const Kokkos::View<const FloatPoint*>& a,
+HOST_FUNCTION float getDistance(const Kokkos::View<const FloatPoint*>& a,
                            const Kokkos::View<const FloatPoint*>& b,
                            const int numInsts)
 {
   if (numInsts <= 0) {
     return 0.0;
   }
 
-  float sumDistance = 0.0;
-  Kokkos::DefaultHostExecutionSpace hostSpace;
-
   auto aPlusbDistance  = Kokkos::View<float*, Kokkos::DefaultExecutionSpace>("aPlusbDistance", numInsts);
   Kokkos::parallel_for(numInsts, KOKKOS_LAMBDA (const int i) {
     const FloatPoint& aPoint = a[i];
@@ -864,11 +855,7 @@ BACKEND_DEPENDENT_FUNCTION float getDistance(const Kokkos::View<const FloatPoint
     aPlusbDistance[i] = aDistance + bDistance;
   });
 
-  auto haPlusbDistance = Kokkos::create_mirror_view_and_copy(hostSpace, aPlusbDistance);
-  for(int i = 0; i<numInsts; ++i) {
-    sumDistance += haPlusbDistance[i];
-  }
-
+  float sumDistance = sumFloats(aPlusbDistance, numInsts);
   return std::sqrt(sumDistance / (2.0 * numInsts));
 }
 
@@ -881,19 +868,6 @@ struct myAbs
   }
 };
 
-BACKEND_DEPENDENT_FUNCTION float getAbsGradSum(const Kokkos::View<const float*>& a, const int numInsts)
-{
-  Kokkos::DefaultHostExecutionSpace hostSpace;
-  auto hA = Kokkos::create_mirror_view_and_copy(hostSpace, a);
-
-  double sumAbs = 0.0;
-  for(int i = 0; i<numInsts; ++i) {
-    double x = hA[i];
-    sumAbs += x;
-  }
-  return sumAbs;
-}
-
 float PlacerBase::getStepLength(const Kokkos::View<const FloatPoint*>& prevSLPCoordi,
                                 const Kokkos::View<const FloatPoint*>& prevSLPSumGrads,
                                 const Kokkos::View<const FloatPoint*>& curSLPCoordi,
@@ -1007,15 +981,14 @@ void PlacerBase::updateGradients(const Kokkos::View<float*>& wireLengthGradients
   densityGradSum_ = 0;
 
   // get the forces on each instance
-  Kokkos::View<float*> wirelenabsGradXPlusY("absGradXPlusY", numInsts_);
-  Kokkos::View<float*> densityabsGradXPlusY("absGradXPlusY", numInsts_);
-  Kokkos::DefaultHostExecutionSpace hostSpace;
+  Kokkos::View<float*> wireLengthGradAbsXPlusY("wireLengthGradAbsXPlusY", numInsts_);
+  Kokkos::View<float*> densityGradAbsXPlusY("densityGradAbsXPlusY", numInsts_);
 
-  getWireLengthGradientWA(wireLengthGradientsX, wireLengthGradientsY, wirelenabsGradXPlusY);
-  getDensityGradient(densityGradientsX, densityGradientsY, densityabsGradXPlusY);
+  getWireLengthGradientWA(wireLengthGradientsX, wireLengthGradientsY, wireLengthGradAbsXPlusY);
+  getDensityGradient(densityGradientsX, densityGradientsY, densityGradAbsXPlusY);
 
-  wireLengthGradSum_ += getAbsGradSum(wirelenabsGradXPlusY, numInsts_);
-  densityGradSum_ += getAbsGradSum(densityabsGradXPlusY, numInsts_);
+  wireLengthGradSum_ += sumFloatsAccurate(wireLengthGradAbsXPlusY, numInsts_);
+  densityGradSum_ += sumFloatsAccurate(densityGradAbsXPlusY, numInsts_);
 
   sumGradientKernel(numInsts_,
                                                densityPenalty_,