 			return 0.08617330350;
 		}
 
+
+		/////////////////////////////////////////////////////////////////
+        // BOILERPLATE CUDA Reductions
+
+        __inline__ __device__
+        scalar warpReduceSum(scalar val)
+        {
+            for (int offset = warpSize/2; offset > 0; offset /= 2) 
+                val += __shfl_down(val, offset);
+            return val;
+        }
+
+        __inline__ __device__
+        scalar blockReduceSum(scalar val)
+        {
+            static __shared__ scalar shared[32]; // Shared mem for 32 partial sums
+            int lane = threadIdx.x % warpSize;
+            int wid = threadIdx.x / warpSize;
+
+            val = warpReduceSum(val);     // Each warp performs partial reduction
+
+            if (lane==0) shared[wid]=val; // Write reduced value to shared memory
+
+            __syncthreads();              // Wait for all partial reductions
+
+            //read from shared memory only if that warp existed
+            val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0;
+
+            if (wid==0) val = warpReduceSum(val); //Final reduce within first warp
+
+            return val;
+        }
+
+        __global__ void cu_sum(const scalar *in, scalar* out, int N)
+        {
+            scalar sum = int(0);
+            for(int i = blockIdx.x * blockDim.x + threadIdx.x; 
+                i < N; 
+                i += blockDim.x * gridDim.x)
+            {
+                sum += in[i];
+            }
+            sum = blockReduceSum(sum);
+            if (threadIdx.x == 0)
+                atomicAdd(out, sum);
+        }
+
+
+
+        __inline__ __device__
+        Vector3 warpReduceSum(Vector3 val)
+        {
+            for (int offset = warpSize/2; offset > 0; offset /= 2)
+            {
+                val[0] += __shfl_down(val[0], offset);
+                val[1] += __shfl_down(val[1], offset);
+                val[2] += __shfl_down(val[2], offset);
+            }
+            return val;
+        }
+
+        __inline__ __device__
+        Vector3 blockReduceSum(Vector3 val)
+        {
+            static __shared__ Vector3 shared[32]; // Shared mem for 32 partial sums
+            int lane = threadIdx.x % warpSize;
+            int wid = threadIdx.x / warpSize;
+
+            val = warpReduceSum(val);     // Each warp performs partial reduction
+
+            if (lane==0) shared[wid]=val; // Write reduced value to shared memory
+
+            __syncthreads();              // Wait for all partial reductions
+
+            // Read from shared memory only if that warp existed
+            val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : Vector3{0,0,0};
+
+            if (wid==0) val = warpReduceSum(val); //Final reduce within first warp
+
+            return val;
+        }
+
+        __global__ void cu_sum(const Vector3 *in, Vector3* out, int N)
+        {
+            Vector3 sum{0,0,0};
+            for(int i = blockIdx.x * blockDim.x + threadIdx.x; 
+                i < N; 
+                i += blockDim.x * gridDim.x)
+            {
+                sum += in[i];
+            }
+            sum = blockReduceSum(sum);
+            if (threadIdx.x == 0)
+            {
+                atomicAdd(&out[0][0], sum[0]);
+                atomicAdd(&out[0][1], sum[1]);
+                atomicAdd(&out[0][2], sum[2]);
+            }
+        }
 
 		/////////////////////////////////////////////////////////////////
 

            cudaDeviceSynchronize();
        }

        __global__ void cu_sum(const Vector3 *sf, Vector3 *out, size_t N)
        {
            int idx = blockIdx.x * blockDim.x + threadIdx.x;
            if(idx < N)
            {
                atomicAdd(&out[0][0], sf[idx][0]);
                atomicAdd(&out[0][1], sf[idx][1]);
                atomicAdd(&out[0][2], sf[idx][2]);
            }
        }
        Vector3 sum(const vectorfield & vf)
		{
			Vector3 ret = { 0,0,0 };
            int n = vf.size();
            cu_sum<<<(n+1023)/1024, 1024>>>(vf.data(), &ret, n);
            cudaDeviceSynchronize();
			return ret;
		}

		Vector3 mean(const vectorfield & vf)
		{
			Vector3 ret = { 0,0,0 };
            int n = vf.size();
            cu_sum<<<(n+1023)/1024, 1024>>>(vf.data(), &ret, n);
            cudaDeviceSynchronize();
            ret /= (scalar)n;
			return ret;
		}



            int idx = blockIdx.x * blockDim.x + threadIdx.x;
            if(idx < N)
            {
                sf[idx] *= s;
            }
        }
		void scale(scalarfield & sf, scalar s)
