1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 #include <ceed.h>
9 #include <ceed/backend.h>
10 #include <stdbool.h>
11 #include <stddef.h>
12 #include <stdint.h>
13
14 #include "ceed-blocked.h"
15
16 //------------------------------------------------------------------------------
17 // Setup Input/Output Fields
18 //------------------------------------------------------------------------------
CeedOperatorSetupFields_Blocked(CeedQFunction qf,CeedOperator op,bool is_input,bool * skip_rstr,CeedInt * e_data_out_indices,bool * apply_add_basis,const CeedInt block_size,CeedElemRestriction * block_rstr,CeedVector * e_vecs_full,CeedVector * e_vecs,CeedVector * q_vecs,CeedInt start_e,CeedInt num_fields,CeedInt Q)19 static int CeedOperatorSetupFields_Blocked(CeedQFunction qf, CeedOperator op, bool is_input, bool *skip_rstr, CeedInt *e_data_out_indices,
20 bool *apply_add_basis, const CeedInt block_size, CeedElemRestriction *block_rstr, CeedVector *e_vecs_full,
21 CeedVector *e_vecs, CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, CeedInt Q) {
22 Ceed ceed;
23 CeedSize e_size, q_size;
24 CeedInt num_comp, size, P;
25 CeedQFunctionField *qf_fields;
26 CeedOperatorField *op_fields;
27
28 {
29 Ceed ceed_parent;
30
31 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
32 CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
33 CeedCallBackend(CeedReferenceCopy(ceed_parent, &ceed));
34 CeedCallBackend(CeedDestroy(&ceed_parent));
35 }
36 if (is_input) {
37 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL));
38 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL));
39 } else {
40 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields));
41 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields));
42 }
43
44 // Loop over fields
45 for (CeedInt i = 0; i < num_fields; i++) {
46 CeedEvalMode eval_mode;
47 CeedBasis basis;
48
49 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
50 if (eval_mode != CEED_EVAL_WEIGHT) {
51 Ceed ceed_rstr;
52 CeedSize l_size;
53 CeedInt num_elem, elem_size, comp_stride;
54 CeedRestrictionType rstr_type;
55 CeedElemRestriction rstr;
56
57 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr));
58 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed_rstr));
59 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr, &num_elem));
60 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr, &elem_size));
61 CeedCallBackend(CeedElemRestrictionGetLVectorSize(rstr, &l_size));
62 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp));
63 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride));
64
65 CeedCallBackend(CeedElemRestrictionGetType(rstr, &rstr_type));
66 switch (rstr_type) {
67 case CEED_RESTRICTION_STANDARD: {
68 const CeedInt *offsets = NULL;
69
70 CeedCallBackend(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets));
71 CeedCallBackend(CeedElemRestrictionCreateBlocked(ceed_rstr, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, CEED_MEM_HOST,
72 CEED_COPY_VALUES, offsets, &block_rstr[i + start_e]));
73 CeedCallBackend(CeedElemRestrictionRestoreOffsets(rstr, &offsets));
74 } break;
75 case CEED_RESTRICTION_ORIENTED: {
76 const bool *orients = NULL;
77 const CeedInt *offsets = NULL;
78
79 CeedCallBackend(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets));
80 CeedCallBackend(CeedElemRestrictionGetOrientations(rstr, CEED_MEM_HOST, &orients));
81 CeedCallBackend(CeedElemRestrictionCreateBlockedOriented(ceed_rstr, num_elem, elem_size, block_size, num_comp, comp_stride, l_size,
82 CEED_MEM_HOST, CEED_COPY_VALUES, offsets, orients, &block_rstr[i + start_e]));
83 CeedCallBackend(CeedElemRestrictionRestoreOffsets(rstr, &offsets));
84 CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr, &orients));
85 } break;
86 case CEED_RESTRICTION_CURL_ORIENTED: {
87 const CeedInt8 *curl_orients = NULL;
88 const CeedInt *offsets = NULL;
89
90 CeedCallBackend(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets));
91 CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr, CEED_MEM_HOST, &curl_orients));
92 CeedCallBackend(CeedElemRestrictionCreateBlockedCurlOriented(ceed_rstr, num_elem, elem_size, block_size, num_comp, comp_stride, l_size,
93 CEED_MEM_HOST, CEED_COPY_VALUES, offsets, curl_orients,
94 &block_rstr[i + start_e]));
95 CeedCallBackend(CeedElemRestrictionRestoreOffsets(rstr, &offsets));
96 CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr, &curl_orients));
97 } break;
98 case CEED_RESTRICTION_STRIDED: {
99 CeedInt strides[3];
100
101 CeedCallBackend(CeedElemRestrictionGetStrides(rstr, strides));
102 CeedCallBackend(CeedElemRestrictionCreateBlockedStrided(ceed_rstr, num_elem, elem_size, block_size, num_comp, l_size, strides,
103 &block_rstr[i + start_e]));
104 } break;
105 case CEED_RESTRICTION_POINTS:
106 // Empty case - won't occur
107 break;
108 }
109 CeedCallBackend(CeedDestroy(&ceed_rstr));
110 CeedCallBackend(CeedElemRestrictionDestroy(&rstr));
111 CeedCallBackend(CeedElemRestrictionCreateVector(block_rstr[i + start_e], NULL, &e_vecs_full[i + start_e]));
112 }
113
114 switch (eval_mode) {
115 case CEED_EVAL_NONE:
116 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
117 q_size = (CeedSize)Q * size * block_size;
118 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
119 break;
120 case CEED_EVAL_INTERP:
121 case CEED_EVAL_GRAD:
122 case CEED_EVAL_DIV:
123 case CEED_EVAL_CURL:
124 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
125 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
126 CeedCallBackend(CeedBasisGetNumNodes(basis, &P));
127 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
128 CeedCallBackend(CeedBasisDestroy(&basis));
129 e_size = (CeedSize)P * num_comp * block_size;
130 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i]));
131 q_size = (CeedSize)Q * size * block_size;
132 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
133 break;
134 case CEED_EVAL_WEIGHT: // Only on input fields
135 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
136 q_size = (CeedSize)Q * block_size;
137 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
138 CeedCallBackend(CeedBasisApply(basis, block_size, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i]));
139 CeedCallBackend(CeedBasisDestroy(&basis));
140 break;
141 }
142 }
143 // Drop duplicate restrictions
144 if (is_input) {
145 for (CeedInt i = 0; i < num_fields; i++) {
146 CeedVector vec_i;
147 CeedElemRestriction rstr_i;
148
149 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i));
150 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i));
151 for (CeedInt j = i + 1; j < num_fields; j++) {
152 CeedVector vec_j;
153 CeedElemRestriction rstr_j;
154
155 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j));
156 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j));
157 if (vec_i == vec_j && rstr_i == rstr_j) {
158 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j]));
159 CeedCallBackend(CeedVectorReferenceCopy(e_vecs_full[i + start_e], &e_vecs_full[j + start_e]));
160 skip_rstr[j] = true;
161 }
162 CeedCallBackend(CeedVectorDestroy(&vec_j));
163 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j));
164 }
165 CeedCallBackend(CeedVectorDestroy(&vec_i));
166 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i));
167 }
168 } else {
169 for (CeedInt i = num_fields - 1; i >= 0; i--) {
170 CeedVector vec_i;
171 CeedElemRestriction rstr_i;
172
173 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i));
174 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i));
175 for (CeedInt j = i - 1; j >= 0; j--) {
176 CeedVector vec_j;
177 CeedElemRestriction rstr_j;
178
179 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j));
180 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j));
181 if (vec_i == vec_j && rstr_i == rstr_j) {
182 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j]));
183 CeedCallBackend(CeedVectorReferenceCopy(e_vecs_full[i + start_e], &e_vecs_full[j + start_e]));
184 skip_rstr[j] = true;
185 apply_add_basis[i] = true;
186 e_data_out_indices[j] = i;
187 }
188 CeedCallBackend(CeedVectorDestroy(&vec_j));
189 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j));
190 }
191 CeedCallBackend(CeedVectorDestroy(&vec_i));
192 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i));
193 }
194 }
195 CeedCallBackend(CeedDestroy(&ceed));
196 return CEED_ERROR_SUCCESS;
197 }
198
199 //------------------------------------------------------------------------------
200 // Setup Operator
201 //------------------------------------------------------------------------------
CeedOperatorSetup_Blocked(CeedOperator op)202 static int CeedOperatorSetup_Blocked(CeedOperator op) {
203 bool is_setup_done;
204 CeedInt Q, num_input_fields, num_output_fields;
205 const CeedInt block_size = 8;
206 CeedQFunctionField *qf_input_fields, *qf_output_fields;
207 CeedQFunction qf;
208 CeedOperatorField *op_input_fields, *op_output_fields;
209 CeedOperator_Blocked *impl;
210
211 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done));
212 if (is_setup_done) return CEED_ERROR_SUCCESS;
213
214 CeedCallBackend(CeedOperatorGetData(op, &impl));
215 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
216 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q));
217 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf));
218 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
219 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
220
221 // Allocate
222 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->block_rstr));
223 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full));
224
225 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_in));
226 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_out));
227 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_data_out_indices));
228 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->apply_add_basis_out));
229 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states));
230 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in));
231 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out));
232 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in));
233 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out));
234
235 impl->num_inputs = num_input_fields;
236 impl->num_outputs = num_output_fields;
237
238 // Set up infield and outfield pointer arrays
239 // Infields
240 CeedCallBackend(CeedOperatorSetupFields_Blocked(qf, op, true, impl->skip_rstr_in, NULL, NULL, block_size, impl->block_rstr, impl->e_vecs_full,
241 impl->e_vecs_in, impl->q_vecs_in, 0, num_input_fields, Q));
242 // Outfields
243 CeedCallBackend(CeedOperatorSetupFields_Blocked(qf, op, false, impl->skip_rstr_out, impl->e_data_out_indices, impl->apply_add_basis_out, block_size,
244 impl->block_rstr, impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields,
245 num_output_fields, Q));
246
247 // Identity QFunctions
248 if (impl->is_identity_qf) {
249 CeedEvalMode in_mode, out_mode;
250 CeedQFunctionField *in_fields, *out_fields;
251
252 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &in_fields, NULL, &out_fields));
253 CeedCallBackend(CeedQFunctionFieldGetEvalMode(in_fields[0], &in_mode));
254 CeedCallBackend(CeedQFunctionFieldGetEvalMode(out_fields[0], &out_mode));
255
256 if (in_mode == CEED_EVAL_NONE && out_mode == CEED_EVAL_NONE) {
257 impl->is_identity_rstr_op = true;
258 } else {
259 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0]));
260 }
261 }
262
263 CeedCallBackend(CeedOperatorSetSetupDone(op));
264 CeedCallBackend(CeedQFunctionDestroy(&qf));
265 return CEED_ERROR_SUCCESS;
266 }
267
268 //------------------------------------------------------------------------------
269 // Setup Operator Inputs
270 //------------------------------------------------------------------------------
CeedOperatorSetupInputs_Blocked(CeedInt num_input_fields,CeedQFunctionField * qf_input_fields,CeedOperatorField * op_input_fields,CeedVector in_vec,bool skip_active,CeedScalar * e_data_full[2* CEED_FIELD_MAX],CeedOperator_Blocked * impl,CeedRequest * request)271 static inline int CeedOperatorSetupInputs_Blocked(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
272 CeedVector in_vec, bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX],
273 CeedOperator_Blocked *impl, CeedRequest *request) {
274 for (CeedInt i = 0; i < num_input_fields; i++) {
275 bool is_active;
276 uint64_t state;
277 CeedEvalMode eval_mode;
278 CeedVector vec;
279
280 // Get input vector
281 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
282 is_active = vec == CEED_VECTOR_ACTIVE;
283 if (is_active) {
284 if (skip_active) continue;
285 else vec = in_vec;
286 }
287
288 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
289 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip
290 } else {
291 // Restrict
292 CeedCallBackend(CeedVectorGetState(vec, &state));
293 if ((state != impl->input_states[i] || vec == in_vec) && !impl->skip_rstr_in[i]) {
294 CeedCallBackend(CeedElemRestrictionApply(impl->block_rstr[i], CEED_NOTRANSPOSE, vec, impl->e_vecs_full[i], request));
295 }
296 impl->input_states[i] = state;
297 // Get evec
298 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_full[i], CEED_MEM_HOST, (const CeedScalar **)&e_data_full[i]));
299 }
300 if (!is_active) CeedCallBackend(CeedVectorDestroy(&vec));
301 }
302 return CEED_ERROR_SUCCESS;
303 }
304
305 //------------------------------------------------------------------------------
306 // Input Basis Action
307 //------------------------------------------------------------------------------
CeedOperatorInputBasis_Blocked(CeedInt e,CeedInt Q,CeedQFunctionField * qf_input_fields,CeedOperatorField * op_input_fields,CeedInt num_input_fields,CeedInt block_size,bool skip_active,CeedScalar * e_data_full[2* CEED_FIELD_MAX],CeedOperator_Blocked * impl)308 static inline int CeedOperatorInputBasis_Blocked(CeedInt e, CeedInt Q, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
309 CeedInt num_input_fields, CeedInt block_size, bool skip_active,
310 CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Blocked *impl) {
311 for (CeedInt i = 0; i < num_input_fields; i++) {
312 CeedInt elem_size, size, num_comp;
313 CeedEvalMode eval_mode;
314 CeedElemRestriction elem_rstr;
315 CeedBasis basis;
316
317 // Skip active input
318 if (skip_active) {
319 bool is_active;
320 CeedVector vec;
321
322 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
323 is_active = vec == CEED_VECTOR_ACTIVE;
324 CeedCallBackend(CeedVectorDestroy(&vec));
325 if (is_active) continue;
326 }
327
328 // Get elem_size, eval_mode, size
329 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
330 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
331 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
332 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
333 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size));
334 // Basis action
335 switch (eval_mode) {
336 case CEED_EVAL_NONE:
337 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][(CeedSize)e * Q * size]));
338 break;
339 case CEED_EVAL_INTERP:
340 case CEED_EVAL_GRAD:
341 case CEED_EVAL_DIV:
342 case CEED_EVAL_CURL:
343 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
344 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
345 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][(CeedSize)e * elem_size * num_comp]));
346 CeedCallBackend(CeedBasisApply(basis, block_size, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i]));
347 CeedCallBackend(CeedBasisDestroy(&basis));
348 break;
349 case CEED_EVAL_WEIGHT:
350 break; // No action
351 }
352 }
353 return CEED_ERROR_SUCCESS;
354 }
355
356 //------------------------------------------------------------------------------
357 // Output Basis Action
358 //------------------------------------------------------------------------------
CeedOperatorOutputBasis_Blocked(CeedInt e,CeedInt Q,CeedQFunctionField * qf_output_fields,CeedOperatorField * op_output_fields,CeedInt block_size,CeedInt num_input_fields,CeedInt num_output_fields,bool * apply_add_basis,CeedOperator op,CeedScalar * e_data_full[2* CEED_FIELD_MAX],CeedOperator_Blocked * impl)359 static inline int CeedOperatorOutputBasis_Blocked(CeedInt e, CeedInt Q, CeedQFunctionField *qf_output_fields, CeedOperatorField *op_output_fields,
360 CeedInt block_size, CeedInt num_input_fields, CeedInt num_output_fields, bool *apply_add_basis,
361 CeedOperator op, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Blocked *impl) {
362 for (CeedInt i = 0; i < num_output_fields; i++) {
363 CeedInt elem_size, num_comp;
364 CeedEvalMode eval_mode;
365 CeedElemRestriction elem_rstr;
366 CeedBasis basis;
367
368 // Get elem_size, eval_mode, size
369 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
370 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
371 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
372 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
373 // Basis action
374 switch (eval_mode) {
375 case CEED_EVAL_NONE:
376 break; // No action
377 case CEED_EVAL_INTERP:
378 case CEED_EVAL_GRAD:
379 case CEED_EVAL_DIV:
380 case CEED_EVAL_CURL:
381 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
382 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
383 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER,
384 &e_data_full[i + num_input_fields][(CeedSize)e * elem_size * num_comp]));
385 if (apply_add_basis[i]) {
386 CeedCallBackend(CeedBasisApplyAdd(basis, block_size, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i]));
387 } else {
388 CeedCallBackend(CeedBasisApply(basis, block_size, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i]));
389 }
390 CeedCallBackend(CeedBasisDestroy(&basis));
391 break;
392 // LCOV_EXCL_START
393 case CEED_EVAL_WEIGHT: {
394 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
395 // LCOV_EXCL_STOP
396 }
397 }
398 }
399 return CEED_ERROR_SUCCESS;
400 }
401
402 //------------------------------------------------------------------------------
403 // Restore Input Vectors
404 //------------------------------------------------------------------------------
CeedOperatorRestoreInputs_Blocked(CeedInt num_input_fields,CeedQFunctionField * qf_input_fields,CeedOperatorField * op_input_fields,bool skip_active,CeedScalar * e_data_full[2* CEED_FIELD_MAX],CeedOperator_Blocked * impl)405 static inline int CeedOperatorRestoreInputs_Blocked(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
406 bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Blocked *impl) {
407 for (CeedInt i = 0; i < num_input_fields; i++) {
408 CeedEvalMode eval_mode;
409
410 // Skip active inputs
411 if (skip_active) {
412 bool is_active;
413 CeedVector vec;
414
415 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
416 is_active = vec == CEED_VECTOR_ACTIVE;
417 CeedCallBackend(CeedVectorDestroy(&vec));
418 if (is_active) continue;
419 }
420 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
421 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip
422 } else {
423 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_full[i], (const CeedScalar **)&e_data_full[i]));
424 }
425 }
426 return CEED_ERROR_SUCCESS;
427 }
428
429 //------------------------------------------------------------------------------
430 // Operator Apply
431 //------------------------------------------------------------------------------
CeedOperatorApplyAdd_Blocked(CeedOperator op,CeedVector in_vec,CeedVector out_vec,CeedRequest * request)432 static int CeedOperatorApplyAdd_Blocked(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) {
433 CeedInt Q, num_input_fields, num_output_fields, num_elem, size;
434 const CeedInt block_size = 8;
435 CeedEvalMode eval_mode;
436 CeedScalar *e_data_full[2 * CEED_FIELD_MAX] = {0};
437 CeedQFunctionField *qf_input_fields, *qf_output_fields;
438 CeedQFunction qf;
439 CeedOperatorField *op_input_fields, *op_output_fields;
440 CeedOperator_Blocked *impl;
441
442 // Setup
443 CeedCallBackend(CeedOperatorSetup_Blocked(op));
444
445 CeedCallBackend(CeedOperatorGetData(op, &impl));
446
447 // Restriction only operator
448 if (impl->is_identity_rstr_op) {
449 CeedCallBackend(CeedElemRestrictionApply(impl->block_rstr[0], CEED_NOTRANSPOSE, in_vec, impl->e_vecs_full[0], request));
450 CeedCallBackend(CeedElemRestrictionApply(impl->block_rstr[1], CEED_TRANSPOSE, impl->e_vecs_full[0], out_vec, request));
451 return CEED_ERROR_SUCCESS;
452 }
453 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
454 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q));
455 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
456 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
457 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
458 const CeedInt num_blocks = (num_elem / block_size) + !!(num_elem % block_size);
459
460 // Input Evecs and Restriction
461 CeedCallBackend(CeedOperatorSetupInputs_Blocked(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data_full, impl, request));
462
463 // Output Evecs
464 for (CeedInt i = num_output_fields - 1; i >= 0; i--) {
465 if (impl->skip_rstr_out[i]) {
466 e_data_full[i + num_input_fields] = e_data_full[impl->e_data_out_indices[i] + num_input_fields];
467 } else {
468 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs_full[i + impl->num_inputs], CEED_MEM_HOST, &e_data_full[i + num_input_fields]));
469 }
470 }
471
472 // Loop through elements
473 for (CeedInt e = 0; e < num_blocks * block_size; e += block_size) {
474 // Output pointers
475 for (CeedInt i = 0; i < num_output_fields; i++) {
476 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
477 if (eval_mode == CEED_EVAL_NONE) {
478 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size));
479 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER,
480 &e_data_full[i + num_input_fields][(CeedSize)e * Q * size]));
481 }
482 }
483
484 // Input basis apply
485 CeedCallBackend(CeedOperatorInputBasis_Blocked(e, Q, qf_input_fields, op_input_fields, num_input_fields, block_size, false, e_data_full, impl));
486
487 // Q function
488 if (!impl->is_identity_qf) {
489 CeedCallBackend(CeedQFunctionApply(qf, Q * block_size, impl->q_vecs_in, impl->q_vecs_out));
490 }
491
492 // Output basis apply
493 CeedCallBackend(CeedOperatorOutputBasis_Blocked(e, Q, qf_output_fields, op_output_fields, block_size, num_input_fields, num_output_fields,
494 impl->apply_add_basis_out, op, e_data_full, impl));
495 }
496
497 // Output restriction
498 for (CeedInt i = 0; i < num_output_fields; i++) {
499 bool is_active;
500 CeedVector vec;
501
502 if (impl->skip_rstr_out[i]) continue;
503 // Restore evec
504 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_full[i + impl->num_inputs], &e_data_full[i + num_input_fields]));
505 // Get output vector
506 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
507 is_active = vec == CEED_VECTOR_ACTIVE;
508 // Active
509 if (is_active) vec = out_vec;
510 // Restrict
511 CeedCallBackend(CeedElemRestrictionApply(impl->block_rstr[i + impl->num_inputs], CEED_TRANSPOSE, impl->e_vecs_full[i + impl->num_inputs], vec,
512 request));
513 if (!is_active) CeedCallBackend(CeedVectorDestroy(&vec));
514 }
515
516 // Restore input arrays
517 CeedCallBackend(CeedOperatorRestoreInputs_Blocked(num_input_fields, qf_input_fields, op_input_fields, false, e_data_full, impl));
518 CeedCallBackend(CeedQFunctionDestroy(&qf));
519 return CEED_ERROR_SUCCESS;
520 }
521
522 //------------------------------------------------------------------------------
523 // Core code for assembling linear QFunction
524 //------------------------------------------------------------------------------
CeedOperatorLinearAssembleQFunctionCore_Blocked(CeedOperator op,bool build_objects,CeedVector * assembled,CeedElemRestriction * rstr,CeedRequest * request)525 static inline int CeedOperatorLinearAssembleQFunctionCore_Blocked(CeedOperator op, bool build_objects, CeedVector *assembled,
526 CeedElemRestriction *rstr, CeedRequest *request) {
527 Ceed ceed;
528 CeedInt qf_size_in, qf_size_out, Q, num_input_fields, num_output_fields, num_elem;
529 const CeedInt block_size = 8;
530 CeedScalar *l_vec_array;
531 CeedScalar *e_data_full[2 * CEED_FIELD_MAX] = {0};
532 CeedQFunctionField *qf_input_fields, *qf_output_fields;
533 CeedQFunction qf;
534 CeedOperatorField *op_input_fields, *op_output_fields;
535 CeedOperator_Blocked *impl;
536
537 CeedCallBackend(CeedOperatorGetData(op, &impl));
538 qf_size_in = impl->qf_size_in;
539 qf_size_out = impl->qf_size_out;
540 CeedVector l_vec = impl->qf_l_vec;
541 CeedElemRestriction block_rstr = impl->qf_block_rstr;
542
543 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
544 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q));
545 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
546 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
547 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
548 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
549 const CeedInt num_blocks = (num_elem / block_size) + !!(num_elem % block_size);
550
551 // Setup
552 CeedCallBackend(CeedOperatorSetup_Blocked(op));
553
554 // Check for restriction only operator
555 CeedCheck(!impl->is_identity_rstr_op, ceed, CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported");
556
557 // Input Evecs and Restriction
558 CeedCallBackend(CeedOperatorSetupInputs_Blocked(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request));
559
560 // Count number of active input fields
561 if (qf_size_in == 0) {
562 for (CeedInt i = 0; i < num_input_fields; i++) {
563 CeedInt field_size;
564 CeedVector vec;
565
566 // Check if active input
567 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
568 if (vec == CEED_VECTOR_ACTIVE) {
569 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size));
570 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0));
571 qf_size_in += field_size;
572 }
573 CeedCallBackend(CeedVectorDestroy(&vec));
574 }
575 CeedCheck(qf_size_in > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs");
576 impl->qf_size_in = qf_size_in;
577 }
578
579 // Count number of active output fields
580 if (qf_size_out == 0) {
581 for (CeedInt i = 0; i < num_output_fields; i++) {
582 CeedInt field_size;
583 CeedVector vec;
584
585 // Check if active output
586 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
587 if (vec == CEED_VECTOR_ACTIVE) {
588 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size));
589 qf_size_out += field_size;
590 }
591 CeedCallBackend(CeedVectorDestroy(&vec));
592 }
593 CeedCheck(qf_size_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs");
594 impl->qf_size_out = qf_size_out;
595 }
596
597 // Setup Lvec
598 if (!l_vec) {
599 const CeedSize l_size = (CeedSize)num_blocks * block_size * Q * qf_size_in * qf_size_out;
600
601 CeedCallBackend(CeedVectorCreate(ceed, l_size, &l_vec));
602 impl->qf_l_vec = l_vec;
603 }
604 CeedCallBackend(CeedVectorGetArrayWrite(l_vec, CEED_MEM_HOST, &l_vec_array));
605
606 // Setup block restriction
607 if (!block_rstr) {
608 const CeedInt strides[3] = {1, Q, qf_size_in * qf_size_out * Q};
609
610 CeedCallBackend(CeedElemRestrictionCreateBlockedStrided(ceed, num_elem, Q, block_size, qf_size_in * qf_size_out,
611 qf_size_in * qf_size_out * num_elem * Q, strides, &block_rstr));
612 impl->qf_block_rstr = block_rstr;
613 }
614
615 // Build objects if needed
616 if (build_objects) {
617 const CeedSize l_size = (CeedSize)num_elem * Q * qf_size_in * qf_size_out;
618 const CeedInt strides[3] = {1, Q, qf_size_in * qf_size_out * Q};
619
620 // Create output restriction
621 CeedCallBackend(CeedElemRestrictionCreateStrided(ceed, num_elem, Q, qf_size_in * qf_size_out,
622 (CeedSize)qf_size_in * (CeedSize)qf_size_out * (CeedSize)num_elem * (CeedSize)Q, strides, rstr));
623 // Create assembled vector
624 CeedCallBackend(CeedVectorCreate(ceed, l_size, assembled));
625 }
626
627 // Loop through elements
628 for (CeedInt e = 0; e < num_blocks * block_size; e += block_size) {
629 // Input basis apply
630 CeedCallBackend(CeedOperatorInputBasis_Blocked(e, Q, qf_input_fields, op_input_fields, num_input_fields, block_size, true, e_data_full, impl));
631
632 // Assemble QFunction
633 for (CeedInt i = 0; i < num_input_fields; i++) {
634 bool is_active;
635 CeedInt field_size;
636 CeedVector vec;
637
638 // Check if active input
639 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
640 is_active = vec == CEED_VECTOR_ACTIVE;
641 CeedCallBackend(CeedVectorDestroy(&vec));
642 if (!is_active) continue;
643 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size));
644 for (CeedInt field = 0; field < field_size; field++) {
645 // Set current portion of input to 1.0
646 {
647 CeedScalar *array;
648
649 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array));
650 for (CeedInt j = 0; j < Q * block_size; j++) array[field * Q * block_size + j] = 1.0;
651 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array));
652 }
653
654 if (!impl->is_identity_qf) {
655 // Set Outputs
656 for (CeedInt out = 0; out < num_output_fields; out++) {
657 CeedInt field_size;
658 CeedVector vec;
659
660 // Get output vector
661 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec));
662 // Check if active output
663 if (vec == CEED_VECTOR_ACTIVE) {
664 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, l_vec_array));
665 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &field_size));
666 l_vec_array += field_size * Q * block_size; // Advance the pointer by the size of the output
667 }
668 CeedCallBackend(CeedVectorDestroy(&vec));
669 }
670 // Apply QFunction
671 CeedCallBackend(CeedQFunctionApply(qf, Q * block_size, impl->q_vecs_in, impl->q_vecs_out));
672 } else {
673 CeedInt field_size;
674 const CeedScalar *array;
675
676 // Copy Identity Outputs
677 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &field_size));
678 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &array));
679 for (CeedInt j = 0; j < field_size * Q * block_size; j++) l_vec_array[j] = array[j];
680 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &array));
681 l_vec_array += field_size * Q * block_size;
682 }
683 // Reset input to 0.0
684 {
685 CeedScalar *array;
686
687 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array));
688 for (CeedInt j = 0; j < Q * block_size; j++) array[field * Q * block_size + j] = 0.0;
689 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array));
690 }
691 }
692 }
693 }
694
695 // Un-set output Qvecs to prevent accidental overwrite of Assembled
696 if (!impl->is_identity_qf) {
697 for (CeedInt out = 0; out < num_output_fields; out++) {
698 CeedVector vec;
699
700 // Check if active output
701 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec));
702 if (vec == CEED_VECTOR_ACTIVE) {
703 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL));
704 }
705 CeedCallBackend(CeedVectorDestroy(&vec));
706 }
707 }
708
709 // Restore input arrays
710 CeedCallBackend(CeedOperatorRestoreInputs_Blocked(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl));
711
712 // Output blocked restriction
713 CeedCallBackend(CeedVectorRestoreArray(l_vec, &l_vec_array));
714 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0));
715 CeedCallBackend(CeedElemRestrictionApply(block_rstr, CEED_TRANSPOSE, l_vec, *assembled, request));
716 CeedCallBackend(CeedDestroy(&ceed));
717 CeedCallBackend(CeedQFunctionDestroy(&qf));
718 return CEED_ERROR_SUCCESS;
719 }
720
721 //------------------------------------------------------------------------------
722 // Assemble Linear QFunction
723 //------------------------------------------------------------------------------
CeedOperatorLinearAssembleQFunction_Blocked(CeedOperator op,CeedVector * assembled,CeedElemRestriction * rstr,CeedRequest * request)724 static int CeedOperatorLinearAssembleQFunction_Blocked(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) {
725 return CeedOperatorLinearAssembleQFunctionCore_Blocked(op, true, assembled, rstr, request);
726 }
727
728 //------------------------------------------------------------------------------
729 // Update Assembled Linear QFunction
730 //------------------------------------------------------------------------------
CeedOperatorLinearAssembleQFunctionUpdate_Blocked(CeedOperator op,CeedVector assembled,CeedElemRestriction rstr,CeedRequest * request)731 static int CeedOperatorLinearAssembleQFunctionUpdate_Blocked(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) {
732 return CeedOperatorLinearAssembleQFunctionCore_Blocked(op, false, &assembled, &rstr, request);
733 }
734
735 //------------------------------------------------------------------------------
736 // Operator Destroy
737 //------------------------------------------------------------------------------
CeedOperatorDestroy_Blocked(CeedOperator op)738 static int CeedOperatorDestroy_Blocked(CeedOperator op) {
739 CeedOperator_Blocked *impl;
740
741 CeedCallBackend(CeedOperatorGetData(op, &impl));
742
743 CeedCallBackend(CeedFree(&impl->skip_rstr_in));
744 CeedCallBackend(CeedFree(&impl->skip_rstr_out));
745 CeedCallBackend(CeedFree(&impl->e_data_out_indices));
746 CeedCallBackend(CeedFree(&impl->apply_add_basis_out));
747 for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) {
748 CeedCallBackend(CeedElemRestrictionDestroy(&impl->block_rstr[i]));
749 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_full[i]));
750 }
751 CeedCallBackend(CeedFree(&impl->block_rstr));
752 CeedCallBackend(CeedFree(&impl->e_vecs_full));
753 CeedCallBackend(CeedFree(&impl->input_states));
754
755 for (CeedInt i = 0; i < impl->num_inputs; i++) {
756 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_in[i]));
757 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i]));
758 }
759 CeedCallBackend(CeedFree(&impl->e_vecs_in));
760 CeedCallBackend(CeedFree(&impl->q_vecs_in));
761
762 for (CeedInt i = 0; i < impl->num_outputs; i++) {
763 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_out[i]));
764 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i]));
765 }
766 CeedCallBackend(CeedFree(&impl->e_vecs_out));
767 CeedCallBackend(CeedFree(&impl->q_vecs_out));
768
769 // QFunction assembly data
770 CeedCallBackend(CeedVectorDestroy(&impl->qf_l_vec));
771 CeedCallBackend(CeedElemRestrictionDestroy(&impl->qf_block_rstr));
772
773 CeedCallBackend(CeedFree(&impl));
774 return CEED_ERROR_SUCCESS;
775 }
776
777 //------------------------------------------------------------------------------
778 // Operator Create
779 //------------------------------------------------------------------------------
CeedOperatorCreate_Blocked(CeedOperator op)780 int CeedOperatorCreate_Blocked(CeedOperator op) {
781 Ceed ceed;
782 CeedOperator_Blocked *impl;
783
784 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
785 CeedCallBackend(CeedCalloc(1, &impl));
786 CeedCallBackend(CeedOperatorSetData(op, impl));
787 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Blocked));
788 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Blocked));
789 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Blocked));
790 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Blocked));
791 CeedCallBackend(CeedDestroy(&ceed));
792 return CEED_ERROR_SUCCESS;
793 }
794
795 //------------------------------------------------------------------------------
796