diff --git a/CHANGELOG.md b/CHANGELOG.md index 62d79f46a2..1184a458d2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,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/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..62f8e6b1c9 100644 --- a/interface/ceed-preconditioning.c +++ b/interface/ceed-preconditioning.c @@ -2375,6 +2375,255 @@ int CeedOperatorLinearAssembleQFunctionBuildOrUpdate(CeedOperator op, CeedVector return CeedOperatorLinearAssembleQFunctionBuildOrUpdate_Core(op, *assembled == NULL, true, assembled, rstr, request); } +/** + @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; + 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 += num_qpts * num_active_inputs * qf_flops; + } + CeedCall(CeedQFunctionDestroy(&qf)); + 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; + + 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) { + CeedInt num_elem, num_points; + CeedMemType mem_type; + CeedElemRestriction rstr_points; + + CeedCall(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); + CeedCall(CeedOperatorGetNumElements(op, &num_elem)); + CeedCall(CeedGetPreferredMemType(CeedOperatorReturnCeed(op), &mem_type)); + if (mem_type == CEED_MEM_DEVICE) { + // Device backends pad every element to the maximum number of points. + CeedCall(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &num_points)); + num_qpts_total = (CeedSize)num_elem * num_points; + } else { + CeedCall(CeedElemRestrictionGetNumPoints(rstr_points, &num_points)); + num_qpts_total = num_points; + } + CeedCall(CeedElemRestrictionDestroy(&rstr_points)); + } else { + 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_qpts_total, 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_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(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 && 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 { + 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)); + } 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 * 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 * 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; +} + +/** + @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) { + CeedCall(CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, false, true, flops)); + return CEED_ERROR_SUCCESS; +} + +/** + @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) { + CeedCall(CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, true, true, flops)); + return CEED_ERROR_SUCCESS; +} + +/** + @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) { + CeedCall(CeedOperatorLinearAssembleGetFlopsEstimate_Core(op, false, false, flops)); + return CEED_ERROR_SUCCESS; +} + /** @brief Assemble the diagonal of a square linear `CeedOperator` diff --git a/tests/t526-operator.c b/tests/t526-operator.c index 8e68ab89b3..0b250f64ee 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 FLOP estimate if (flop_estimate != 3042) printf("Incorrect FLOP estimate computed, %" CeedSize_FMT " != 3042\n", flop_estimate); + // 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); + 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); + + // Check QFunction data reuse reduces FLOPs estimate on second assembly + 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..f9ea6f1b91 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; + 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 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,14 +101,30 @@ int main(int argc, char **argv) { // Estimate FLOPs CeedQFunctionSetUserFlopsEstimate(qf_mass, 1); + CeedGetPreferredMemType(ceed, &mem_type); + // 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 != 16317) { + if (flop_estimate != expected_flops_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_flops_apply); // LCOV_EXCL_STOP } + // Check assembly FLOP estimates + CeedOperatorLinearAssembleGetFlopsEstimate(op_mass, &flop_estimate); + 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_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_flops_diagonal) + printf("Incorrect AtPoints point-block diagonal assembly FLOP estimate, %" CeedSize_FMT " != %" CeedSize_FMT "\n", flop_estimate, + expected_flops_diagonal); CeedVectorDestroy(&x_points); CeedVectorDestroy(&q_data);