Adding filter function

include/ceed-impl.h

@@ -159,6 +159,7 @@ struct CeedVector_private {
   int (*RestoreArrayRead)(CeedVector);
   int (*Norm)(CeedVector, CeedNormType, CeedScalar *);
   int (*Scale)(CeedVector, CeedScalar);
+  int (*Filter)(CeedVector, CeedScalar);
   int (*AXPY)(CeedVector, CeedScalar, CeedVector);
   int (*AXPBY)(CeedVector, CeedScalar, CeedScalar, CeedVector);
   int (*PointwiseMult)(CeedVector, CeedVector, CeedVector);

include/ceed/ceed.h

@@ -212,6 +212,7 @@ CEED_EXTERN int  CeedVectorRestoreArray(CeedVector vec, CeedScalar **array);
 CEED_EXTERN int  CeedVectorRestoreArrayRead(CeedVector vec, const CeedScalar **array);
 CEED_EXTERN int  CeedVectorNorm(CeedVector vec, CeedNormType type, CeedScalar *norm);
 CEED_EXTERN int  CeedVectorScale(CeedVector x, CeedScalar alpha);
+CEED_EXTERN int  CeedVectorFilter(CeedVector x, CeedScalar threshold);
 CEED_EXTERN int  CeedVectorAXPY(CeedVector y, CeedScalar alpha, CeedVector x);
 CEED_EXTERN int  CeedVectorAXPBY(CeedVector y, CeedScalar alpha, CeedScalar beta, CeedVector x);
 CEED_EXTERN int  CeedVectorPointwiseMult(CeedVector w, CeedVector x, CeedVector y);

interface/ceed-vector.c

@@ -756,6 +756,46 @@ int CeedVectorScale(CeedVector x, CeedScalar alpha) {
   return CEED_ERROR_SUCCESS;
 }
 
+/**
+  @brief Filters or clips a `CeedVector` using a threshold value. 
+  All entries in `x` with an absolute value less than or equal to `threshold` are set to `0.0`.
+
+  @param[in,out] x         `CeedVector` to filter
+  @param[in]     threshold clipping threshold or tolerance
+
+  @return An error code: 0 - success, otherwise - failure
+
+  @ref User
+**/
+int CeedVectorFilter(CeedVector x, CeedScalar threshold) {
+  bool        has_valid_array = true;
+  CeedSize    length;
+  CeedScalar *x_array = NULL;
+
+  CeedCall(CeedVectorHasValidArray(x, &has_valid_array));
+  CeedCheck(has_valid_array, CeedVectorReturnCeed(x), CEED_ERROR_BACKEND,
+            "CeedVector has no valid data to filter, must set data with CeedVectorSetValue or CeedVectorSetArray");
+
+  // Return early for empty vector
+  CeedCall(CeedVectorGetLength(x, &length));
+  if (length == 0) return CEED_ERROR_SUCCESS;
+
+  // Backend implementation
+  if (x->Filter) {
+    CeedCall(x->Filter(x, threshold));
+    return CEED_ERROR_SUCCESS;
+  }
+
+  // Default implementation
+  CeedCall(CeedVectorGetArray(x, CEED_MEM_HOST, &x_array));
+  assert(x_array);
+  CeedPragmaSIMD for (CeedSize i = 0; i < length; i++) {
+    if (fabs(x_array[i]) <= threshold) x_array[i] = 0.0;
+  }
+  CeedCall(CeedVectorRestoreArray(x, &x_array));
+  return CEED_ERROR_SUCCESS;
+}
+
 /**
   @brief Compute `y = alpha x + y`
 

interface/ceed.c

@@ -1306,6 +1306,7 @@ int CeedInit(const char *resource, Ceed *ceed) {
       CEED_FTABLE_ENTRY(CeedVector, RestoreArrayRead),
       CEED_FTABLE_ENTRY(CeedVector, Norm),
       CEED_FTABLE_ENTRY(CeedVector, Scale),
+      CEED_FTABLE_ENTRY(CeedVector, Filter),
       CEED_FTABLE_ENTRY(CeedVector, AXPY),
       CEED_FTABLE_ENTRY(CeedVector, AXPBY),
       CEED_FTABLE_ENTRY(CeedVector, PointwiseMult),

tests/t129-vector.c

@@ -0,0 +1,85 @@
+/// @file
+/// Test filtering a vector using a threshold or tolerance
+/// \test Test filtering a vector using a threshold or tolerance
+
+//TESTARGS(name="vector of length 10") {ceed_resource} 10
+//TESTARGS(name="empty vector") {ceed_resource} 0
+#include <ceed.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+// Test builds the vector [1.0, -2.0, 3.0, -4.0, 5.0, -6.0, 7.0, -8.0, 9.0, -10.0]
+// and filters it using different tolerances
+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);
+  return 0;
+}
+
+static int VerifyFilter(CeedVector x, CeedInt len, CeedScalar tolerance) {
+  const CeedScalar *read_array;
+
+  CeedVectorGetArrayRead(x, CEED_MEM_HOST, &read_array);
+  for (CeedInt i = 0; i < len; i++) {
+    CeedScalar initial_value = (1.0 + i) * pow(-1, i);
+    CeedScalar expected_value = (fabs(initial_value) <= tolerance) ? 0.0 : initial_value;
+
+    if (fabs(read_array[i] - expected_value) > 1e-14) {
+      // LCOV_EXCL_START
+      printf("Error in filtered vector at index %" CeedInt_FMT ", computed: %f actual: %f\n", i, read_array[i], expected_value);
+      // LCOV_EXCL_STOP
+    }
+    }
+  CeedVectorRestoreArrayRead(x, &read_array);
+  return 0;
+}
+
+int main(int argc, char **argv) {
+  Ceed       ceed;
+  CeedVector x;
+  CeedInt    len       = 10;
+  CeedScalar tolerance = 1e-10; 
+
+  CeedInit(argv[1], &ceed);
+  len = argc > 2 ? atoi(argv[2]) : len;
+
+  CeedVectorCreate(ceed, len, &x);
+
+  // Test Case 1 - tolerance between vector values
+  InitVector(x, len);
+  tolerance = 3.5;
+  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;
+  CeedVectorFilter(x, tolerance);
+  VerifyFilter(x, len, tolerance);
+
+  // Test Case 3 - filter no values
+  InitVector(x, len);
+  tolerance = 0.0;
+  CeedVectorFilter(x, tolerance);
+  VerifyFilter(x, len, tolerance);
+
+  // Test Case 4 - filter all values
+  InitVector(x, len);
+  tolerance = 100.0;
+  CeedVectorFilter(x, tolerance);
+  VerifyFilter(x, len, tolerance);
+
+  CeedVectorDestroy(&x);
+  CeedDestroy(&ceed);
+  return 0;
+}