Skip to content

Commit

Permalink
Merge pull request #970 from CEED/jeremy/interface-audit
Browse files Browse the repository at this point in the history 
Audit object encapsulation
  • Loading branch information
@jeremylt
jeremylt authored May 26, 2022
2 parents e3d63e5 + 21d1e94 commit fe329bd
Show file tree
Hide file tree
Showing 8 changed files with 119 additions and 92 deletions.
1 change: 1 addition & 0 deletions include/ceed/backend.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,7 @@ CEED_EXTERN int CeedOperatorGetActiveBasis(CeedOperator op,
CEED_EXTERN int CeedOperatorGetActiveElemRestriction(CeedOperator op, CeedElemRestriction *active_rstr);
CEED_EXTERN int CeedGetOperatorFallbackResource(Ceed ceed,
const char **resource);
CEED_EXTERN int CeedGetOperatorFallbackCeed(Ceed ceed, Ceed *fallback_ceed);
CEED_EXTERN int CeedSetOperatorFallbackResource(Ceed ceed,
const char *resource);
CEED_EXTERN int CeedGetOperatorFallbackParentCeed(Ceed ceed, Ceed *parent);
Expand Down
20 changes: 10 additions & 10 deletions interface/ceed-basis.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,12 +158,12 @@ static int CeedGivensRotation(CeedScalar *A, CeedScalar c, CeedScalar s,
**/
static int CeedScalarView(const char *name, const char *fp_fmt, CeedInt m,
CeedInt n, const CeedScalar *a, FILE *stream) {
for (int i=0; i<m; i++) {
for (CeedInt i=0; i<m; i++) {
if (m > 1)
fprintf(stream, "%12s[%d]:", name, i);
else
fprintf(stream, "%12s:", name);
for (int j=0; j<n; j++)
for (CeedInt j=0; j<n; j++)
fprintf(stream, fp_fmt, fabs(a[i*n+j]) > 1E-14 ? a[i*n+j] : 0);
fputs("\n", stream);
}
Expand Down Expand Up @@ -1366,14 +1366,14 @@ int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d,
// Allocate
CeedScalar P0, P1, P2, dP2, xi, wi, PI = 4.0*atan(1.0);
// Build q_ref_1d, q_weight_1d
for (int i = 0; i <= Q/2; i++) {
for (CeedInt i = 0; i <= Q/2; i++) {
// Guess
xi = cos(PI*(CeedScalar)(2*i+1)/((CeedScalar)(2*Q)));
// Pn(xi)
P0 = 1.0;
P1 = xi;
P2 = 0.0;
for (int j = 2; j <= Q; j++) {
for (CeedInt j = 2; j <= Q; j++) {
P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j));
P0 = P1;
P1 = P2;
Expand All @@ -1382,10 +1382,10 @@ int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d,
dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0);
xi = xi-P2/dP2;
// Newton to convergence
for (int k=0; k<100 && fabs(P2)>10*CEED_EPSILON; k++) {
for (CeedInt k=0; k<100 && fabs(P2)>10*CEED_EPSILON; k++) {
P0 = 1.0;
P1 = xi;
for (int j = 2; j <= Q; j++) {
for (CeedInt j = 2; j <= Q; j++) {
P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j));
P0 = P1;
P1 = P2;
Expand Down Expand Up @@ -1434,14 +1434,14 @@ int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d,
q_ref_1d[0] = -1.0;
q_ref_1d[Q-1] = 1.0;
// Interior
for (int i = 1; i <= (Q-1)/2; i++) {
for (CeedInt i = 1; i <= (Q-1)/2; i++) {
// Guess
xi = cos(PI*(CeedScalar)(i)/(CeedScalar)(Q-1));
// Pn(xi)
P0 = 1.0;
P1 = xi;
P2 = 0.0;
for (int j = 2; j < Q; j++) {
for (CeedInt j = 2; j < Q; j++) {
P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j));
P0 = P1;
P1 = P2;
Expand All @@ -1451,10 +1451,10 @@ int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d,
d2P2 = (2*xi*dP2 - (CeedScalar)(Q*(Q-1))*P2)/(1.0-xi*xi);
xi = xi-dP2/d2P2;
// Newton to convergence
for (int k=0; k<100 && fabs(dP2)>10*CEED_EPSILON; k++) {
for (CeedInt k=0; k<100 && fabs(dP2)>10*CEED_EPSILON; k++) {
P0 = 1.0;
P1 = xi;
for (int j = 2; j < Q; j++) {
for (CeedInt j = 2; j < Q; j++) {
P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j));
P0 = P1;
P1 = P2;
Expand Down
14 changes: 7 additions & 7 deletions interface/ceed-elemrestriction.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,8 +43,8 @@ int CeedPermutePadOffsets(const CeedInt *offsets, CeedInt *blk_offsets,
CeedInt num_blk, CeedInt num_elem, CeedInt blk_size,
CeedInt elem_size) {
for (CeedInt e=0; e<num_blk*blk_size; e+=blk_size)
for (int j=0; j<blk_size; j++)
for (int k=0; k<elem_size; k++)
for (CeedInt j=0; j<blk_size; j++)
for (CeedInt k=0; k<elem_size; k++)
blk_offsets[e*elem_size + k*blk_size + j]
= offsets[CeedIntMin(e+j,num_elem-1)*elem_size + k];
return CEED_ERROR_SUCCESS;
Expand Down Expand Up @@ -77,7 +77,7 @@ int CeedElemRestrictionGetStrides(CeedElemRestriction rstr,
"ElemRestriction has no stride data");
// LCOV_EXCL_STOP

for (int i=0; i<3; i++)
for (CeedInt i=0; i<3; i++)
(*strides)[i] = rstr->strides[i];
return CEED_ERROR_SUCCESS;
}
Expand Down Expand Up @@ -205,7 +205,7 @@ int CeedElemRestrictionGetELayout(CeedElemRestriction rstr,
"ElemRestriction has no layout data");
// LCOV_EXCL_STOP

for (int i=0; i<3; i++)
for (CeedInt i=0; i<3; i++)
(*layout)[i] = rstr->layout[i];
return CEED_ERROR_SUCCESS;
}
Expand All @@ -227,7 +227,7 @@ int CeedElemRestrictionGetELayout(CeedElemRestriction rstr,
**/
int CeedElemRestrictionSetELayout(CeedElemRestriction rstr,
CeedInt layout[3]) {
for (int i = 0; i<3; i++)
for (CeedInt i = 0; i<3; i++)
rstr->layout[i] = layout[i];
return CEED_ERROR_SUCCESS;
}
Expand Down Expand Up @@ -570,7 +570,7 @@ int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt num_elem,
(*rstr)->blk_size = 1;
(*rstr)->is_oriented = 0;
ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr);
for (int i=0; i<3; i++)
for (CeedInt i=0; i<3; i++)
(*rstr)->strides[i] = strides[i];
ierr = ceed->ElemRestrictionCreate(CEED_MEM_HOST, CEED_OWN_POINTER, NULL,
*rstr);
Expand Down Expand Up @@ -765,7 +765,7 @@ int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt num_elem,
(*rstr)->blk_size = blk_size;
(*rstr)->is_oriented = 0;
ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr);
for (int i=0; i<3; i++)
for (CeedInt i=0; i<3; i++)
(*rstr)->strides[i] = strides[i];
ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER,
NULL, *rstr); CeedChk(ierr);
Expand Down
4 changes: 2 additions & 2 deletions interface/ceed-operator.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,7 +213,7 @@ int CeedOperatorSingleView(CeedOperator op, bool sub, FILE *stream) {
int CeedOperatorGetActiveBasis(CeedOperator op, CeedBasis *active_basis) {
*active_basis = NULL;
if (op->is_composite) return CEED_ERROR_SUCCESS;
for (int i = 0; i < op->qf->num_input_fields; i++)
for (CeedInt i = 0; i < op->qf->num_input_fields; i++)
if (op->input_fields[i]->vec == CEED_VECTOR_ACTIVE) {
*active_basis = op->input_fields[i]->basis;
break;
Expand Down Expand Up @@ -245,7 +245,7 @@ int CeedOperatorGetActiveElemRestriction(CeedOperator op,
CeedElemRestriction *active_rstr) {
*active_rstr = NULL;
if (op->is_composite) return CEED_ERROR_SUCCESS;
for (int i = 0; i < op->qf->num_input_fields; i++)
for (CeedInt i = 0; i < op->qf->num_input_fields; i++)
if (op->input_fields[i]->vec == CEED_VECTOR_ACTIVE) {
*active_rstr = op->input_fields[i]->elem_restr;
break;
Expand Down
115 changes: 51 additions & 64 deletions interface/ceed-preconditioning.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,42 +34,25 @@
**/
int CeedOperatorCreateFallback(CeedOperator op) {
int ierr;
Ceed fallback_ceed;

// Check not already created
if (op->op_fallback) return CEED_ERROR_SUCCESS;

// Fallback Ceed
const char *resource, *fallback_resource;
ierr = CeedGetResource(op->ceed, &resource); CeedChk(ierr);
ierr = CeedGetOperatorFallbackResource(op->ceed, &fallback_resource);
CeedChk(ierr);
if (!strcmp(resource, fallback_resource))
// LCOV_EXCL_START
return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED,
"Backend %s cannot create an operator"
"fallback to resource %s", resource, fallback_resource);
// LCOV_EXCL_STOP

// Fallback Ceed
if (!op->ceed->op_fallback_ceed) {
Ceed ceed_ref;
ierr = CeedInit(fallback_resource, &ceed_ref); CeedChk(ierr);
ceed_ref->op_fallback_parent = op->ceed;
ceed_ref->Error = op->ceed->Error;
op->ceed->op_fallback_ceed = ceed_ref;
}
ierr = CeedGetOperatorFallbackCeed(op->ceed, &fallback_ceed); CeedChk(ierr);

// Clone Op
CeedOperator op_fallback;
if (op->is_composite) {
ierr = CeedCompositeOperatorCreate(op->ceed->op_fallback_ceed, &op_fallback);
ierr = CeedCompositeOperatorCreate(fallback_ceed, &op_fallback);
CeedChk(ierr);
for (CeedInt i = 0; i < op->num_suboperators; i++) {
ierr = CeedCompositeOperatorAddSub(op_fallback, op->sub_operators[i]);
CeedChk(ierr);
}
} else {
ierr = CeedOperatorCreate(op->ceed->op_fallback_ceed, op->qf, op->dqf, op->dqfT,
ierr = CeedOperatorCreate(fallback_ceed, op->qf, op->dqf, op->dqfT,
&op_fallback); CeedChk(ierr);
for (CeedInt i = 0; i < op->qf->num_input_fields; i++) {
ierr = CeedOperatorSetField(op_fallback, op->input_fields[i]->field_name,
Expand Down Expand Up @@ -500,11 +483,11 @@ static int CeedSingleOperatorAssembleSymbolic(CeedOperator op, CeedInt offset,

// Determine i, j locations for element matrices
CeedInt count = 0;
for (int e = 0; e < num_elem; ++e) {
for (int comp_in = 0; comp_in < num_comp; ++comp_in) {
for (int comp_out = 0; comp_out < num_comp; ++comp_out) {
for (int i = 0; i < elem_size; ++i) {
for (int j = 0; j < elem_size; ++j) {
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) {
for (CeedInt i = 0; i < elem_size; ++i) {
for (CeedInt j = 0; j < elem_size; ++j) {
const CeedInt elem_dof_index_row = (i)*layout_er[0] +
(comp_out)*layout_er[1] + e*layout_er[2];
const CeedInt elem_dof_index_col = (j)*layout_er[0] +
Expand Down Expand Up @@ -685,31 +668,32 @@ static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset,
num_qpts]; // logically 3-tensor
CeedScalar BTD[elem_size * num_qpts*num_eval_mode_in];
CeedScalar elem_mat[elem_size * elem_size];
int count = 0;
CeedInt count = 0;
CeedScalar *vals;
ierr = CeedVectorGetArrayWrite(values, CEED_MEM_HOST, &vals); CeedChk(ierr);
for (int e = 0; e < num_elem; ++e) {
for (int comp_in = 0; comp_in < num_comp; ++comp_in) {
for (int comp_out = 0; comp_out < num_comp; ++comp_out) {
for (int ell = 0; ell < (num_qpts * num_eval_mode_in) * elem_size; ++ell) {
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) {
for (CeedInt ell = 0; ell < (num_qpts * num_eval_mode_in) * elem_size; ++ell) {
B_mat_in[ell] = 0.0;
}
for (int ell = 0; ell < (num_qpts * num_eval_mode_out) * elem_size; ++ell) {
for (CeedInt ell = 0; ell < (num_qpts * num_eval_mode_out) * elem_size; ++ell) {
B_mat_out[ell] = 0.0;
}
// Store block-diagonal D matrix as collection of small dense blocks
for (int ell = 0; ell < num_eval_mode_in*num_eval_mode_out*num_qpts; ++ell) {
for (CeedInt ell = 0; ell < num_eval_mode_in*num_eval_mode_out*num_qpts;
++ell) {
D_mat[ell] = 0.0;
}
// form element matrix itself (for each block component)
for (int ell = 0; ell < elem_size*elem_size; ++ell) {
for (CeedInt ell = 0; ell < elem_size*elem_size; ++ell) {
elem_mat[ell] = 0.0;
}
for (int q = 0; q < num_qpts; ++q) {
for (int n = 0; n < elem_size; ++n) {
for (CeedInt q = 0; q < num_qpts; ++q) {
for (CeedInt n = 0; n < elem_size; ++n) {
CeedInt d_in = -1;
for (int e_in = 0; e_in < num_eval_mode_in; ++e_in) {
const int qq = num_eval_mode_in*q;
for (CeedInt e_in = 0; e_in < num_eval_mode_in; ++e_in) {
const CeedInt qq = num_eval_mode_in*q;
if (eval_mode_in[e_in] == CEED_EVAL_INTERP) {
B_mat_in[(qq+e_in)*elem_size + n] += interp_in[q * elem_size + n];
} else if (eval_mode_in[e_in] == CEED_EVAL_GRAD) {
Expand All @@ -723,8 +707,8 @@ static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset,
}
}
CeedInt d_out = -1;
for (int e_out = 0; e_out < num_eval_mode_out; ++e_out) {
const int qq = num_eval_mode_out*q;
for (CeedInt e_out = 0; e_out < num_eval_mode_out; ++e_out) {
const CeedInt qq = num_eval_mode_out*q;
if (eval_mode_out[e_out] == CEED_EVAL_INTERP) {
B_mat_out[(qq+e_out)*elem_size + n] += interp_in[q * elem_size + n];
} else if (eval_mode_out[e_out] == CEED_EVAL_GRAD) {
Expand All @@ -738,25 +722,26 @@ static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset,
}
}
}
for (int ei = 0; ei < num_eval_mode_out; ++ei) {
for (int ej = 0; ej < num_eval_mode_in; ++ej) {
const int eval_mode_index = ((ei*num_comp+comp_in)*num_eval_mode_in+ej)*num_comp
+comp_out;
const int index = q*layout_qf[0] + eval_mode_index*layout_qf[1] +
e*layout_qf[2];
for (CeedInt ei = 0; ei < num_eval_mode_out; ++ei) {
for (CeedInt ej = 0; ej < num_eval_mode_in; ++ej) {
const CeedInt eval_mode_index = ((ei*num_comp+comp_in)*num_eval_mode_in+ej)
*num_comp
+comp_out;
const CeedInt index = q*layout_qf[0] + eval_mode_index*layout_qf[1] +
e*layout_qf[2];
D_mat[(ei*num_eval_mode_in+ej)*num_qpts + q] += assembled_qf_array[index];
}
}
}
// Compute B^T*D
for (int ell = 0; ell < elem_size*num_qpts*num_eval_mode_in; ++ell) {
for (CeedInt ell = 0; ell < elem_size*num_qpts*num_eval_mode_in; ++ell) {
BTD[ell] = 0.0;
}
for (int j = 0; j<elem_size; ++j) {
for (int q = 0; q<num_qpts; ++q) {
int qq = num_eval_mode_out*q;
for (int ei = 0; ei < num_eval_mode_in; ++ei) {
for (int ej = 0; ej < num_eval_mode_out; ++ej) {
for (CeedInt j = 0; j<elem_size; ++j) {
for (CeedInt q = 0; q<num_qpts; ++q) {
const CeedInt qq = num_eval_mode_out*q;
for (CeedInt ei = 0; ei < num_eval_mode_in; ++ei) {
for (CeedInt ej = 0; ej < num_eval_mode_out; ++ej) {
BTD[j*(num_qpts*num_eval_mode_in) + (qq+ei)] +=
B_mat_out[(qq+ej)*elem_size + j] * D_mat[(ei*num_eval_mode_in+ej)*num_qpts + q];
}
Expand All @@ -768,8 +753,8 @@ static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset,
elem_size, num_qpts*num_eval_mode_in); CeedChk(ierr);

// put element matrix in coordinate data structure
for (int i = 0; i < elem_size; ++i) {
for (int j = 0; j < elem_size; ++j) {
for (CeedInt i = 0; i < elem_size; ++i) {
for (CeedInt j = 0; j < elem_size; ++j) {
vals[offset + count] = elem_mat[i*elem_size + j];
count++;
}
Expand Down Expand Up @@ -861,7 +846,7 @@ static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine,
op_coarse); CeedChk(ierr);
CeedElemRestriction rstr_fine = NULL;
// -- Clone input fields
for (int i = 0; i < op_fine->qf->num_input_fields; i++) {
for (CeedInt i = 0; i < op_fine->qf->num_input_fields; i++) {
if (op_fine->input_fields[i]->vec == CEED_VECTOR_ACTIVE) {
rstr_fine = op_fine->input_fields[i]->elem_restr;
ierr = CeedOperatorSetField(*op_coarse, op_fine->input_fields[i]->field_name,
Expand All @@ -875,7 +860,7 @@ static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine,
}
}
// -- Clone output fields
for (int i = 0; i < op_fine->qf->num_output_fields; i++) {
for (CeedInt i = 0; i < op_fine->qf->num_output_fields; i++) {
if (op_fine->output_fields[i]->vec == CEED_VECTOR_ACTIVE) {
ierr = CeedOperatorSetField(*op_coarse, op_fine->output_fields[i]->field_name,
rstr_coarse, basis_coarse, CEED_VECTOR_ACTIVE);
Expand Down Expand Up @@ -1722,7 +1707,7 @@ int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries,
if (is_composite) {
ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr);
ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr);
for (int k = 0; k < num_suboperators; ++k) {
for (CeedInt k = 0; k < num_suboperators; ++k) {
ierr = CeedSingleOperatorAssemblyCountEntries(sub_operators[k],
&single_entries); CeedChk(ierr);
*num_entries += single_entries;
Expand All @@ -1740,7 +1725,7 @@ int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries,
if (is_composite) {
ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr);
ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr);
for (int k = 0; k < num_suboperators; ++k) {
for (CeedInt k = 0; k < num_suboperators; ++k) {
ierr = CeedSingleOperatorAssembleSymbolic(sub_operators[k], offset, *rows,
*cols); CeedChk(ierr);
ierr = CeedSingleOperatorAssemblyCountEntries(sub_operators[k],
Expand Down Expand Up @@ -1812,7 +1797,7 @@ int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) {
if (is_composite) {
ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChk(ierr);
ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChk(ierr);
for (int k = 0; k < num_suboperators; ++k) {
for (CeedInt k = 0; k < num_suboperators; ++k) {
ierr = CeedSingleOperatorAssemble(sub_operators[k], offset, values);
CeedChk(ierr);
ierr = CeedSingleOperatorAssemblyCountEntries(sub_operators[k],
Expand Down Expand Up @@ -1893,14 +1878,16 @@ int CeedOperatorMultigridLevelCreate(CeedOperator op_fine,
ierr = CeedMalloc(Q*P_c, &interp_c); CeedChk(ierr);
ierr = CeedCalloc(P_c*P_f, &interp_c_to_f); CeedChk(ierr);
ierr = CeedMalloc(Q, &tau); CeedChk(ierr);
const CeedScalar *interp_f_source = NULL, *interp_c_source = NULL;
if (is_tensor_f) {
memcpy(interp_f, basis_fine->interp_1d, Q*P_f*sizeof basis_fine->interp_1d[0]);
memcpy(interp_c, basis_coarse->interp_1d,
Q*P_c*sizeof basis_coarse->interp_1d[0]);
ierr = CeedBasisGetInterp1D(basis_fine, &interp_f_source); CeedChk(ierr);
ierr = CeedBasisGetInterp1D(basis_coarse, &interp_c_source); CeedChk(ierr);
} else {
memcpy(interp_f, basis_fine->interp, Q*P_f*sizeof basis_fine->interp[0]);
memcpy(interp_c, basis_coarse->interp, Q*P_c*sizeof basis_coarse->interp[0]);
ierr = CeedBasisGetInterp(basis_fine, &interp_f_source); CeedChk(ierr);
ierr = CeedBasisGetInterp(basis_coarse, &interp_c_source); CeedChk(ierr);
}
memcpy(interp_f, interp_f_source, Q*P_f*sizeof interp_f_source[0]);
memcpy(interp_c, interp_c_source, Q*P_c*sizeof interp_c_source[0]);

// -- QR Factorization, interp_f = Q R
ierr = CeedQRFactorization(ceed, interp_f, tau, Q, P_f); CeedChk(ierr);
Expand Down
Loading

0 comments on commit fe329bd

Please sign in to comment.