Adding CUDA GPU support for Filter function

backends/cuda-ref/ceed-cuda-ref-vector.c

@@ -717,6 +717,36 @@ static int CeedVectorScale_Cuda(CeedVector x, CeedScalar alpha) {
   return CEED_ERROR_SUCCESS;
 }
 
+//------------------------------------------------------------------------------
+// Filter or clip a vector using a threshold value on the host
+//------------------------------------------------------------------------------
+static int CeedHostFilter_Cuda(CeedScalar *x_array, CeedScalar threshold, CeedSize length) {
+  CeedPragmaSIMD for (CeedSize i = 0; i < length; i++) {
+    if (fabs(x_array[i]) <= threshold) x_array[i] = 0.0;
+  }
+  return CEED_ERROR_SUCCESS;
+}
+
+//------------------------------------------------------------------------------
+// Filter or clip a vector using a threshold value on device (impl in .cu file)
+//------------------------------------------------------------------------------
+int CeedDeviceFilter_Cuda(CeedScalar *x_array, CeedScalar threshold, CeedSize length);
+
+//------------------------------------------------------------------------------
+// Filter or clip a vector using a threshold value
+//------------------------------------------------------------------------------
+static int CeedVectorFilter_Cuda(CeedVector vec, CeedScalar threshold) {
+  CeedSize         length;
+  CeedVector_Cuda *impl;
+
+  CeedCallBackend(CeedVectorGetData(vec, &impl));
+  CeedCallBackend(CeedVectorGetLength(vec, &length));
+  // Set value for synced device/host array
+  if (impl->d_array) CeedCallBackend(CeedDeviceFilter_Cuda(impl->d_array, threshold, length));
+  if (impl->h_array) CeedCallBackend(CeedHostFilter_Cuda(impl->h_array, threshold, length));
+  return CEED_ERROR_SUCCESS;
+}
+
 //------------------------------------------------------------------------------
 // Compute y = alpha x + y on the host
 //------------------------------------------------------------------------------
@@ -875,6 +905,7 @@ int CeedVectorCreate_Cuda(CeedSize n, CeedVector vec) {
   CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "Norm", CeedVectorNorm_Cuda));
   CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "Reciprocal", CeedVectorReciprocal_Cuda));
   CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "Scale", CeedVectorScale_Cuda));
+  CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "Filter", CeedVectorFilter_Cuda));
   CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "AXPY", CeedVectorAXPY_Cuda));
   CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "AXPBY", CeedVectorAXPBY_Cuda));
   CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "PointwiseMult", CeedVectorPointwiseMult_Cuda));

backends/cuda-ref/kernels/cuda-ref-vector.cu

@@ -124,6 +124,30 @@ extern "C" int CeedDeviceScale_Cuda(CeedScalar *x_array, CeedScalar alpha, CeedS
   return 0;
 }
 
+//------------------------------------------------------------------------------
+// Kernel for filter
+//------------------------------------------------------------------------------
+__global__ static void filterValueK(CeedScalar *__restrict__ x, CeedScalar threshold, CeedSize size) {
+  const CeedSize index = threadIdx.x + (CeedSize)blockDim.x * blockIdx.x;
+
+  if (index < size) {
+    if (fabs(x[index]) <= threshold) x[index] = 0.0;
+  }
+}
+
+//------------------------------------------------------------------------------
+// Filter or clip vector components to zero if their absolute value is less than or equal to threshold on device
+//------------------------------------------------------------------------------
+extern "C" int CeedDeviceFilter_Cuda(CeedScalar *x_array, CeedScalar threshold, CeedSize length) {
+  const int      block_size = 512;
+  const CeedSize vec_size   = length;
+  int            grid_size  = vec_size / block_size;
+
+  if (block_size * grid_size < vec_size) grid_size += 1;
+  filterValueK<<<grid_size, block_size>>>(x_array, threshold, length);
+  return 0;
+}
+
 //------------------------------------------------------------------------------
 // Kernel for axpy
 //------------------------------------------------------------------------------

tests/t129-vector.c

@@ -14,8 +14,16 @@ static int InitVector(CeedVector x, CeedInt len) {
   if (len <= 0) return 0;  // Nothing to set for an empty vector
 
   CeedScalar array[len];
+
   for (CeedInt i = 0; i < len; i++) array[i] = (1.0 + i) * pow(-1, i);
   CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, array);
+  {
+    // Sync memtype to device for GPU backends
+    CeedMemType type = CEED_MEM_HOST;
+
+    CeedGetPreferredMemType(CeedVectorReturnCeed(x), &type);
+    CeedVectorSyncArray(x, type);
+  }
   return 0;
 }
 
@@ -53,13 +61,6 @@ int main(int argc, char **argv) {
   CeedVectorFilter(x, tolerance);
   VerifyFilter(x, len, tolerance);
 
-  {
-    // Sync memtype to device for GPU backends
-    CeedMemType type = CEED_MEM_HOST;
-    CeedGetPreferredMemType(ceed, &type);
-    CeedVectorSyncArray(x, type);
-  }
-
   // Test Case 2 - tolerance equal to a vector value
   InitVector(x, len);
   tolerance = 7.0;