From f30b1135883847d0dd65c557f1d0ec6354b82b2b Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Mon, 17 Apr 2023 15:33:38 -0700 Subject: [PATCH 1/6] Add CeedElemRestrictionApplyUnsigned for operator diagonal assembly in the case of an oriented element restriction --- backends/ref/ceed-ref-restriction.c | 100 ++++++++++++++++++---------- backends/ref/ceed-ref.h | 2 +- include/ceed-impl.h | 1 + include/ceed/ceed.h | 2 + interface/ceed-elemrestriction.c | 46 +++++++++++++ interface/ceed-preconditioning.c | 4 +- interface/ceed.c | 1 + 7 files changed, 120 insertions(+), 36 deletions(-) diff --git a/backends/ref/ceed-ref-restriction.c b/backends/ref/ceed-ref-restriction.c index deabe0b0f1..3f0c282649 100644 --- a/backends/ref/ceed-ref-restriction.c +++ b/backends/ref/ceed-ref-restriction.c @@ -16,7 +16,7 @@ // Core ElemRestriction Apply Code //------------------------------------------------------------------------------ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { CeedElemRestriction_Ref *impl; CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); @@ -79,8 +79,15 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, const for (CeedInt e = start * blk_size; e < stop * blk_size; e += blk_size) { CeedPragmaSIMD for (CeedInt k = 0; k < num_comp; k++) { CeedPragmaSIMD for (CeedInt i = 0; i < elem_size * blk_size; i++) { - vv[elem_size * (k * blk_size + num_comp * e) + i - v_offset] = - uu[impl->offsets[i + elem_size * e] + k * comp_stride] * (is_oriented && impl->orient[i + elem_size * e] ? -1. : 1.); + if (!use_orient || !impl->orient) { + // Unsigned restriction + vv[elem_size * (k * blk_size + num_comp * e) + i - v_offset] = + uu[impl->offsets[i + elem_size * e] + k * comp_stride]; + } else { + // Signed restriction + vv[elem_size * (k * blk_size + num_comp * e) + i - v_offset] = + uu[impl->offsets[i + elem_size * e] + k * comp_stride] * (impl->orient[i + elem_size * e] ? -1.0 : 1.0); + } } } } @@ -129,8 +136,15 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, const for (CeedInt i = 0; i < elem_size * blk_size; i += blk_size) { // Iteration bound set to discard padding elements for (CeedInt j = i; j < i + CeedIntMin(blk_size, num_elem - e); j++) { - vv[impl->offsets[j + e * elem_size] + k * comp_stride] += - uu[elem_size * (k * blk_size + num_comp * e) + j - v_offset] * (is_oriented && impl->orient[j + e * elem_size] ? -1. : 1.); + if (!use_orient || !impl->orient) { + // Unsigned restriction + vv[impl->offsets[j + e * elem_size] + k * comp_stride] += + uu[elem_size * (k * blk_size + num_comp * e) + j - v_offset]; + } else { + // Signed restriction + vv[impl->offsets[j + e * elem_size] + k * comp_stride] += + uu[elem_size * (k * blk_size + num_comp * e) + j - v_offset] * (impl->orient[j + e * elem_size] ? -1.0 : 1.0); + } } } } @@ -147,67 +161,67 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, const // ElemRestriction Apply - Common Sizes //------------------------------------------------------------------------------ static int CeedElemRestrictionApply_Ref_110(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 1, 1, comp_stride, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 1, 1, comp_stride, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_111(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 1, 1, 1, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 1, 1, 1, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_180(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 1, 8, comp_stride, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 1, 8, comp_stride, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_181(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 1, 8, 1, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 1, 8, 1, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_310(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 3, 1, comp_stride, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 3, 1, comp_stride, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_311(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 3, 1, 1, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 3, 1, 1, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_380(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 3, 8, comp_stride, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 3, 8, comp_stride, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_381(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 3, 8, 1, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 3, 8, 1, start, stop, use_orient, t_mode, u, v, request); } // LCOV_EXCL_START static int CeedElemRestrictionApply_Ref_510(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 5, 1, comp_stride, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 5, 1, comp_stride, start, stop, use_orient, t_mode, u, v, request); } // LCOV_EXCL_STOP static int CeedElemRestrictionApply_Ref_511(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 5, 1, 1, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 5, 1, 1, start, stop, use_orient, t_mode, u, v, request); } // LCOV_EXCL_START static int CeedElemRestrictionApply_Ref_580(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 5, 8, comp_stride, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 5, 8, comp_stride, start, stop, use_orient, t_mode, u, v, request); } // LCOV_EXCL_STOP static int CeedElemRestrictionApply_Ref_581(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { - return CeedElemRestrictionApply_Ref_Core(r, 5, 8, 1, start, stop, t_mode, u, v, request); + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + return CeedElemRestrictionApply_Ref_Core(r, 5, 8, 1, start, stop, use_orient, t_mode, u, v, request); } //------------------------------------------------------------------------------ @@ -222,7 +236,22 @@ static int CeedElemRestrictionApply_Ref(CeedElemRestriction r, CeedTransposeMode CeedElemRestriction_Ref *impl; CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); - return impl->Apply(r, num_comp, blk_size, comp_stride, 0, num_blk, t_mode, u, v, request); + return impl->Apply(r, num_comp, blk_size, comp_stride, 0, num_blk, true, t_mode, u, v, request); +} + +//------------------------------------------------------------------------------ +// ElemRestriction Apply Unsigned +//------------------------------------------------------------------------------ +static int CeedElemRestrictionApplyUnsigned_Ref(CeedElemRestriction r, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt num_blk, blk_size, num_comp, comp_stride; + CeedCallBackend(CeedElemRestrictionGetNumBlocks(r, &num_blk)); + CeedCallBackend(CeedElemRestrictionGetBlockSize(r, &blk_size)); + CeedCallBackend(CeedElemRestrictionGetNumComponents(r, &num_comp)); + CeedCallBackend(CeedElemRestrictionGetCompStride(r, &comp_stride)); + CeedElemRestriction_Ref *impl; + CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); + + return impl->Apply(r, num_comp, blk_size, comp_stride, 0, num_blk, false, t_mode, u, v, request); } //------------------------------------------------------------------------------ @@ -237,7 +266,7 @@ static int CeedElemRestrictionApplyBlock_Ref(CeedElemRestriction r, CeedInt bloc CeedElemRestriction_Ref *impl; CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); - return impl->Apply(r, num_comp, blk_size, comp_stride, block, block + 1, t_mode, u, v, request); + return impl->Apply(r, num_comp, blk_size, comp_stride, block, block + 1, true, t_mode, u, v, request); } //------------------------------------------------------------------------------ @@ -329,6 +358,7 @@ int CeedElemRestrictionCreate_Ref(CeedMemType mem_type, CeedCopyMode copy_mode, CeedInt layout[3] = {1, elem_size, elem_size * num_comp}; CeedCallBackend(CeedElemRestrictionSetELayout(r, layout)); CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", CeedElemRestrictionApply_Ref)); + CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyUnsigned", CeedElemRestrictionApplyUnsigned_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyBlock", CeedElemRestrictionApplyBlock_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "GetOffsets", CeedElemRestrictionGetOffsets_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", CeedElemRestrictionDestroy_Ref)); @@ -390,15 +420,16 @@ int CeedElemRestrictionCreate_Ref(CeedMemType mem_type, CeedCopyMode copy_mode, //------------------------------------------------------------------------------ int CeedElemRestrictionCreateOriented_Ref(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const bool *orient, CeedElemRestriction r) { - CeedElemRestriction_Ref *impl; - CeedInt num_elem, elem_size; // Set up for normal restriction with explicit offsets. This sets up dispatch to // CeedElemRestrictionApply_Ref_* and manages the impl->offsets array copy/allocation. + CeedElemRestriction_Ref *impl; + CeedInt num_elem, elem_size; CeedCallBackend(CeedElemRestrictionCreate_Ref(mem_type, copy_mode, offsets, r)); - CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); CeedCallBackend(CeedElemRestrictionGetNumElements(r, &num_elem)); CeedCallBackend(CeedElemRestrictionGetElementSize(r, &elem_size)); + + // Copy data switch (copy_mode) { case CEED_COPY_VALUES: CeedCallBackend(CeedMalloc(num_elem * elem_size, &impl->orient_allocated)); @@ -412,6 +443,7 @@ int CeedElemRestrictionCreateOriented_Ref(CeedMemType mem_type, CeedCopyMode cop case CEED_USE_POINTER: impl->orient = orient; } + return CEED_ERROR_SUCCESS; } diff --git a/backends/ref/ceed-ref.h b/backends/ref/ceed-ref.h index 168ba273be..41acf27771 100644 --- a/backends/ref/ceed-ref.h +++ b/backends/ref/ceed-ref.h @@ -30,7 +30,7 @@ typedef struct { // Orientation, if it exists, is true when the face must be flipped (multiplies by -1.). const bool *orient; bool *orient_allocated; - int (*Apply)(CeedElemRestriction, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, CeedTransposeMode, CeedVector, CeedVector, CeedRequest *); + int (*Apply)(CeedElemRestriction, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt, bool, CeedTransposeMode, CeedVector, CeedVector, CeedRequest *); } CeedElemRestriction_Ref; typedef struct { diff --git a/include/ceed-impl.h b/include/ceed-impl.h index ba2f465bf3..9101814d82 100644 --- a/include/ceed-impl.h +++ b/include/ceed-impl.h @@ -151,6 +151,7 @@ struct CeedVector_private { struct CeedElemRestriction_private { Ceed ceed; int (*Apply)(CeedElemRestriction, CeedTransposeMode, CeedVector, CeedVector, CeedRequest *); + int (*ApplyUnsigned)(CeedElemRestriction, CeedTransposeMode, CeedVector, CeedVector, CeedRequest *); int (*ApplyBlock)(CeedElemRestriction, CeedInt, CeedTransposeMode, CeedVector, CeedVector, CeedRequest *); int (*GetOffsets)(CeedElemRestriction, CeedMemType, const CeedInt **); int (*Destroy)(CeedElemRestriction); diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h index df23ff2798..6f5185df1d 100644 --- a/include/ceed/ceed.h +++ b/include/ceed/ceed.h @@ -243,6 +243,8 @@ CEED_EXTERN int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_e CEED_EXTERN int CeedElemRestrictionReferenceCopy(CeedElemRestriction rstr, CeedElemRestriction *rstr_copy); CEED_EXTERN int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *lvec, CeedVector *evec); CEED_EXTERN int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request); +CEED_EXTERN int CeedElemRestrictionApplyUnsigned(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, + CeedRequest *request); CEED_EXTERN int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request); CEED_EXTERN int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed); diff --git a/interface/ceed-elemrestriction.c b/interface/ceed-elemrestriction.c index 9deb75da3d..11bf0c34c1 100644 --- a/interface/ceed-elemrestriction.c +++ b/interface/ceed-elemrestriction.c @@ -631,6 +631,52 @@ int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, return CEED_ERROR_SUCCESS; } +/** + @brief Restrict an L-vector to an E-vector or apply its transpose ignoring any + provided orientations + + @param[in] rstr CeedElemRestriction + @param[in] t_mode Apply restriction or transpose + @param[in] u Input vector (of size @a l_size when t_mode=@ref CEED_NOTRANSPOSE) + @param[out] ru Output vector (of shape [@a num_elem * @a elem_size] when t_mode=@ref CEED_NOTRANSPOSE). + Ordering of the e-vector is decided by the backend. + @param[in] request Request or @ref CEED_REQUEST_IMMEDIATE + + @return An error code: 0 - success, otherwise - failure + + @ref User +**/ +int CeedElemRestrictionApplyUnsigned(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) { + CeedInt m, n; + + if (t_mode == CEED_NOTRANSPOSE) { + m = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; + n = rstr->l_size; + } else { + m = rstr->l_size; + n = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; + } + if (n != u->length) { + // LCOV_EXCL_START + return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, + "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", u->length, m, + n); + // LCOV_EXCL_STOP + } + if (m != ru->length) { + // LCOV_EXCL_START + return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, + "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", ru->length, + m, n); + // LCOV_EXCL_STOP + } + if (rstr->num_elem > 0) { + if (rstr->ApplyUnsigned) CeedCall(rstr->ApplyUnsigned(rstr, t_mode, u, ru, request)); + else CeedCall(rstr->Apply(rstr, t_mode, u, ru, request)); + } + return CEED_ERROR_SUCCESS; +} + /** @brief Restrict an L-vector to a block of an E-vector or apply its transpose diff --git a/interface/ceed-preconditioning.c b/interface/ceed-preconditioning.c index d10eaceded..ded55ddca5 100644 --- a/interface/ceed-preconditioning.c +++ b/interface/ceed-preconditioning.c @@ -396,7 +396,7 @@ static inline int CeedSingleOperatorAssembleAddDiagonal_Core(CeedOperator op, Ce CeedCall(CeedVectorRestoreArray(elem_diag, &elem_diag_array)); // Assemble local operator diagonal - CeedCall(CeedElemRestrictionApply(diag_elem_rstr, CEED_TRANSPOSE, elem_diag, assembled, request)); + CeedCall(CeedElemRestrictionApplyUnsigned(diag_elem_rstr, CEED_TRANSPOSE, elem_diag, assembled, request)); // Cleanup if (is_pointblock) CeedCall(CeedElemRestrictionDestroy(&diag_elem_rstr)); @@ -442,6 +442,8 @@ static inline int CeedCompositeOperatorLinearAssembleAddDiagonal(CeedOperator op Users should generally use CeedOperatorLinearAssembleSymbolic() + Note: For operators using oriented element restrictions, entries in rows or cols may be negative indicating the assembled value at this nonzero should be negated + @param[in] op CeedOperator to assemble nonzero pattern @param[in] offset Offset for number of entries @param[out] rows Row number for each entry diff --git a/interface/ceed.c b/interface/ceed.c index b8fae07d5e..852d52c359 100644 --- a/interface/ceed.c +++ b/interface/ceed.c @@ -814,6 +814,7 @@ int CeedInit(const char *resource, Ceed *ceed) { CEED_FTABLE_ENTRY(CeedVector, Reciprocal), CEED_FTABLE_ENTRY(CeedVector, Destroy), CEED_FTABLE_ENTRY(CeedElemRestriction, Apply), + CEED_FTABLE_ENTRY(CeedElemRestriction, ApplyUnsigned), CEED_FTABLE_ENTRY(CeedElemRestriction, ApplyBlock), CEED_FTABLE_ENTRY(CeedElemRestriction, GetOffsets), CEED_FTABLE_ENTRY(CeedElemRestriction, Destroy), From 7e2def383022a8dc258779d80e0e212f0ebe40f9 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Mon, 17 Apr 2023 15:34:32 -0700 Subject: [PATCH 2/6] Cherry-pick H(div) mass matrix diagonal assembly test from 02773612 --- tests/README.md | 3 +- tests/t570-operator.c | 146 ++++++++++++++++++++++++++++++++++++++++++ tests/t570-operator.h | 73 +++++++++++++++++++++ 3 files changed, 221 insertions(+), 1 deletion(-) create mode 100644 tests/t570-operator.c create mode 100644 tests/t570-operator.h diff --git a/tests/README.md b/tests/README.md index e461403046..a7a5600a7b 100644 --- a/tests/README.md +++ b/tests/README.md @@ -25,4 +25,5 @@ The tests are organized by API object, and some tests are further organized, as 5.3. CeedOperator diagonal and CeedQFunction assembly tests 5.4. CeedOperator element inverse tests 5.5. CeedOperator multigrid level tests - 5.6. CeedOperator full assembly tests + 5.6. CeedOperator full assembly tests + 5.7. CeedOperator extra diagonal and full assembly tests diff --git a/tests/t570-operator.c b/tests/t570-operator.c new file mode 100644 index 0000000000..20eb3b4336 --- /dev/null +++ b/tests/t570-operator.c @@ -0,0 +1,146 @@ +/// @file +/// Test assembly of H(div) mass matrix operator diagonal +/// \test Test assembly of H(div) mass matrix operator diagonal +#include +#include +#include +#include + +#include "t330-basis.h" +#include "t570-operator.h" + +int main(int argc, char **argv) { + Ceed ceed; + CeedElemRestriction elem_restriction_x, elem_restriction_u; + CeedBasis basis_x, basis_u; + CeedQFunction qf_mass; + CeedOperator op_mass; + CeedVector x, assembled, u, v; + CeedInt dim = 2, p = 8, q = 3, px = 2; + CeedInt n_x = 1, n_y = 1; // Currently only implemented for single element + CeedInt row, column, offset; + CeedInt num_elem = n_x * n_y, num_faces = (n_x + 1) * n_y + (n_y + 1) * n_x; + CeedInt num_dofs_x = (n_x + 1) * (n_y + 1), num_dofs_u = num_faces * 2, num_qpts = q * q; + CeedInt ind_x[num_elem * px * px], ind_u[num_elem * p]; + bool orient_u[num_elem * p]; + CeedScalar assembled_true[num_dofs_u]; + CeedScalar q_ref[dim * num_qpts], q_weight[num_qpts]; + CeedScalar interp[dim * p * num_qpts], div[p * num_qpts]; + + CeedInit(argv[1], &ceed); + + // Vectors + CeedVectorCreate(ceed, dim * num_dofs_x, &x); + { + CeedScalar x_array[dim * num_dofs_x]; + + for (CeedInt i = 0; i < n_x + 1; i++) { + for (CeedInt j = 0; j < n_y + 1; j++) { + x_array[i + j * (n_x + 1) + 0 * num_dofs_x] = i / (CeedScalar)n_x; + x_array[i + j * (n_x + 1) + 1 * num_dofs_x] = j / (CeedScalar)n_y; + } + } + CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); + } + CeedVectorCreate(ceed, num_dofs_u, &u); + CeedVectorCreate(ceed, num_dofs_u, &v); + + // Restrictions + for (CeedInt i = 0; i < num_elem; i++) { + column = i % n_x; + row = i / n_x; + offset = column * (px - 1) + row * (n_x + 1) * (px - 1); + + for (CeedInt j = 0; j < px; j++) { + for (CeedInt k = 0; k < px; k++) { + ind_x[px * (px * i + k) + j] = offset + k * (n_x + 1) + j; + } + } + } + bool orient_u_local[8] = {false, false, false, false, true, true, true, true}; + CeedInt ind_u_local[8] = {0, 1, 6, 7, 2, 3, 4, 5}; + for (CeedInt j = 0; j < n_y; j++) { + for (CeedInt i = 0; i < n_x; i++) { + for (CeedInt k = 0; k < p; k++) { + ind_u[k] = ind_u_local[k]; + orient_u[k] = orient_u_local[k]; + } + } + } + CeedElemRestrictionCreate(ceed, num_elem, px * px, dim, num_dofs_x, dim * num_dofs_x, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x); + CeedElemRestrictionCreateOriented(ceed, num_elem, p, 1, 1, num_dofs_u, CEED_MEM_HOST, CEED_COPY_VALUES, ind_u, orient_u, &elem_restriction_u); + + // Bases + CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, px, q, CEED_GAUSS, &basis_x); + + BuildHdivQuadrilateral(q, q_ref, q_weight, interp, div, CEED_GAUSS); + CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, 1, p, num_qpts, interp, div, q_ref, q_weight, &basis_u); + + // QFunctions + CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass); + CeedQFunctionAddInput(qf_mass, "weight", 1, CEED_EVAL_WEIGHT); + CeedQFunctionAddInput(qf_mass, "dx", dim * dim, CEED_EVAL_GRAD); + CeedQFunctionAddInput(qf_mass, "u", dim, CEED_EVAL_INTERP); + CeedQFunctionAddOutput(qf_mass, "v", dim, CEED_EVAL_INTERP); + + // Operators + CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass); + CeedOperatorSetField(op_mass, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); + CeedOperatorSetField(op_mass, "dx", elem_restriction_x, basis_x, x); + CeedOperatorSetField(op_mass, "u", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_mass, "v", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); + + // Assemble diagonal + CeedVectorCreate(ceed, num_dofs_u, &assembled); + CeedOperatorLinearAssembleDiagonal(op_mass, assembled, CEED_REQUEST_IMMEDIATE); + + // Manually assemble diagonal + CeedVectorSetValue(u, 0.0); + for (int i = 0; i < num_dofs_u; i++) { + CeedScalar *u_array; + const CeedScalar *v_array; + + // Set input + CeedVectorGetArray(u, CEED_MEM_HOST, &u_array); + u_array[i] = 1.0; + if (i) u_array[i - 1] = 0.0; + CeedVectorRestoreArray(u, &u_array); + + // Compute diag entry for DoF i + CeedOperatorApply(op_mass, u, v, CEED_REQUEST_IMMEDIATE); + + // Retrieve entry + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + assembled_true[i] = v_array[i]; + CeedVectorRestoreArrayRead(v, &v_array); + } + + // Check output + { + const CeedScalar *assembled_array; + + CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array); + for (int i = 0; i < num_dofs_u; i++) { + if (fabs(assembled_array[i] - assembled_true[i]) > 100. * CEED_EPSILON) { + // LCOV_EXCL_START + printf("[%" CeedInt_FMT "] Error in assembly: %f != %f\n", i, assembled_array[i], assembled_true[i]); + // LCOV_EXCL_STOP + } + } + CeedVectorRestoreArrayRead(assembled, &assembled_array); + } + + // Cleanup + CeedVectorDestroy(&x); + CeedVectorDestroy(&assembled); + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedElemRestrictionDestroy(&elem_restriction_u); + CeedElemRestrictionDestroy(&elem_restriction_x); + CeedBasisDestroy(&basis_u); + CeedBasisDestroy(&basis_x); + CeedQFunctionDestroy(&qf_mass); + CeedOperatorDestroy(&op_mass); + CeedDestroy(&ceed); + return 0; +} diff --git a/tests/t570-operator.h b/tests/t570-operator.h new file mode 100644 index 0000000000..cab7f4865b --- /dev/null +++ b/tests/t570-operator.h @@ -0,0 +1,73 @@ +// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. +// All Rights Reserved. See the top-level LICENSE and NOTICE files for details. +// +// SPDX-License-Identifier: BSD-2-Clause +// +// This file is part of CEED: http://github.com/ceed + +#include + +// Compute det(A) +CEED_QFUNCTION_HELPER CeedScalar MatDet2x2(const CeedScalar A[2][2]) { + return A[0][0] * A[1][1] - A[1][0] * A[0][1]; +} + +// Compute alpha * A^T * B = C +CEED_QFUNCTION_HELPER int AlphaMatTransposeMatMult2x2(const CeedScalar alpha, const CeedScalar A[2][2], const CeedScalar B[2][2], + CeedScalar C[2][2]) { + for (CeedInt j = 0; j < 2; j++) { + for (CeedInt k = 0; k < 2; k++) { + C[j][k] = 0; + for (CeedInt m = 0; m < 2; m++) { + C[j][k] += alpha * A[m][j] * B[m][k]; + } + } + } + + return 0; +} + +// Compute matrix-vector product: alpha*A*u +CEED_QFUNCTION_HELPER int AlphaMatVecMult2x2(const CeedScalar alpha, const CeedScalar A[2][2], const CeedScalar u[2], CeedScalar v[2]) { + // Compute v = alpha*A*u + for (CeedInt k = 0; k < 2; k++) { + v[k] = 0; + for (CeedInt m = 0; m < 2; m++) v[k] += A[k][m] * u[m] * alpha; + } + + return 0; +}; + +CEED_QFUNCTION(mass)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { + // Inputs + const CeedScalar(*w) = in[0], (*dxdX)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[1], + (*u)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; + + // Outputs + CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; + + // Quadrature Point Loop + CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { + // Setup, J = dx/dX + const CeedScalar J[2][2] = { + {dxdX[0][0][i], dxdX[1][0][i]}, + {dxdX[0][1][i], dxdX[1][1][i]} + }; + const CeedScalar det_J = MatDet2x2(J); + + CeedScalar u1[2] = {u[0][i], u[1][i]}, v1[2]; + // *INDENT-ON* + // Piola map: J^T*J*u*w/detJ + // 1) Compute J^T * J + CeedScalar JT_J[2][2]; + AlphaMatTransposeMatMult2x2(1., J, J, JT_J); + + // 2) Compute J^T*J*u * w /detJ + AlphaMatVecMult2x2(w[i] / det_J, JT_J, u1, v1); + for (CeedInt k = 0; k < 2; k++) { + v[k][i] = v1[k]; + } + } // End of Quadrature Point Loop + + return 0; +} From 3bdd4e5a84903d9ef10177c79e65525d59684ccd Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Mon, 17 Apr 2023 15:44:01 -0700 Subject: [PATCH 3/6] clang-format --- backends/ref/ceed-ref-restriction.c | 46 ++++++++++++++++++----------- interface/ceed-preconditioning.c | 3 +- tests/t570-operator.c | 3 +- tests/t570-operator.h | 4 +-- 4 files changed, 33 insertions(+), 23 deletions(-) diff --git a/backends/ref/ceed-ref-restriction.c b/backends/ref/ceed-ref-restriction.c index 3f0c282649..5e1a33ed8a 100644 --- a/backends/ref/ceed-ref-restriction.c +++ b/backends/ref/ceed-ref-restriction.c @@ -16,8 +16,8 @@ // Core ElemRestriction Apply Code //------------------------------------------------------------------------------ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, - CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, + CeedVector v, CeedRequest *request) { CeedElemRestriction_Ref *impl; CeedCallBackend(CeedElemRestrictionGetData(r, &impl)); const CeedScalar *uu; @@ -81,8 +81,7 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, const CeedPragmaSIMD for (CeedInt i = 0; i < elem_size * blk_size; i++) { if (!use_orient || !impl->orient) { // Unsigned restriction - vv[elem_size * (k * blk_size + num_comp * e) + i - v_offset] = - uu[impl->offsets[i + elem_size * e] + k * comp_stride]; + vv[elem_size * (k * blk_size + num_comp * e) + i - v_offset] = uu[impl->offsets[i + elem_size * e] + k * comp_stride]; } else { // Signed restriction vv[elem_size * (k * blk_size + num_comp * e) + i - v_offset] = @@ -138,8 +137,7 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, const for (CeedInt j = i; j < i + CeedIntMin(blk_size, num_elem - e); j++) { if (!use_orient || !impl->orient) { // Unsigned restriction - vv[impl->offsets[j + e * elem_size] + k * comp_stride] += - uu[elem_size * (k * blk_size + num_comp * e) + j - v_offset]; + vv[impl->offsets[j + e * elem_size] + k * comp_stride] += uu[elem_size * (k * blk_size + num_comp * e) + j - v_offset]; } else { // Signed restriction vv[impl->offsets[j + e * elem_size] + k * comp_stride] += @@ -161,66 +159,78 @@ static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, const // ElemRestriction Apply - Common Sizes //------------------------------------------------------------------------------ static int CeedElemRestrictionApply_Ref_110(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 1, 1, comp_stride, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_111(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 1, 1, 1, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_180(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 1, 8, comp_stride, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_181(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 1, 8, 1, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_310(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 3, 1, comp_stride, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_311(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 3, 1, 1, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_380(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 3, 8, comp_stride, start, stop, use_orient, t_mode, u, v, request); } static int CeedElemRestrictionApply_Ref_381(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 3, 8, 1, start, stop, use_orient, t_mode, u, v, request); } // LCOV_EXCL_START static int CeedElemRestrictionApply_Ref_510(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 5, 1, comp_stride, start, stop, use_orient, t_mode, u, v, request); } // LCOV_EXCL_STOP static int CeedElemRestrictionApply_Ref_511(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 5, 1, 1, start, stop, use_orient, t_mode, u, v, request); } // LCOV_EXCL_START static int CeedElemRestrictionApply_Ref_580(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 5, 8, comp_stride, start, stop, use_orient, t_mode, u, v, request); } // LCOV_EXCL_STOP static int CeedElemRestrictionApply_Ref_581(CeedElemRestriction r, const CeedInt num_comp, const CeedInt blk_size, const CeedInt comp_stride, - CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { + CeedInt start, CeedInt stop, bool use_orient, CeedTransposeMode t_mode, CeedVector u, CeedVector v, + CeedRequest *request) { return CeedElemRestrictionApply_Ref_Core(r, 5, 8, 1, start, stop, use_orient, t_mode, u, v, request); } diff --git a/interface/ceed-preconditioning.c b/interface/ceed-preconditioning.c index ded55ddca5..edf41d67bb 100644 --- a/interface/ceed-preconditioning.c +++ b/interface/ceed-preconditioning.c @@ -442,7 +442,8 @@ static inline int CeedCompositeOperatorLinearAssembleAddDiagonal(CeedOperator op Users should generally use CeedOperatorLinearAssembleSymbolic() - Note: For operators using oriented element restrictions, entries in rows or cols may be negative indicating the assembled value at this nonzero should be negated + Note: For operators using oriented element restrictions, entries in rows or cols may be negative indicating the assembled value at this nonzero +should be negated @param[in] op CeedOperator to assemble nonzero pattern @param[in] offset Offset for number of entries diff --git a/tests/t570-operator.c b/tests/t570-operator.c index 20eb3b4336..3da2208be3 100644 --- a/tests/t570-operator.c +++ b/tests/t570-operator.c @@ -1,13 +1,14 @@ /// @file /// Test assembly of H(div) mass matrix operator diagonal /// \test Test assembly of H(div) mass matrix operator diagonal +#include "t570-operator.h" + #include #include #include #include #include "t330-basis.h" -#include "t570-operator.h" int main(int argc, char **argv) { Ceed ceed; diff --git a/tests/t570-operator.h b/tests/t570-operator.h index cab7f4865b..e5cc2eb467 100644 --- a/tests/t570-operator.h +++ b/tests/t570-operator.h @@ -8,9 +8,7 @@ #include // Compute det(A) -CEED_QFUNCTION_HELPER CeedScalar MatDet2x2(const CeedScalar A[2][2]) { - return A[0][0] * A[1][1] - A[1][0] * A[0][1]; -} +CEED_QFUNCTION_HELPER CeedScalar MatDet2x2(const CeedScalar A[2][2]) { return A[0][0] * A[1][1] - A[1][0] * A[0][1]; } // Compute alpha * A^T * B = C CEED_QFUNCTION_HELPER int AlphaMatTransposeMatMult2x2(const CeedScalar alpha, const CeedScalar A[2][2], const CeedScalar B[2][2], From c7be5f818cd9fc8bddb81c9f9576a7155fc22c9f Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Mon, 24 Apr 2023 14:42:23 -0700 Subject: [PATCH 4/6] Rebase fix for new CeedCheck macro --- interface/ceed-elemrestriction.c | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/interface/ceed-elemrestriction.c b/interface/ceed-elemrestriction.c index 11bf0c34c1..88dc307a97 100644 --- a/interface/ceed-elemrestriction.c +++ b/interface/ceed-elemrestriction.c @@ -656,20 +656,10 @@ int CeedElemRestrictionApplyUnsigned(CeedElemRestriction rstr, CeedTransposeMode m = rstr->l_size; n = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; } - if (n != u->length) { - // LCOV_EXCL_START - return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, - "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", u->length, m, - n); - // LCOV_EXCL_STOP - } - if (m != ru->length) { - // LCOV_EXCL_START - return CeedError(rstr->ceed, CEED_ERROR_DIMENSION, - "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", ru->length, - m, n); - // LCOV_EXCL_STOP - } + CeedCheck(n == u->length, rstr->ceed, CEED_ERROR_DIMENSION, + "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", u->length, m, n); + CeedCheck(m == ru->length, rstr->ceed, CEED_ERROR_DIMENSION, + "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", ru->length, m, n); if (rstr->num_elem > 0) { if (rstr->ApplyUnsigned) CeedCall(rstr->ApplyUnsigned(rstr, t_mode, u, ru, request)); else CeedCall(rstr->Apply(rstr, t_mode, u, ru, request)); From b17517ee98e70bbde11768f0b2de0644e9c68f9a Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Tue, 25 Apr 2023 10:18:21 -0700 Subject: [PATCH 5/6] Remove default fallback from CeedElemRestrictionApplyUnsigned to CeedElemRestrictionApply and use backend registration instead --- backends/cuda-ref/ceed-cuda-restriction.c | 1 + backends/hip-ref/ceed-hip-ref-restriction.c | 1 + backends/magma/ceed-magma-restriction.c | 1 + backends/occa/ceed-occa-elem-restriction.cpp | 1 + interface/ceed-elemrestriction.c | 7 +++---- 5 files changed, 7 insertions(+), 4 deletions(-) diff --git a/backends/cuda-ref/ceed-cuda-restriction.c b/backends/cuda-ref/ceed-cuda-restriction.c index 6104b10be4..4126ae4245 100644 --- a/backends/cuda-ref/ceed-cuda-restriction.c +++ b/backends/cuda-ref/ceed-cuda-restriction.c @@ -322,6 +322,7 @@ int CeedElemRestrictionCreate_Cuda(CeedMemType mem_type, CeedCopyMode copy_mode, // Register backend functions CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", CeedElemRestrictionApply_Cuda)); + CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyUnsigned", CeedElemRestrictionApply_Cuda)); CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "GetOffsets", CeedElemRestrictionGetOffsets_Cuda)); CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", CeedElemRestrictionDestroy_Cuda)); return CEED_ERROR_SUCCESS; diff --git a/backends/hip-ref/ceed-hip-ref-restriction.c b/backends/hip-ref/ceed-hip-ref-restriction.c index cf7f05d6a3..dee2082309 100644 --- a/backends/hip-ref/ceed-hip-ref-restriction.c +++ b/backends/hip-ref/ceed-hip-ref-restriction.c @@ -320,6 +320,7 @@ int CeedElemRestrictionCreate_Hip(CeedMemType mem_type, CeedCopyMode copy_mode, // Register backend functions CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", CeedElemRestrictionApply_Hip)); + CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyUnsigned", CeedElemRestrictionApply_Hip)); CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "GetOffsets", CeedElemRestrictionGetOffsets_Hip)); CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", CeedElemRestrictionDestroy_Hip)); return CEED_ERROR_SUCCESS; diff --git a/backends/magma/ceed-magma-restriction.c b/backends/magma/ceed-magma-restriction.c index bdd68edeb6..4f5a92fa06 100644 --- a/backends/magma/ceed-magma-restriction.c +++ b/backends/magma/ceed-magma-restriction.c @@ -260,6 +260,7 @@ int CeedElemRestrictionCreate_Magma(CeedMemType mtype, CeedCopyMode cmode, const CeedInt layout[3] = {1, elemsize * nelem, elemsize}; CeedCallBackend(CeedElemRestrictionSetELayout(r, layout)); CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", CeedElemRestrictionApply_Magma)); + CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyUnsigned", CeedElemRestrictionApply_Magma)); CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "GetOffsets", CeedElemRestrictionGetOffsets_Magma)); CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", CeedElemRestrictionDestroy_Magma)); CeedCallBackend(CeedFree(&restriction_kernel_path)); diff --git a/backends/occa/ceed-occa-elem-restriction.cpp b/backends/occa/ceed-occa-elem-restriction.cpp index 11983a8a72..b7580d2b36 100644 --- a/backends/occa/ceed-occa-elem-restriction.cpp +++ b/backends/occa/ceed-occa-elem-restriction.cpp @@ -316,6 +316,7 @@ int ElemRestriction::ceedCreate(CeedMemType memType, CeedCopyMode copyMode, cons CeedChkBackend(CeedElemRestrictionSetELayout(r, defaultLayout)); CeedOccaRegisterFunction(r, "Apply", ElemRestriction::ceedApply); + CeedOccaRegisterFunction(r, "ApplyUnsigned", ElemRestriction::ceedApply); CeedOccaRegisterFunction(r, "ApplyBlock", ElemRestriction::ceedApplyBlock); CeedOccaRegisterFunction(r, "GetOffsets", ElemRestriction::ceedGetOffsets); CeedOccaRegisterFunction(r, "Destroy", ElemRestriction::ceedDestroy); diff --git a/interface/ceed-elemrestriction.c b/interface/ceed-elemrestriction.c index 88dc307a97..a0be65f9ad 100644 --- a/interface/ceed-elemrestriction.c +++ b/interface/ceed-elemrestriction.c @@ -649,6 +649,8 @@ int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, int CeedElemRestrictionApplyUnsigned(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) { CeedInt m, n; + CeedCheck(rstr->ApplyUnsigned, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionApplyUnsigned"); + if (t_mode == CEED_NOTRANSPOSE) { m = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp; n = rstr->l_size; @@ -660,10 +662,7 @@ int CeedElemRestrictionApplyUnsigned(CeedElemRestriction rstr, CeedTransposeMode "Input vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", u->length, m, n); CeedCheck(m == ru->length, rstr->ceed, CEED_ERROR_DIMENSION, "Output vector size %" CeedInt_FMT " not compatible with element restriction (%" CeedInt_FMT ", %" CeedInt_FMT ")", ru->length, m, n); - if (rstr->num_elem > 0) { - if (rstr->ApplyUnsigned) CeedCall(rstr->ApplyUnsigned(rstr, t_mode, u, ru, request)); - else CeedCall(rstr->Apply(rstr, t_mode, u, ru, request)); - } + if (rstr->num_elem > 0) CeedCall(rstr->ApplyUnsigned(rstr, t_mode, u, ru, request)); return CEED_ERROR_SUCCESS; } From b04eb3d9dcd03054fbe981a83f15e207a4bc97bc Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Tue, 25 Apr 2023 11:18:29 -0700 Subject: [PATCH 6/6] Consistency for 'does not support' --- interface/ceed-elemrestriction.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/interface/ceed-elemrestriction.c b/interface/ceed-elemrestriction.c index a0be65f9ad..1fc10b115d 100644 --- a/interface/ceed-elemrestriction.c +++ b/interface/ceed-elemrestriction.c @@ -409,6 +409,7 @@ int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt num_elem, CeedInt elem_s CeedElemRestriction *rstr) { if (!ceed->ElemRestrictionCreate) { Ceed delegate; + CeedCall(CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction")); CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support ElemRestrictionCreateStrided"); CeedCall(CeedElemRestrictionCreateStrided(delegate, num_elem, elem_size, num_comp, l_size, strides, rstr)); @@ -649,7 +650,7 @@ int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode t_mode, int CeedElemRestrictionApplyUnsigned(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request) { CeedInt m, n; - CeedCheck(rstr->ApplyUnsigned, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement ElemRestrictionApplyUnsigned"); + CeedCheck(rstr->ApplyUnsigned, rstr->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support ElemRestrictionApplyUnsigned"); if (t_mode == CEED_NOTRANSPOSE) { m = rstr->num_blk * rstr->blk_size * rstr->elem_size * rstr->num_comp;