From 9f46b6c6e7587225052a2699c244c4dd386fc0bf Mon Sep 17 00:00:00 2001 From: alvaroborras <13698600+alvaroborras@users.noreply.github.com> Date: Fri, 28 Aug 2026 17:20:52 +0200 Subject: [PATCH 1/4] Add FLOP estimates for operator assembly Signed-off-by: alvaroborras <13698600+alvaroborras@users.noreply.github.com> --- include/ceed/ceed.h | 3 + interface/ceed-preconditioning.c | 161 +++++++++++++++++++++++++++++++ tests/t526-operator.c | 34 ++++++- tests/t537-operator.c | 22 ++++- tests/t570-operator.c | 32 +++++- tests/t595-operator.c | 6 ++ 6 files changed, 248 insertions(+), 10 deletions(-) diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h index f72ea9c69b..615a0a5e3a 100644 --- a/include/ceed/ceed.h +++ b/include/ceed/ceed.h @@ -457,13 +457,16 @@ CEED_EXTERN int CeedOperatorSetQFunctionAssemblyDataUpdateNeeded(CeedOperator o CEED_EXTERN int CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request); CEED_EXTERN int CeedOperatorLinearAssembleQFunctionBuildOrUpdate(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request); +CEED_EXTERN int CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(CeedOperator op, CeedSize *flops); CEED_EXTERN int CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request); CEED_EXTERN int CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request); +CEED_EXTERN int CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(CeedOperator op, CeedSize *flops); CEED_EXTERN int CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request); CEED_EXTERN int CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request); CEED_EXTERN int CeedOperatorLinearAssemblePointBlockDiagonalSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols); CEED_EXTERN int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols); CEED_EXTERN int CeedOperatorLinearAssembleGetNumEntries(CeedOperator op, CeedSize *num_entries); +CEED_EXTERN int CeedOperatorLinearAssembleGetFlopsEstimate(CeedOperator op, CeedSize *flops); CEED_EXTERN int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values); CEED_EXTERN int CeedOperatorCompositeGetMultiplicity(CeedOperator op, CeedInt num_skip_indices, CeedInt *skip_indices, CeedVector mult); CEED_EXTERN int CeedOperatorMultigridLevelCreate(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, diff --git a/interface/ceed-preconditioning.c b/interface/ceed-preconditioning.c index d3376aec0f..08f668eba5 100644 --- a/interface/ceed-preconditioning.c +++ b/interface/ceed-preconditioning.c @@ -2375,6 +2375,167 @@ int CeedOperatorLinearAssembleQFunctionBuildOrUpdate(CeedOperator op, CeedVector return CeedOperatorLinearAssembleQFunctionBuildOrUpdate_Core(op, *assembled == NULL, true, assembled, rstr, request); } +// Estimate the work to rebuild the assembled QFunction data, without changing its state. +static int CeedOperatorLinearAssembleQFunctionGetFlopsEstimate(CeedOperator op, CeedInt num_qpts, CeedSize *flops) { + CeedInt num_elem, num_input_fields, num_active_inputs = 0; + CeedQFunction qf; + CeedQFunctionField *qf_input_fields; + CeedOperatorField *op_input_fields; + CeedQFunctionAssemblyData data; + bool is_setup, update_needed; + + *flops = 0; + CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data)); + CeedCall(CeedQFunctionAssemblyDataIsSetup(data, &is_setup)); + if (is_setup) { + CeedCall(CeedQFunctionAssemblyDataIsUpdateNeeded(data, &update_needed)); + } else { + update_needed = true; + } + if (!update_needed) return CEED_ERROR_SUCCESS; + + CeedCall(CeedOperatorGetNumElements(op, &num_elem)); + CeedCall(CeedOperatorGetQFunction(op, &qf)); + CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_input_fields, NULL, NULL)); + CeedCall(CeedOperatorGetFields(op, NULL, &op_input_fields, NULL, NULL)); + for (CeedInt i = 0; i < num_input_fields; i++) { + CeedEvalMode eval_mode; + CeedVector vec; + + CeedCall(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); + if (vec == CEED_VECTOR_ACTIVE) { + CeedInt size; + + CeedCall(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); + num_active_inputs += size; + } else { + CeedBasis basis; + CeedElemRestriction rstr; + CeedSize basis_flops, rstr_flops; + + CeedCall(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); + if (eval_mode != CEED_EVAL_WEIGHT) { + CeedCall(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr)); + CeedCall(CeedElemRestrictionGetFlopsEstimate(rstr, CEED_NOTRANSPOSE, &rstr_flops)); + CeedCall(CeedElemRestrictionDestroy(&rstr)); + CeedCall(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); + basis_flops = 0; + if (basis != CEED_BASIS_NONE) CeedCall(CeedBasisGetFlopsEstimate(basis, CEED_NOTRANSPOSE, eval_mode, false, 0, &basis_flops)); + CeedCall(CeedBasisDestroy(&basis)); + *flops += rstr_flops + basis_flops * num_elem; + } + } + CeedCall(CeedVectorDestroy(&vec)); + } + { + CeedSize qf_flops; + + CeedCall(CeedQFunctionGetFlopsEstimate(qf, &qf_flops)); + CeedCheck(qf_flops > -1, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, + "Must set CeedQFunction FLOPs estimate with CeedQFunctionSetUserFlopsEstimate"); + *flops += (CeedSize)num_elem * num_qpts * num_active_inputs * qf_flops; + } + CeedCall(CeedQFunctionDestroy(&qf)); + return CEED_ERROR_SUCCESS; +} + +static int CeedOperatorLinearAssembleGetFlopsEstimate_Core(CeedOperator op, bool is_point_block, bool is_diagonal, CeedSize *flops) { + bool is_at_points, is_composite; + CeedInt num_points; + + CeedCall(CeedOperatorCheckReady(op)); + *flops = 0; + CeedCall(CeedOperatorIsComposite(op, &is_composite)); + if (is_composite) { + CeedInt num_sub; + CeedOperator *suboperators; + + CeedCall(CeedOperatorCompositeGetNumSub(op, &num_sub)); + CeedCall(CeedOperatorCompositeGetSubList(op, &suboperators)); + for (CeedInt i = 0; i < num_sub; i++) { + CeedSize sub_flops; + + CeedCall(CeedOperatorLinearAssembleGetFlopsEstimate_Core(suboperators[i], is_point_block, is_diagonal, &sub_flops)); + *flops += sub_flops; + } + return CEED_ERROR_SUCCESS; + } + CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); + if (is_at_points) { + CeedElemRestriction rstr_points; + + CeedCall(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); + CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &num_points)); + CeedCall(CeedElemRestrictionDestroy(&rstr_points)); + } else { + CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_points)); + } + + CeedCall(CeedOperatorLinearAssembleQFunctionGetFlopsEstimate(op, num_points, flops)); + { + CeedOperatorAssemblyData data; + CeedBasis *bases_in, *bases_out; + CeedElemRestriction *rstrs_in, *rstrs_out; + CeedInt num_bases_in, num_bases_out, *num_eval_modes_in, *num_eval_modes_out; + + CeedCall(CeedOperatorGetOperatorAssemblyData(op, &data)); + CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_bases_in, &num_eval_modes_in, NULL, NULL, &num_bases_out, &num_eval_modes_out, + NULL, NULL, NULL)); + CeedCall(CeedOperatorAssemblyDataGetBases(data, NULL, &bases_in, NULL, NULL, &bases_out, NULL)); + CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, NULL, &rstrs_in, NULL, &rstrs_out)); + for (CeedInt b_in = 0; b_in < num_bases_in; b_in++) { + for (CeedInt b_out = 0; b_out < num_bases_out; b_out++) { + CeedInt num_elem, num_nodes_in, num_nodes_out, num_comp_in, num_comp_out, num_qpts; + + if (is_diagonal && bases_in[b_in] != bases_out[b_out]) continue; + CeedCall(CeedElemRestrictionGetNumElements(rstrs_in[b_in], &num_elem)); + CeedCall(CeedElemRestrictionGetNumComponents(rstrs_in[b_in], &num_comp_in)); + CeedCall(CeedElemRestrictionGetNumComponents(rstrs_out[b_out], &num_comp_out)); + if (bases_in[b_in] == CEED_BASIS_NONE) { + CeedCall(CeedElemRestrictionGetElementSize(rstrs_in[b_in], &num_nodes_in)); + } else { + CeedCall(CeedBasisGetNumNodes(bases_in[b_in], &num_nodes_in)); + } + if (is_at_points) { + num_qpts = num_points; + } else if (bases_in[b_in] == CEED_BASIS_NONE) { + num_qpts = num_nodes_in; + } else { + CeedCall(CeedBasisGetNumQuadraturePoints(bases_in[b_in], &num_qpts)); + } + if (bases_out[b_out] == CEED_BASIS_NONE) { + CeedCall(CeedElemRestrictionGetElementSize(rstrs_out[b_out], &num_nodes_out)); + } else { + CeedCall(CeedBasisGetNumNodes(bases_out[b_out], &num_nodes_out)); + } + if (is_diagonal) { + const CeedSize num_comp = is_point_block ? num_comp_in * num_comp_out : num_comp_in; + + *flops += 3 * (CeedSize)num_elem * num_nodes_in * num_qpts * num_eval_modes_in[b_in] * num_eval_modes_out[b_out] * num_comp; + *flops += (CeedSize)num_elem * num_nodes_in * num_comp; + if (is_point_block) *flops += (CeedSize)num_elem * num_nodes_in * num_comp * (num_comp_out - 1); + } else { + *flops += 2 * (CeedSize)num_elem * num_comp_in * num_comp_out * num_nodes_out * num_qpts * num_eval_modes_in[b_in] * + (num_eval_modes_out[b_out] + num_nodes_in); + } + } + } + } + return CEED_ERROR_SUCCESS; +} + +int CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(CeedOperator op, CeedSize *flops) { + return CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, false, true, flops); +} + +int CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(CeedOperator op, CeedSize *flops) { + return CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, true, true, flops); +} + +int CeedOperatorLinearAssembleGetFlopsEstimate(CeedOperator op, CeedSize *flops) { + return CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, false, false, flops); +} + /** @brief Assemble the diagonal of a square linear `CeedOperator` diff --git a/tests/t526-operator.c b/tests/t526-operator.c index 8e68ab89b3..e383769d3f 100644 --- a/tests/t526-operator.c +++ b/tests/t526-operator.c @@ -1,6 +1,6 @@ /// @file -/// Test FLOP estimation for composite mass matrix operator -/// \test Test FLOP estimation for composite mass matrix operator +/// Test apply and assembly FLOP estimation for composite mass matrix operator +/// \test Test apply and assembly FLOP estimation for composite mass matrix operator #include #include #include @@ -19,13 +19,13 @@ int main(int argc, char **argv) { Ceed ceed; - CeedSize flop_estimate; + CeedSize flop_estimate, num_entries; CeedElemRestriction elem_restriction_x_tet, elem_restriction_u_tet, elem_restriction_q_data_tet, elem_restriction_x_hex, elem_restriction_u_hex, elem_restriction_q_data_hex; CeedBasis basis_x_tet, basis_u_tet, basis_x_hex, basis_u_hex; CeedQFunction qf_mass; CeedOperator op_mass_tet, op_mass_hex, op_mass; - CeedVector q_data_tet, q_data_hex; + CeedVector q_data_tet, q_data_hex, assembled; CeedInt num_elem_tet = 6, p_tet = 6, q_tet = 4, num_elem_hex = 6, p_hex = 3, q_hex = 4, dim = 2; CeedInt n_x = 3, n_y = 3, n_x_tet = 3, n_y_tet = 1, n_x_hex = 3; CeedInt row, col, offset; @@ -39,6 +39,8 @@ int main(int argc, char **argv) { // Qdata Vectors CeedVectorCreate(ceed, num_qpts_tet, &q_data_tet); CeedVectorCreate(ceed, num_qpts_hex, &q_data_hex); + CeedVectorSetValue(q_data_tet, 1.0); + CeedVectorSetValue(q_data_hex, 1.0); // Set up Tet Elements // -- Restrictions @@ -123,12 +125,34 @@ int main(int argc, char **argv) { CeedQFunctionSetUserFlopsEstimate(qf_mass, 1); CeedOperatorGetFlopsEstimate(op_mass, &flop_estimate); - // Check output + // Check apply estimate if (flop_estimate != 3042) printf("Incorrect FLOP estimate computed, %" CeedSize_FMT " != 3042\n", flop_estimate); + // Check assembly estimates with stale QFunction data. Repeating the query must not update it. + CeedOperatorLinearAssembleGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 19416) printf("Incorrect full assembly FLOP estimate, %" CeedSize_FMT " != 19416\n", flop_estimate); + CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 3234) printf("Incorrect diagonal assembly FLOP estimate, %" CeedSize_FMT " != 3234\n", flop_estimate); + CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 3234) printf("Incorrect point-block diagonal assembly FLOP estimate, %" CeedSize_FMT " != 3234\n", flop_estimate); + + // Assemble once with reuse enabled, then check the cached-QFunction path + CeedOperatorSetQFunctionAssemblyReuse(op_mass, true); + CeedOperatorLinearAssembleGetNumEntries(op_mass, &num_entries); + CeedVectorCreate(ceed, num_entries, &assembled); + CeedVectorSetValue(assembled, 0.0); + CeedOperatorLinearAssemble(op_mass, assembled); + CeedOperatorLinearAssembleGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 19296) printf("Incorrect cached full assembly FLOP estimate, %" CeedSize_FMT " != 19296\n", flop_estimate); + CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 3114) printf("Incorrect cached diagonal assembly FLOP estimate, %" CeedSize_FMT " != 3114\n", flop_estimate); + CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 3114) printf("Incorrect cached point-block diagonal assembly FLOP estimate, %" CeedSize_FMT " != 3114\n", flop_estimate); + // Cleanup CeedVectorDestroy(&q_data_tet); CeedVectorDestroy(&q_data_hex); + CeedVectorDestroy(&assembled); CeedElemRestrictionDestroy(&elem_restriction_u_tet); CeedElemRestrictionDestroy(&elem_restriction_x_tet); CeedElemRestrictionDestroy(&elem_restriction_q_data_tet); diff --git a/tests/t537-operator.c b/tests/t537-operator.c index d7e0ffcc4b..63d9bbcbdb 100644 --- a/tests/t537-operator.c +++ b/tests/t537-operator.c @@ -1,6 +1,6 @@ /// @file -/// Test assembly of mass matrix operator point block diagonal -/// \test Test assembly of mass matrix operator point block diagonal +/// Test point-block-diagonal assembly and FLOP estimation of mass matrix operator +/// \test Test point-block-diagonal assembly and FLOP estimation of mass matrix operator #include "t537-operator.h" #include @@ -20,7 +20,7 @@ int main(int argc, char **argv) { CeedInt num_dofs = (n_x * 2 + 1) * (n_y * 2 + 1), num_qpts = num_elem * q * q; CeedInt ind_x[num_elem * p * p]; CeedInt *rows, *cols; - CeedSize num_entries; + CeedSize num_entries, flop_estimate; CeedScalar assembled_true[num_comp * num_comp * num_dofs]; CeedScalar assembled_full_true[num_comp * num_dofs][num_comp * num_dofs]; @@ -87,10 +87,26 @@ int main(int argc, char **argv) { // Apply Setup Operator CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE); + // Estimate assembly FLOPs before assembling QFunction data + CeedQFunctionSetUserFlopsEstimate(qf_mass, 1); + CeedOperatorLinearAssembleGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 69312) printf("Incorrect full assembly FLOP estimate, %" CeedSize_FMT " != 69312\n", flop_estimate); + CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 5484) printf("Incorrect diagonal assembly FLOP estimate, %" CeedSize_FMT " != 5484\n", flop_estimate); + CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 10992) printf("Incorrect point-block diagonal assembly FLOP estimate, %" CeedSize_FMT " != 10992\n", flop_estimate); + // Assemble diagonal + CeedOperatorSetQFunctionAssemblyReuse(op_mass, true); CeedVectorCreate(ceed, num_comp * num_comp * num_dofs, &assembled); CeedOperatorLinearAssemblePointBlockDiagonal(op_mass, assembled, CEED_REQUEST_IMMEDIATE); CeedOperatorLinearAssemblePointBlockDiagonalSymbolic(op_mass, &num_entries, &rows, &cols); + CeedOperatorLinearAssembleGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 69120) printf("Incorrect cached full assembly FLOP estimate, %" CeedSize_FMT " != 69120\n", flop_estimate); + CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 5292) printf("Incorrect cached diagonal assembly FLOP estimate, %" CeedSize_FMT " != 5292\n", flop_estimate); + CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 10800) printf("Incorrect cached point-block diagonal assembly FLOP estimate, %" CeedSize_FMT " != 10800\n", flop_estimate); // Manually assemble point-block diagonal CeedVectorCreate(ceed, num_comp * num_dofs, &u); diff --git a/tests/t570-operator.c b/tests/t570-operator.c index c627c0f5bb..213c949e02 100644 --- a/tests/t570-operator.c +++ b/tests/t570-operator.c @@ -1,6 +1,6 @@ /// @file -/// Test full assembly of an identity operator (see t509) -/// \test Test full assembly of an identity operator +/// Test full assembly and FLOP estimation of an identity operator (see t509) +/// \test Test full assembly and FLOP estimation of an identity operator #include #include #include @@ -46,6 +46,34 @@ int main(int argc, char **argv) { CeedOperatorSetField(op_identity, "input", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); CeedOperatorSetField(op_identity, "output", elem_restriction_u_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); + // Estimate assembly FLOPs with a basis-none active output + { + CeedSize flop_estimate; + + CeedQFunctionSetUserFlopsEstimate(qf_identity, 0); + CeedOperatorLinearAssembleGetFlopsEstimate(op_identity, &flop_estimate); + if (flop_estimate != 11520) printf("Incorrect full assembly FLOP estimate, %" CeedSize_FMT " != 11520\n", flop_estimate); + } + { + CeedOperator op_identity_none; + CeedQFunction qf_identity_none; + CeedSize flop_estimate; + + CeedQFunctionCreateIdentity(ceed, 1, CEED_EVAL_NONE, CEED_EVAL_NONE, &qf_identity_none); + CeedQFunctionSetUserFlopsEstimate(qf_identity_none, 0); + CeedOperatorCreate(ceed, qf_identity_none, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_identity_none); + CeedOperatorSetField(op_identity_none, "input", elem_restriction_u_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_identity_none, "output", elem_restriction_u_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); + CeedOperatorLinearAssembleGetFlopsEstimate(op_identity_none, &flop_estimate); + if (flop_estimate != 17280) printf("Incorrect basis-none full assembly FLOP estimate, %" CeedSize_FMT " != 17280\n", flop_estimate); + CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(op_identity_none, &flop_estimate); + if (flop_estimate != 3000) printf("Incorrect basis-none diagonal assembly FLOP estimate, %" CeedSize_FMT " != 3000\n", flop_estimate); + CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(op_identity_none, &flop_estimate); + if (flop_estimate != 3000) printf("Incorrect basis-none point-block diagonal assembly FLOP estimate, %" CeedSize_FMT " != 3000\n", flop_estimate); + CeedOperatorDestroy(&op_identity_none); + CeedQFunctionDestroy(&qf_identity_none); + } + // Fully assemble operator CeedSize num_entries; CeedInt *rows; diff --git a/tests/t595-operator.c b/tests/t595-operator.c index e874ccb2ba..cb620e1bed 100644 --- a/tests/t595-operator.c +++ b/tests/t595-operator.c @@ -110,6 +110,12 @@ int main(int argc, char **argv) { printf("Incorrect FLOP estimate computed, %ld != 16317\n", flop_estimate); // LCOV_EXCL_STOP } + CeedOperatorLinearAssembleGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 6516) printf("Incorrect AtPoints full assembly FLOP estimate, %" CeedSize_FMT " != 6516\n", flop_estimate); + CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 1089) printf("Incorrect AtPoints diagonal assembly FLOP estimate, %" CeedSize_FMT " != 1089\n", flop_estimate); + CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(op_mass, &flop_estimate); + if (flop_estimate != 1089) printf("Incorrect AtPoints point-block diagonal assembly FLOP estimate, %" CeedSize_FMT " != 1089\n", flop_estimate); CeedVectorDestroy(&x_points); CeedVectorDestroy(&q_data); From 0239b02c11349715c3120ca34a23892b0c9efb3a Mon Sep 17 00:00:00 2001 From: alvaroborras <13698600+alvaroborras@users.noreply.github.com> Date: Tue, 1 Sep 2026 08:22:38 +0200 Subject: [PATCH 2/4] Fix AtPoints assembly FLOP estimates Signed-off-by: alvaroborras <13698600+alvaroborras@users.noreply.github.com> --- CHANGELOG.md | 2 + interface/ceed-preconditioning.c | 124 ++++++++++++++++++++++++------- tests/t526-operator.c | 6 +- tests/t595-operator.c | 52 +++++++------ 4 files changed, 132 insertions(+), 52 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 62d79f46a2..36b6568c33 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,8 @@ On this page we provide a summary of the main API changes, new features and exam ### Interface changes +- Add `CeedOperatorLinearAssembleGetFlopsEstimate`, `CeedOperatorLinearAssembleDiagonalGetFlopsEstimate`, and + `CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate` to estimate FLOPs for linear operator assembly. - Add `bool` field type for `CeedQFunctionContext` and related interfaces to use `bool` fields. - `CEED_BASIS_COLLOCATED` removed; users should only use `CEED_BASIS_NONE`. - Remove unneeded pointer for `CeedElemRestrictionGetELayout`. diff --git a/interface/ceed-preconditioning.c b/interface/ceed-preconditioning.c index 08f668eba5..f0a5debb3c 100644 --- a/interface/ceed-preconditioning.c +++ b/interface/ceed-preconditioning.c @@ -2376,13 +2376,13 @@ int CeedOperatorLinearAssembleQFunctionBuildOrUpdate(CeedOperator op, CeedVector } // Estimate the work to rebuild the assembled QFunction data, without changing its state. -static int CeedOperatorLinearAssembleQFunctionGetFlopsEstimate(CeedOperator op, CeedInt num_qpts, CeedSize *flops) { - CeedInt num_elem, num_input_fields, num_active_inputs = 0; - CeedQFunction qf; - CeedQFunctionField *qf_input_fields; - CeedOperatorField *op_input_fields; - CeedQFunctionAssemblyData data; - bool is_setup, update_needed; +static int CeedOperatorLinearAssembleQFunctionGetFlopsEstimate(CeedOperator op, CeedSize num_qpts, CeedSize *flops) { + CeedInt num_elem, num_input_fields, num_active_inputs = 0; + CeedQFunction qf; + CeedQFunctionField *qf_input_fields; + CeedOperatorField *op_input_fields; + CeedQFunctionAssemblyData data; + bool is_setup, update_needed; *flops = 0; CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data)); @@ -2433,15 +2433,15 @@ static int CeedOperatorLinearAssembleQFunctionGetFlopsEstimate(CeedOperator op, CeedCall(CeedQFunctionGetFlopsEstimate(qf, &qf_flops)); CeedCheck(qf_flops > -1, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPLETE, "Must set CeedQFunction FLOPs estimate with CeedQFunctionSetUserFlopsEstimate"); - *flops += (CeedSize)num_elem * num_qpts * num_active_inputs * qf_flops; + *flops += num_qpts * num_active_inputs * qf_flops; } CeedCall(CeedQFunctionDestroy(&qf)); return CEED_ERROR_SUCCESS; } static int CeedOperatorLinearAssembleGetFlopsEstimate_Core(CeedOperator op, bool is_point_block, bool is_diagonal, CeedSize *flops) { - bool is_at_points, is_composite; - CeedInt num_points; + bool is_at_points, is_composite; + CeedSize num_qpts_total; CeedCall(CeedOperatorCheckReady(op)); *flops = 0; @@ -2462,16 +2462,38 @@ static int CeedOperatorLinearAssembleGetFlopsEstimate_Core(CeedOperator op, bool } CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); if (is_at_points) { + CeedInt num_elem; + CeedMemType mem_type; CeedElemRestriction rstr_points; CeedCall(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); - CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &num_points)); + CeedCall(CeedOperatorGetNumElements(op, &num_elem)); + CeedCall(CeedGetPreferredMemType(CeedOperatorReturnCeed(op), &mem_type)); + if (mem_type == CEED_MEM_DEVICE) { + CeedInt max_points; + + // Device backends pad every element to the maximum number of points. + CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_points)); + num_qpts_total = (CeedSize)num_elem * max_points; + } else { + num_qpts_total = 0; + for (CeedInt i = 0; i < num_elem; i++) { + CeedInt points_in_elem; + + CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr_points, i, &points_in_elem)); + num_qpts_total += points_in_elem; + } + } CeedCall(CeedElemRestrictionDestroy(&rstr_points)); } else { - CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_points)); + CeedInt num_elem, num_qpts_per_elem; + + CeedCall(CeedOperatorGetNumElements(op, &num_elem)); + CeedCall(CeedOperatorGetNumQuadraturePoints(op, &num_qpts_per_elem)); + num_qpts_total = (CeedSize)num_elem * num_qpts_per_elem; } - CeedCall(CeedOperatorLinearAssembleQFunctionGetFlopsEstimate(op, num_points, flops)); + CeedCall(CeedOperatorLinearAssembleQFunctionGetFlopsEstimate(op, num_qpts_total, flops)); { CeedOperatorAssemblyData data; CeedBasis *bases_in, *bases_out; @@ -2479,16 +2501,16 @@ static int CeedOperatorLinearAssembleGetFlopsEstimate_Core(CeedOperator op, bool CeedInt num_bases_in, num_bases_out, *num_eval_modes_in, *num_eval_modes_out; CeedCall(CeedOperatorGetOperatorAssemblyData(op, &data)); - CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_bases_in, &num_eval_modes_in, NULL, NULL, &num_bases_out, &num_eval_modes_out, - NULL, NULL, NULL)); + CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_bases_in, &num_eval_modes_in, NULL, NULL, &num_bases_out, &num_eval_modes_out, NULL, + NULL, NULL)); CeedCall(CeedOperatorAssemblyDataGetBases(data, NULL, &bases_in, NULL, NULL, &bases_out, NULL)); CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, NULL, &rstrs_in, NULL, &rstrs_out)); for (CeedInt b_in = 0; b_in < num_bases_in; b_in++) { for (CeedInt b_out = 0; b_out < num_bases_out; b_out++) { - CeedInt num_elem, num_nodes_in, num_nodes_out, num_comp_in, num_comp_out, num_qpts; + CeedInt num_nodes_in, num_nodes_out, num_comp_in, num_comp_out; + CeedSize num_qpts = num_qpts_total; if (is_diagonal && bases_in[b_in] != bases_out[b_out]) continue; - CeedCall(CeedElemRestrictionGetNumElements(rstrs_in[b_in], &num_elem)); CeedCall(CeedElemRestrictionGetNumComponents(rstrs_in[b_in], &num_comp_in)); CeedCall(CeedElemRestrictionGetNumComponents(rstrs_out[b_out], &num_comp_out)); if (bases_in[b_in] == CEED_BASIS_NONE) { @@ -2496,12 +2518,20 @@ static int CeedOperatorLinearAssembleGetFlopsEstimate_Core(CeedOperator op, bool } else { CeedCall(CeedBasisGetNumNodes(bases_in[b_in], &num_nodes_in)); } - if (is_at_points) { - num_qpts = num_points; - } else if (bases_in[b_in] == CEED_BASIS_NONE) { - num_qpts = num_nodes_in; + if (!is_at_points && bases_in[b_in] == CEED_BASIS_NONE) { + CeedInt num_elem, num_qpts_per_elem; + + CeedCall(CeedElemRestrictionGetNumElements(rstrs_in[b_in], &num_elem)); + num_qpts_per_elem = num_nodes_in; + num_qpts = (CeedSize)num_elem * num_qpts_per_elem; } else { - CeedCall(CeedBasisGetNumQuadraturePoints(bases_in[b_in], &num_qpts)); + if (!is_at_points) { + CeedInt num_elem, num_qpts_per_elem; + + CeedCall(CeedElemRestrictionGetNumElements(rstrs_in[b_in], &num_elem)); + CeedCall(CeedBasisGetNumQuadraturePoints(bases_in[b_in], &num_qpts_per_elem)); + num_qpts = (CeedSize)num_elem * num_qpts_per_elem; + } } if (bases_out[b_out] == CEED_BASIS_NONE) { CeedCall(CeedElemRestrictionGetElementSize(rstrs_out[b_out], &num_nodes_out)); @@ -2511,12 +2541,16 @@ static int CeedOperatorLinearAssembleGetFlopsEstimate_Core(CeedOperator op, bool if (is_diagonal) { const CeedSize num_comp = is_point_block ? num_comp_in * num_comp_out : num_comp_in; - *flops += 3 * (CeedSize)num_elem * num_nodes_in * num_qpts * num_eval_modes_in[b_in] * num_eval_modes_out[b_out] * num_comp; - *flops += (CeedSize)num_elem * num_nodes_in * num_comp; - if (is_point_block) *flops += (CeedSize)num_elem * num_nodes_in * num_comp * (num_comp_out - 1); + *flops += 3 * num_qpts * num_nodes_in * num_eval_modes_in[b_in] * num_eval_modes_out[b_out] * num_comp; + { + CeedInt num_elem; + + CeedCall(CeedElemRestrictionGetNumElements(rstrs_in[b_in], &num_elem)); + *flops += (CeedSize)num_elem * num_nodes_in * num_comp; + if (is_point_block) *flops += (CeedSize)num_elem * num_nodes_in * num_comp * (num_comp_out - 1); + } } else { - *flops += 2 * (CeedSize)num_elem * num_comp_in * num_comp_out * num_nodes_out * num_qpts * num_eval_modes_in[b_in] * - (num_eval_modes_out[b_out] + num_nodes_in); + *flops += 2 * num_comp_in * num_comp_out * num_nodes_out * num_qpts * num_eval_modes_in[b_in] * (num_eval_modes_out[b_out] + num_nodes_in); } } } @@ -2524,14 +2558,50 @@ static int CeedOperatorLinearAssembleGetFlopsEstimate_Core(CeedOperator op, bool return CEED_ERROR_SUCCESS; } +/** + @brief Estimate the number of FLOPs required to assemble the diagonal of a linear `CeedOperator`. + + This estimate accounts for rebuilding stale QFunction assembly data without modifying its state. + + @param[in] op `CeedOperator` to estimate FLOPs for + @param[out] flops Address of variable to hold FLOPs estimate + + @return An error code: 0 - success, otherwise - failure + + @ref User +**/ int CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(CeedOperator op, CeedSize *flops) { return CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, false, true, flops); } +/** + @brief Estimate the number of FLOPs required to assemble the point-block diagonal of a linear `CeedOperator`. + + This estimate accounts for rebuilding stale QFunction assembly data without modifying its state. + + @param[in] op `CeedOperator` to estimate FLOPs for + @param[out] flops Address of variable to hold FLOPs estimate + + @return An error code: 0 - success, otherwise - failure + + @ref User +**/ int CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(CeedOperator op, CeedSize *flops) { return CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, true, true, flops); } +/** + @brief Estimate the number of FLOPs required to fully assemble a linear `CeedOperator`. + + This estimate accounts for rebuilding stale QFunction assembly data without modifying its state. + + @param[in] op `CeedOperator` to estimate FLOPs for + @param[out] flops Address of variable to hold FLOPs estimate + + @return An error code: 0 - success, otherwise - failure + + @ref User +**/ int CeedOperatorLinearAssembleGetFlopsEstimate(CeedOperator op, CeedSize *flops) { return CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, false, false, flops); } diff --git a/tests/t526-operator.c b/tests/t526-operator.c index e383769d3f..0b250f64ee 100644 --- a/tests/t526-operator.c +++ b/tests/t526-operator.c @@ -125,10 +125,10 @@ int main(int argc, char **argv) { CeedQFunctionSetUserFlopsEstimate(qf_mass, 1); CeedOperatorGetFlopsEstimate(op_mass, &flop_estimate); - // Check apply estimate + // Check FLOP estimate if (flop_estimate != 3042) printf("Incorrect FLOP estimate computed, %" CeedSize_FMT " != 3042\n", flop_estimate); - // Check assembly estimates with stale QFunction data. Repeating the query must not update it. + // Check assembly FLOP estimates with stale QFunction data. Repeating the query must not update it. CeedOperatorLinearAssembleGetFlopsEstimate(op_mass, &flop_estimate); if (flop_estimate != 19416) printf("Incorrect full assembly FLOP estimate, %" CeedSize_FMT " != 19416\n", flop_estimate); CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(op_mass, &flop_estimate); @@ -136,7 +136,7 @@ int main(int argc, char **argv) { CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(op_mass, &flop_estimate); if (flop_estimate != 3234) printf("Incorrect point-block diagonal assembly FLOP estimate, %" CeedSize_FMT " != 3234\n", flop_estimate); - // Assemble once with reuse enabled, then check the cached-QFunction path + // Check QFunction data reuse reduces FLOPs estimate on second assembly CeedOperatorSetQFunctionAssemblyReuse(op_mass, true); CeedOperatorLinearAssembleGetNumEntries(op_mass, &num_entries); CeedVectorCreate(ceed, num_entries, &assembled); diff --git a/tests/t595-operator.c b/tests/t595-operator.c index cb620e1bed..5b86654871 100644 --- a/tests/t595-operator.c +++ b/tests/t595-operator.c @@ -1,6 +1,6 @@ /// @file -/// Test FLOP estimation for mass matrix operator at points -/// \test Test FLOP estimation for mass matrix operator at points +/// Test apply and assembly FLOP estimation for mass matrix operator at points +/// \test Test apply and assembly FLOP estimation for mass matrix operator at points #include "t595-operator.h" #include @@ -9,10 +9,11 @@ #include int main(int argc, char **argv) { - Ceed ceed; - CeedInt num_elem_1d = 3, num_elem = num_elem_1d * num_elem_1d, dim = 2, p = 3, q = 5; - CeedInt num_nodes = (num_elem_1d * (p - 1) + 1) * (num_elem_1d * (p - 1) + 1), num_points_per_elem = 4, num_points = num_elem * num_points_per_elem; - CeedSize flop_estimate = 0; + Ceed ceed; + CeedInt num_elem_1d = 3, num_elem = num_elem_1d * num_elem_1d, dim = 2, p = 3, q = 5; + CeedInt num_nodes = (num_elem_1d * (p - 1) + 1) * (num_elem_1d * (p - 1) + 1), num_points = 36; + CeedSize expected_apply, expected_full, expected_diagonal, flop_estimate = 0; + CeedMemType mem_type; CeedVector x_points, q_data; CeedElemRestriction elem_restriction_x_points, elem_restriction_q_data, elem_restriction_u; CeedBasis basis_x, basis_u; @@ -27,21 +28,19 @@ int main(int argc, char **argv) { { CeedScalar x_array[dim * num_points]; - for (CeedInt e = 0; e < num_elem; e++) { - for (CeedInt d = 0; d < dim; d++) { - x_array[num_points_per_elem * (e * dim + d) + 0] = 0.25; - x_array[num_points_per_elem * (e * dim + d) + 1] = d == 0 ? -0.25 : 0.25; - x_array[num_points_per_elem * (e * dim + d) + 2] = d == 0 ? 0.25 : -0.25; - x_array[num_points_per_elem * (e * dim + d) + 3] = 0.25; - } - } + for (CeedInt i = 0; i < dim * num_points; i++) x_array[i] = 0.25; CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); } { - CeedInt ind_x[num_elem + 1 + num_points]; + const CeedInt num_points_per_elem[] = {1, 2, 3, 4, 5, 6, 7, 4, 4}; + CeedInt ind_x[num_elem + 1 + num_points], offset = num_elem + 1; - for (CeedInt i = 0; i <= num_elem; i++) ind_x[i] = num_elem + 1 + i * num_points_per_elem; - for (CeedInt i = 0; i < num_points; i++) ind_x[num_elem + 1 + i] = i; + for (CeedInt e = 0; e < num_elem; e++) { + ind_x[e] = offset; + for (CeedInt i = 0; i < num_points_per_elem[e]; i++) ind_x[offset + i] = offset - num_elem - 1 + i; + offset += num_points_per_elem[e]; + } + ind_x[num_elem] = offset; CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x, &elem_restriction_x_points); CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, 1, num_points, CEED_MEM_HOST, CEED_COPY_VALUES, ind_x, &elem_restriction_q_data); @@ -102,20 +101,29 @@ int main(int argc, char **argv) { // Estimate FLOPs CeedQFunctionSetUserFlopsEstimate(qf_mass, 1); + CeedGetPreferredMemType(ceed, &mem_type); + expected_apply = mem_type == CEED_MEM_DEVICE ? 22824 : 16317; + expected_full = mem_type == CEED_MEM_DEVICE ? 11403 : 6516; + expected_diagonal = mem_type == CEED_MEM_DEVICE ? 1845 : 1089; CeedOperatorGetFlopsEstimate(op_mass, &flop_estimate); // Check output - if (flop_estimate != 16317) { + if (flop_estimate != expected_apply) { // LCOV_EXCL_START - printf("Incorrect FLOP estimate computed, %ld != 16317\n", flop_estimate); + printf("Incorrect FLOP estimate computed, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, expected_apply); // LCOV_EXCL_STOP } + // Check assembly FLOP estimates. Device backends pad each element to seven points. CeedOperatorLinearAssembleGetFlopsEstimate(op_mass, &flop_estimate); - if (flop_estimate != 6516) printf("Incorrect AtPoints full assembly FLOP estimate, %" CeedSize_FMT " != 6516\n", flop_estimate); + if (flop_estimate != expected_full) + printf("Incorrect AtPoints full assembly FLOP estimate, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, expected_full); CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(op_mass, &flop_estimate); - if (flop_estimate != 1089) printf("Incorrect AtPoints diagonal assembly FLOP estimate, %" CeedSize_FMT " != 1089\n", flop_estimate); + if (flop_estimate != expected_diagonal) + printf("Incorrect AtPoints diagonal assembly FLOP estimate, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, expected_diagonal); CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(op_mass, &flop_estimate); - if (flop_estimate != 1089) printf("Incorrect AtPoints point-block diagonal assembly FLOP estimate, %" CeedSize_FMT " != 1089\n", flop_estimate); + if (flop_estimate != expected_diagonal) + printf("Incorrect AtPoints point-block diagonal assembly FLOP estimate, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, + expected_diagonal); CeedVectorDestroy(&x_points); CeedVectorDestroy(&q_data); From a1abbe4f40e55f76c1f428a1ec66369ab9f77b2a Mon Sep 17 00:00:00 2001 From: alvaroborras <13698600+alvaroborras@users.noreply.github.com> Date: Tue, 1 Sep 2026 08:50:22 +0200 Subject: [PATCH 3/4] Tidy AtPoints FLOP estimates test Signed-off-by: alvaroborras <13698600+alvaroborras@users.noreply.github.com> --- tests/t595-operator.c | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/tests/t595-operator.c b/tests/t595-operator.c index 5b86654871..f9ea6f1b91 100644 --- a/tests/t595-operator.c +++ b/tests/t595-operator.c @@ -12,7 +12,7 @@ int main(int argc, char **argv) { Ceed ceed; CeedInt num_elem_1d = 3, num_elem = num_elem_1d * num_elem_1d, dim = 2, p = 3, q = 5; CeedInt num_nodes = (num_elem_1d * (p - 1) + 1) * (num_elem_1d * (p - 1) + 1), num_points = 36; - CeedSize expected_apply, expected_full, expected_diagonal, flop_estimate = 0; + CeedSize flop_estimate = 0; CeedMemType mem_type; CeedVector x_points, q_data; CeedElemRestriction elem_restriction_x_points, elem_restriction_q_data, elem_restriction_u; @@ -102,28 +102,29 @@ int main(int argc, char **argv) { // Estimate FLOPs CeedQFunctionSetUserFlopsEstimate(qf_mass, 1); CeedGetPreferredMemType(ceed, &mem_type); - expected_apply = mem_type == CEED_MEM_DEVICE ? 22824 : 16317; - expected_full = mem_type == CEED_MEM_DEVICE ? 11403 : 6516; - expected_diagonal = mem_type == CEED_MEM_DEVICE ? 1845 : 1089; + // Device backends pad each element to seven points. + const CeedSize expected_flops_apply = mem_type == CEED_MEM_DEVICE ? 22824 : 16317; + const CeedSize expected_flops_full = mem_type == CEED_MEM_DEVICE ? 11403 : 6516; + const CeedSize expected_flops_diagonal = mem_type == CEED_MEM_DEVICE ? 1845 : 1089; CeedOperatorGetFlopsEstimate(op_mass, &flop_estimate); // Check output - if (flop_estimate != expected_apply) { + if (flop_estimate != expected_flops_apply) { // LCOV_EXCL_START - printf("Incorrect FLOP estimate computed, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, expected_apply); + printf("Incorrect FLOP estimate computed, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, expected_flops_apply); // LCOV_EXCL_STOP } - // Check assembly FLOP estimates. Device backends pad each element to seven points. + // Check assembly FLOP estimates CeedOperatorLinearAssembleGetFlopsEstimate(op_mass, &flop_estimate); - if (flop_estimate != expected_full) - printf("Incorrect AtPoints full assembly FLOP estimate, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, expected_full); + if (flop_estimate != expected_flops_full) + printf("Incorrect AtPoints full assembly FLOP estimate, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, expected_flops_full); CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(op_mass, &flop_estimate); - if (flop_estimate != expected_diagonal) - printf("Incorrect AtPoints diagonal assembly FLOP estimate, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, expected_diagonal); + if (flop_estimate != expected_flops_diagonal) + printf("Incorrect AtPoints diagonal assembly FLOP estimate, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, expected_flops_diagonal); CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(op_mass, &flop_estimate); - if (flop_estimate != expected_diagonal) + if (flop_estimate != expected_flops_diagonal) printf("Incorrect AtPoints point-block diagonal assembly FLOP estimate, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, - expected_diagonal); + expected_flops_diagonal); CeedVectorDestroy(&x_points); CeedVectorDestroy(&q_data); From 76da6ea45b81f18e66264fbb167e74364eaca383 Mon Sep 17 00:00:00 2001 From: alvaroborras <13698600+alvaroborras@users.noreply.github.com> Date: Sun, 6 Sep 2026 09:00:45 +0200 Subject: [PATCH 4/4] interface - clarify assembly FLOP estimates Document that QFunction FLOP estimates depend on rebuild status without modifying internal state. Reuse the point-count temporary across host and device paths, follow interface error-handling conventions, and keep the CHANGELOG chronological. Signed-off-by: alvaroborras <13698600+alvaroborras@users.noreply.github.com> --- CHANGELOG.md | 4 +-- interface/ceed-preconditioning.c | 50 ++++++++++++++++++++++---------- 2 files changed, 36 insertions(+), 18 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 36b6568c33..1184a458d2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,8 +8,6 @@ On this page we provide a summary of the main API changes, new features and exam ### Interface changes -- Add `CeedOperatorLinearAssembleGetFlopsEstimate`, `CeedOperatorLinearAssembleDiagonalGetFlopsEstimate`, and - `CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate` to estimate FLOPs for linear operator assembly. - Add `bool` field type for `CeedQFunctionContext` and related interfaces to use `bool` fields. - `CEED_BASIS_COLLOCATED` removed; users should only use `CEED_BASIS_NONE`. - Remove unneeded pointer for `CeedElemRestrictionGetELayout`. @@ -22,6 +20,8 @@ On this page we provide a summary of the main API changes, new features and exam - Add `build_objects` parameter to `CeedOperatorLinearAssembleQFunctionBuildOrUpdateFallback` to allow for passing uninitialized vectors and restrictions - Move JiT helper functions only required by SYCL backends to `ceed/jit-tools-deprecated.h`. These functions will be removed when the SYCL backends are updated to reflect the improvements in the CUDA and HIP backends. +- Add `CeedOperatorLinearAssembleGetFlopsEstimate`, `CeedOperatorLinearAssembleDiagonalGetFlopsEstimate`, and + `CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate` to estimate FLOPs for linear operator assembly. ### New features diff --git a/interface/ceed-preconditioning.c b/interface/ceed-preconditioning.c index f0a5debb3c..62f8e6b1c9 100644 --- a/interface/ceed-preconditioning.c +++ b/interface/ceed-preconditioning.c @@ -2375,7 +2375,17 @@ int CeedOperatorLinearAssembleQFunctionBuildOrUpdate(CeedOperator op, CeedVector return CeedOperatorLinearAssembleQFunctionBuildOrUpdate_Core(op, *assembled == NULL, true, assembled, rstr, request); } -// Estimate the work to rebuild the assembled QFunction data, without changing its state. +/** + @brief Estimate FLOPs based on QFunction rebuild status without modifying internal state + + @param[in] op `CeedOperator` to estimate FLOPs for + @param[in] num_qpts Total number of quadrature points + @param[out] flops Address of variable to hold FLOPs estimate + + @return An error code: 0 - success, otherwise - failure + + @ref Developer +**/ static int CeedOperatorLinearAssembleQFunctionGetFlopsEstimate(CeedOperator op, CeedSize num_qpts, CeedSize *flops) { CeedInt num_elem, num_input_fields, num_active_inputs = 0; CeedQFunction qf; @@ -2439,6 +2449,18 @@ static int CeedOperatorLinearAssembleQFunctionGetFlopsEstimate(CeedOperator op, return CEED_ERROR_SUCCESS; } +/** + @brief Estimate the FLOPs required to assemble a linear `CeedOperator` + + @param[in] op `CeedOperator` to estimate FLOPs for + @param[in] is_point_block Boolean flag indicating point-block diagonal assembly + @param[in] is_diagonal Boolean flag indicating diagonal assembly + @param[out] flops Address of variable to hold FLOPs estimate + + @return An error code: 0 - success, otherwise - failure + + @ref Developer +**/ static int CeedOperatorLinearAssembleGetFlopsEstimate_Core(CeedOperator op, bool is_point_block, bool is_diagonal, CeedSize *flops) { bool is_at_points, is_composite; CeedSize num_qpts_total; @@ -2462,7 +2484,7 @@ static int CeedOperatorLinearAssembleGetFlopsEstimate_Core(CeedOperator op, bool } CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); if (is_at_points) { - CeedInt num_elem; + CeedInt num_elem, num_points; CeedMemType mem_type; CeedElemRestriction rstr_points; @@ -2470,19 +2492,12 @@ static int CeedOperatorLinearAssembleGetFlopsEstimate_Core(CeedOperator op, bool CeedCall(CeedOperatorGetNumElements(op, &num_elem)); CeedCall(CeedGetPreferredMemType(CeedOperatorReturnCeed(op), &mem_type)); if (mem_type == CEED_MEM_DEVICE) { - CeedInt max_points; - // Device backends pad every element to the maximum number of points. - CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_points)); - num_qpts_total = (CeedSize)num_elem * max_points; + CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &num_points)); + num_qpts_total = (CeedSize)num_elem * num_points; } else { - num_qpts_total = 0; - for (CeedInt i = 0; i < num_elem; i++) { - CeedInt points_in_elem; - - CeedCall(CeedElemRestrictionGetNumPointsInElement(rstr_points, i, &points_in_elem)); - num_qpts_total += points_in_elem; - } + CeedCall(CeedElemRestrictionGetNumPoints(rstr_points, &num_points)); + num_qpts_total = num_points; } CeedCall(CeedElemRestrictionDestroy(&rstr_points)); } else { @@ -2571,7 +2586,8 @@ static int CeedOperatorLinearAssembleGetFlopsEstimate_Core(CeedOperator op, bool @ref User **/ int CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(CeedOperator op, CeedSize *flops) { - return CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, false, true, flops); + CeedCall(CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, false, true, flops)); + return CEED_ERROR_SUCCESS; } /** @@ -2587,7 +2603,8 @@ int CeedOperatorLinearAssembleDiagonalGetFlopsEstimate(CeedOperator op, CeedSize @ref User **/ int CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(CeedOperator op, CeedSize *flops) { - return CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, true, true, flops); + CeedCall(CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, true, true, flops)); + return CEED_ERROR_SUCCESS; } /** @@ -2603,7 +2620,8 @@ int CeedOperatorLinearAssemblePointBlockDiagonalGetFlopsEstimate(CeedOperator op @ref User **/ int CeedOperatorLinearAssembleGetFlopsEstimate(CeedOperator op, CeedSize *flops) { - return CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, false, false, flops); + CeedCall(CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, false, false, flops)); + return CEED_ERROR_SUCCESS; } /**