diff --git a/include/ceed-impl.h b/include/ceed-impl.h index e5f8773f37..a627cf2528 100644 --- a/include/ceed-impl.h +++ b/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); diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h index a76b9238c3..3ff8293c12 100644 --- a/include/ceed/ceed.h +++ b/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); diff --git a/interface/ceed-vector.c b/interface/ceed-vector.c index eb9f4fc85b..90b9662b39 100644 --- a/interface/ceed-vector.c +++ b/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` diff --git a/interface/ceed.c b/interface/ceed.c index 2088980c8e..c1b7180faf 100644 --- a/interface/ceed.c +++ b/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), diff --git a/tests/t129-vector.c b/tests/t129-vector.c new file mode 100644 index 0000000000..02f257d4bf --- /dev/null +++ b/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 +#include +#include +#include +// 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; +}