Skip to content

Add FLOP estimates for operator assembly #2018

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
alvaroborras wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Add FLOP estimates for operator assembly #2018

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
9f46b6c
Add FLOP estimates for operator assembly
alvaroborras Aug 28, 2026
0239b02
Fix AtPoints assembly FLOP estimates
alvaroborras Sep 1, 2026
a1abbe4
Tidy AtPoints FLOP estimates test
alvaroborras Sep 1, 2026
76da6ea
interface - clarify assembly FLOP estimates
alvaroborras Sep 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
3 changes: 3 additions & 0 deletions include/ceed/ceed.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
249 changes: 249 additions & 0 deletions interface/ceed-preconditioning.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Comment thread
jeremylt marked this conversation as resolved.
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) {
Comment thread
jeremylt marked this conversation as resolved.
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));
Comment thread
jeremylt marked this conversation as resolved.
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`

Expand Down
34 changes: 29 additions & 5 deletions tests/t526-operator.c
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 <ceed.h>
#include <math.h>
#include <stdio.h>
Expand All @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -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);
Expand Down
22 changes: 19 additions & 3 deletions tests/t537-operator.c
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 <ceed.h>
Expand All @@ -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];

Expand Down Expand Up @@ -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);
Expand Down
Loading