Skip to content

Commit

Permalink
CeedVector/Preconditioning: fix CeedInt loop vars to CeedSize
Browse files Browse the repository at this point in the history 
While 32-bit is sufficient for CeedElemRestriction, a Vector is used
to store matrix entries and the number of entries can overflow 32-bit
even for a small number of dofs. For example, 85k Q3 fluid elements is
enough to overflow.

Reported-by: Ken Jansen
  • Loading branch information
@jedbrown
jedbrown committed Jun 23, 2023
1 parent f95bab6 commit e911285
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 25 deletions.
32 changes: 17 additions & 15 deletions interface/ceed-preconditioning.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -340,7 +340,7 @@ static inline int CeedSingleOperatorAssembleAddDiagonal_Core(CeedOperator op, Ce

// Compute the diagonal of B^T D B
// Each element
for (CeedInt e = 0; e < num_elem; e++) {
for (CeedSize e = 0; e < num_elem; e++) {
// Each basis eval mode pair
CeedInt d_out = 0, q_comp_out;
CeedEvalMode eval_mode_out_prev = CEED_EVAL_NONE;
Expand Down Expand Up @@ -373,7 +373,7 @@ static inline int CeedSingleOperatorAssembleAddDiagonal_Core(CeedOperator op, Ce
if (is_pointblock) {
// Point Block Diagonal
for (CeedInt c_in = 0; c_in < num_components; c_in++) {
const CeedInt c_offset = (eval_mode_offsets_in[b][e_in] + c_in) * num_output_components + eval_mode_offsets_out[b][e_out] + c_out;
const CeedSize c_offset = (eval_mode_offsets_in[b][e_in] + c_in) * num_output_components + eval_mode_offsets_out[b][e_out] + c_out;
const CeedScalar qf_value = assembled_qf_array[q * layout[0] + c_offset * layout[1] + e * layout[2]];
for (CeedInt n = 0; n < num_nodes; n++) {
elem_diag_array[((e * num_components + c_out) * num_components + c_in) * num_nodes + n] +=
Expand Down Expand Up @@ -491,7 +491,7 @@ static int CeedSingleOperatorAssembleSymbolic(CeedOperator op, CeedInt offset, C
CeedCall(CeedVectorDestroy(&index_vec));

// Determine i, j locations for element matrices
CeedInt count = 0;
CeedSize count = 0;
for (CeedInt e = 0; e < num_elem; e++) {
for (CeedInt comp_in = 0; comp_in < num_comp; comp_in++) {
for (CeedInt comp_out = 0; comp_out < num_comp; comp_out++) {
Expand Down Expand Up @@ -614,22 +614,22 @@ static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset, CeedVecto
const CeedScalar *B_mat_in = B_mats_in[0], *B_mat_out = B_mats_out[0];
CeedScalar BTD_mat[elem_size * num_qpts * num_eval_modes_in[0]];
CeedScalar elem_mat[elem_size * elem_size];
CeedInt count = 0;
CeedSize count = 0;
CeedScalar *vals;
CeedCall(CeedVectorGetArray(values, CEED_MEM_HOST, &vals));
for (CeedInt e = 0; e < num_elem; e++) {
for (CeedSize e = 0; e < num_elem; e++) {
for (CeedInt comp_in = 0; comp_in < num_comp; comp_in++) {
for (CeedInt comp_out = 0; comp_out < num_comp; comp_out++) {
// Compute B^T*D
for (CeedInt n = 0; n < elem_size; n++) {
for (CeedInt q = 0; q < num_qpts; q++) {
for (CeedSize n = 0; n < elem_size; n++) {
for (CeedSize q = 0; q < num_qpts; q++) {
for (CeedInt e_in = 0; e_in < num_eval_modes_in[0]; e_in++) {
const CeedInt btd_index = n * (num_qpts * num_eval_modes_in[0]) + (num_eval_modes_in[0] * q + e_in);
const CeedSize btd_index = n * (num_qpts * num_eval_modes_in[0]) + (num_eval_modes_in[0] * q + e_in);
CeedScalar sum = 0.0;
for (CeedInt e_out = 0; e_out < num_eval_modes_out[0]; e_out++) {
const CeedInt b_out_index = (num_eval_modes_out[0] * q + e_out) * elem_size + n;
const CeedInt eval_mode_index = ((e_in * num_comp + comp_in) * num_eval_modes_out[0] + e_out) * num_comp + comp_out;
const CeedInt qf_index = q * layout_qf[0] + eval_mode_index * layout_qf[1] + e * layout_qf[2];
const CeedSize b_out_index = (num_eval_modes_out[0] * q + e_out) * elem_size + n;
const CeedSize eval_mode_index = ((e_in * num_comp + comp_in) * num_eval_modes_out[0] + e_out) * num_comp + comp_out;
const CeedSize qf_index = q * layout_qf[0] + eval_mode_index * layout_qf[1] + e * layout_qf[2];
sum += B_mat_out[b_out_index] * assembled_qf_array[qf_index];
}
BTD_mat[btd_index] = sum;
Expand Down Expand Up @@ -668,7 +668,7 @@ static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset, CeedVecto
@ref Utility
**/
static int CeedSingleOperatorAssemblyCountEntries(CeedOperator op, CeedInt *num_entries) {
static int CeedSingleOperatorAssemblyCountEntries(CeedOperator op, CeedSize *num_entries) {
bool is_composite;
CeedElemRestriction rstr;
CeedInt num_elem, elem_size, num_comp;
Expand All @@ -679,7 +679,7 @@ static int CeedSingleOperatorAssemblyCountEntries(CeedOperator op, CeedInt *num_
CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem));
CeedCall(CeedElemRestrictionGetElementSize(rstr, &elem_size));
CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp));
*num_entries = elem_size * num_comp * elem_size * num_comp * num_elem;
*num_entries = (CeedSize)elem_size * num_comp * elem_size * num_comp * num_elem;

return CEED_ERROR_SUCCESS;
}
Expand Down Expand Up @@ -1839,7 +1839,8 @@ matrix in entry (i, j).
@ref User
**/
int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols) {
CeedInt num_suboperators, single_entries;
CeedInt num_suboperators;
CeedSize single_entries;
CeedOperator *sub_operators;
bool is_composite;
CeedCall(CeedOperatorCheckReady(op));
Expand Down Expand Up @@ -1915,7 +1916,8 @@ matrix in entry (i, j).
@ref User
**/
int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) {
CeedInt num_suboperators, single_entries = 0;
CeedInt num_suboperators;
CeedSize single_entries = 0;
CeedOperator *sub_operators;
bool is_composite;
CeedCall(CeedOperatorCheckReady(op));
Expand Down
20 changes: 10 additions & 10 deletions interface/ceed-vector.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -274,7 +274,7 @@ int CeedVectorSetValue(CeedVector vec, CeedScalar value) {
} else {
CeedScalar *array;
CeedCall(CeedVectorGetArrayWrite(vec, CEED_MEM_HOST, &array));
for (CeedInt i = 0; i < vec->length; i++) array[i] = value;
for (CeedSize i = 0; i < vec->length; i++) array[i] = value;
CeedCall(CeedVectorRestoreArray(vec, &array));
}
vec->state += 2;
Expand Down Expand Up @@ -503,17 +503,17 @@ int CeedVectorNorm(CeedVector vec, CeedNormType norm_type, CeedScalar *norm) {
*norm = 0.;
switch (norm_type) {
case CEED_NORM_1:
for (CeedInt i = 0; i < vec->length; i++) {
for (CeedSize i = 0; i < vec->length; i++) {
*norm += fabs(array[i]);
}
break;
case CEED_NORM_2:
for (CeedInt i = 0; i < vec->length; i++) {
for (CeedSize i = 0; i < vec->length; i++) {
*norm += fabs(array[i]) * fabs(array[i]);
}
break;
case CEED_NORM_MAX:
for (CeedInt i = 0; i < vec->length; i++) {
for (CeedSize i = 0; i < vec->length; i++) {
const CeedScalar abs_v_i = fabs(array[i]);
*norm = *norm > abs_v_i ? *norm : abs_v_i;
}
Expand Down Expand Up @@ -550,7 +550,7 @@ int CeedVectorScale(CeedVector x, CeedScalar alpha) {

// Default implementation
CeedCall(CeedVectorGetArrayWrite(x, CEED_MEM_HOST, &x_array));
for (CeedInt i = 0; i < n_x; i++) x_array[i] *= alpha;
for (CeedSize i = 0; i < n_x; i++) x_array[i] *= alpha;
CeedCall(CeedVectorRestoreArray(x, &x_array));

return CEED_ERROR_SUCCESS;
Expand Down Expand Up @@ -603,7 +603,7 @@ int CeedVectorAXPY(CeedVector y, CeedScalar alpha, CeedVector x) {
assert(x_array);
assert(y_array);

for (CeedInt i = 0; i < n_y; i++) y_array[i] += alpha * x_array[i];
for (CeedSize i = 0; i < n_y; i++) y_array[i] += alpha * x_array[i];

CeedCall(CeedVectorRestoreArray(y, &y_array));
CeedCall(CeedVectorRestoreArrayRead(x, &x_array));
Expand Down Expand Up @@ -659,7 +659,7 @@ int CeedVectorAXPBY(CeedVector y, CeedScalar alpha, CeedScalar beta, CeedVector
assert(x_array);
assert(y_array);

for (CeedInt i = 0; i < n_y; i++) y_array[i] += alpha * x_array[i] + beta * y_array[i];
for (CeedSize i = 0; i < n_y; i++) y_array[i] += alpha * x_array[i] + beta * y_array[i];

CeedCall(CeedVectorRestoreArray(y, &y_array));
CeedCall(CeedVectorRestoreArrayRead(x, &x_array));
Expand Down Expand Up @@ -732,7 +732,7 @@ int CeedVectorPointwiseMult(CeedVector w, CeedVector x, CeedVector y) {
assert(x_array);
assert(y_array);

for (CeedInt i = 0; i < n_w; i++) w_array[i] = x_array[i] * y_array[i];
for (CeedSize i = 0; i < n_w; i++) w_array[i] = x_array[i] * y_array[i];

if (y != w && y != x) CeedCall(CeedVectorRestoreArrayRead(y, &y_array));
if (x != w) CeedCall(CeedVectorRestoreArrayRead(x, &x_array));
Expand Down Expand Up @@ -769,7 +769,7 @@ int CeedVectorReciprocal(CeedVector vec) {
CeedCall(CeedVectorGetLength(vec, &len));
CeedScalar *array;
CeedCall(CeedVectorGetArrayWrite(vec, CEED_MEM_HOST, &array));
for (CeedInt i = 0; i < len; i++) {
for (CeedSize i = 0; i < len; i++) {
if (fabs(array[i]) > CEED_EPSILON) array[i] = 1. / array[i];
}

Expand Down Expand Up @@ -809,7 +809,7 @@ int CeedVectorViewRange(CeedVector vec, CeedSize start, CeedSize stop, CeedInt s

snprintf(fmt, sizeof fmt, " %s\n", fp_fmt ? fp_fmt : "%g");
CeedCall(CeedVectorGetArrayRead(vec, CEED_MEM_HOST, &x));
for (CeedInt i = start; step > 0 ? (i < stop) : (i > stop); i += step) fprintf(stream, fmt, x[i]);
for (CeedSize i = start; step > 0 ? (i < stop) : (i > stop); i += step) fprintf(stream, fmt, x[i]);
CeedCall(CeedVectorRestoreArrayRead(vec, &x));
if (stop != vec->length) fprintf(stream, " ...\n");

Expand Down

0 comments on commit e911285

Please sign in to comment.