xref: /libCEED/backends/hip-ref/ceed-hip-ref-operator.c (revision 24a65d3da2f623912f26b42c0b9ba6f37de25307)
1 // Copyright (c) 2017-2022, 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 <ceed/jit-tools.h>
11 #include <assert.h>
12 #include <stdbool.h>
13 #include <string.h>
14 #include <hip/hip_runtime.h>
15 
16 #include "../hip/ceed-hip-common.h"
17 #include "../hip/ceed-hip-compile.h"
18 #include "ceed-hip-ref.h"
19 
20 //------------------------------------------------------------------------------
21 // Destroy operator
22 //------------------------------------------------------------------------------
23 static int CeedOperatorDestroy_Hip(CeedOperator op) {
24   CeedOperator_Hip *impl;
25 
26   CeedCallBackend(CeedOperatorGetData(op, &impl));
27 
28   // Apply data
29   for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) {
30     CeedCallBackend(CeedVectorDestroy(&impl->e_vecs[i]));
31   }
32   CeedCallBackend(CeedFree(&impl->e_vecs));
33 
34   for (CeedInt i = 0; i < impl->num_inputs; i++) {
35     CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i]));
36   }
37   CeedCallBackend(CeedFree(&impl->q_vecs_in));
38 
39   for (CeedInt i = 0; i < impl->num_outputs; i++) {
40     CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i]));
41   }
42   CeedCallBackend(CeedFree(&impl->q_vecs_out));
43 
44   // QFunction assembly data
45   for (CeedInt i = 0; i < impl->num_active_in; i++) {
46     CeedCallBackend(CeedVectorDestroy(&impl->qf_active_in[i]));
47   }
48   CeedCallBackend(CeedFree(&impl->qf_active_in));
49 
50   // Diag data
51   if (impl->diag) {
52     Ceed ceed;
53 
54     CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
55     if (impl->diag->module) {
56       CeedCallHip(ceed, hipModuleUnload(impl->diag->module));
57     }
58     if (impl->diag->module_point_block) {
59       CeedCallHip(ceed, hipModuleUnload(impl->diag->module_point_block));
60     }
61     CeedCallHip(ceed, hipFree(impl->diag->d_eval_modes_in));
62     CeedCallHip(ceed, hipFree(impl->diag->d_eval_modes_out));
63     CeedCallHip(ceed, hipFree(impl->diag->d_identity));
64     CeedCallHip(ceed, hipFree(impl->diag->d_interp_in));
65     CeedCallHip(ceed, hipFree(impl->diag->d_interp_out));
66     CeedCallHip(ceed, hipFree(impl->diag->d_grad_in));
67     CeedCallHip(ceed, hipFree(impl->diag->d_grad_out));
68     CeedCallHip(ceed, hipFree(impl->diag->d_div_in));
69     CeedCallHip(ceed, hipFree(impl->diag->d_div_out));
70     CeedCallHip(ceed, hipFree(impl->diag->d_curl_in));
71     CeedCallHip(ceed, hipFree(impl->diag->d_curl_out));
72     CeedCallBackend(CeedElemRestrictionDestroy(&impl->diag->diag_rstr));
73     CeedCallBackend(CeedElemRestrictionDestroy(&impl->diag->point_block_diag_rstr));
74     CeedCallBackend(CeedVectorDestroy(&impl->diag->elem_diag));
75     CeedCallBackend(CeedVectorDestroy(&impl->diag->point_block_elem_diag));
76   }
77   CeedCallBackend(CeedFree(&impl->diag));
78 
79   if (impl->asmb) {
80     Ceed ceed;
81 
82     CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
83     CeedCallHip(ceed, hipModuleUnload(impl->asmb->module));
84     CeedCallHip(ceed, hipFree(impl->asmb->d_B_in));
85     CeedCallHip(ceed, hipFree(impl->asmb->d_B_out));
86   }
87   CeedCallBackend(CeedFree(&impl->asmb));
88 
89   CeedCallBackend(CeedFree(&impl));
90   return CEED_ERROR_SUCCESS;
91 }
92 
93 //------------------------------------------------------------------------------
94 // Setup infields or outfields
95 //------------------------------------------------------------------------------
96 static int CeedOperatorSetupFields_Hip(CeedQFunction qf, CeedOperator op, bool is_input, CeedVector *e_vecs, CeedVector *q_vecs, CeedInt start_e,
97                                        CeedInt num_fields, CeedInt Q, CeedInt num_elem) {
98   Ceed                ceed;
99   CeedQFunctionField *qf_fields;
100   CeedOperatorField  *op_fields;
101 
102   CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
103   if (is_input) {
104     CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL));
105     CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL));
106   } else {
107     CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields));
108     CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields));
109   }
110 
111   // Loop over fields
112   for (CeedInt i = 0; i < num_fields; i++) {
113     bool         is_strided = false, skip_restriction = false;
114     CeedSize     q_size;
115     CeedInt      size;
116     CeedEvalMode eval_mode;
117     CeedBasis    basis;
118 
119     CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
120     if (eval_mode != CEED_EVAL_WEIGHT) {
121       CeedElemRestriction elem_rstr;
122 
123       // Check whether this field can skip the element restriction:
124       // Must be passive input, with eval_mode NONE, and have a strided restriction with CEED_STRIDES_BACKEND.
125       CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr));
126 
127       // First, check whether the field is input or output:
128       if (is_input) {
129         CeedVector vec;
130 
131         // Check for passive input
132         CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec));
133         if (vec != CEED_VECTOR_ACTIVE) {
134           // Check eval_mode
135           if (eval_mode == CEED_EVAL_NONE) {
136             // Check for strided restriction
137             CeedCallBackend(CeedElemRestrictionIsStrided(elem_rstr, &is_strided));
138             if (is_strided) {
139               // Check if vector is already in preferred backend ordering
140               CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &skip_restriction));
141             }
142           }
143         }
144       }
145       if (skip_restriction) {
146         // We do not need an E-Vector, but will use the input field vector's data directly in the operator application.
147         e_vecs[i + start_e] = NULL;
148       } else {
149         CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs[i + start_e]));
150       }
151     }
152 
153     switch (eval_mode) {
154       case CEED_EVAL_NONE:
155         CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
156         q_size = (CeedSize)num_elem * Q * size;
157         CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
158         break;
159       case CEED_EVAL_INTERP:
160       case CEED_EVAL_GRAD:
161       case CEED_EVAL_DIV:
162       case CEED_EVAL_CURL:
163         CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
164         q_size = (CeedSize)num_elem * Q * size;
165         CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
166         break;
167       case CEED_EVAL_WEIGHT:  // Only on input fields
168         CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
169         q_size = (CeedSize)num_elem * Q;
170         CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
171         CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i]));
172         break;
173     }
174   }
175   return CEED_ERROR_SUCCESS;
176 }
177 
178 //------------------------------------------------------------------------------
179 // CeedOperator needs to connect all the named fields (be they active or passive) to the named inputs and outputs of its CeedQFunction.
180 //------------------------------------------------------------------------------
181 static int CeedOperatorSetup_Hip(CeedOperator op) {
182   Ceed                ceed;
183   bool                is_setup_done;
184   CeedInt             Q, num_elem, num_input_fields, num_output_fields;
185   CeedQFunctionField *qf_input_fields, *qf_output_fields;
186   CeedQFunction       qf;
187   CeedOperatorField  *op_input_fields, *op_output_fields;
188   CeedOperator_Hip   *impl;
189 
190   CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done));
191   if (is_setup_done) return CEED_ERROR_SUCCESS;
192 
193   CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
194   CeedCallBackend(CeedOperatorGetData(op, &impl));
195   CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
196   CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q));
197   CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
198   CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
199   CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
200 
201   // Allocate
202   CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs));
203   CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in));
204   CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out));
205   impl->num_inputs  = num_input_fields;
206   impl->num_outputs = num_output_fields;
207 
208   // Set up infield and outfield e_vecs and q_vecs
209   // Infields
210   CeedCallBackend(CeedOperatorSetupFields_Hip(qf, op, true, impl->e_vecs, impl->q_vecs_in, 0, num_input_fields, Q, num_elem));
211   // Outfields
212   CeedCallBackend(CeedOperatorSetupFields_Hip(qf, op, false, impl->e_vecs, impl->q_vecs_out, num_input_fields, num_output_fields, Q, num_elem));
213 
214   CeedCallBackend(CeedOperatorSetSetupDone(op));
215   return CEED_ERROR_SUCCESS;
216 }
217 
218 //------------------------------------------------------------------------------
219 // Setup Operator Inputs
220 //------------------------------------------------------------------------------
221 static inline int CeedOperatorSetupInputs_Hip(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
222                                               CeedVector in_vec, const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX],
223                                               CeedOperator_Hip *impl, CeedRequest *request) {
224   for (CeedInt i = 0; i < num_input_fields; i++) {
225     CeedEvalMode        eval_mode;
226     CeedVector          vec;
227     CeedElemRestriction elem_rstr;
228 
229     // Get input vector
230     CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
231     if (vec == CEED_VECTOR_ACTIVE) {
232       if (skip_active) continue;
233       else vec = in_vec;
234     }
235 
236     CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
237     if (eval_mode == CEED_EVAL_WEIGHT) {  // Skip
238     } else {
239       // Get input vector
240       CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
241       // Get input element restriction
242       CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
243       if (vec == CEED_VECTOR_ACTIVE) vec = in_vec;
244       // Restrict, if necessary
245       if (!impl->e_vecs[i]) {
246         // No restriction for this field; read data directly from vec.
247         CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, (const CeedScalar **)&e_data[i]));
248       } else {
249         CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, vec, impl->e_vecs[i], request));
250         // Get evec
251         CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs[i], CEED_MEM_DEVICE, (const CeedScalar **)&e_data[i]));
252       }
253     }
254   }
255   return CEED_ERROR_SUCCESS;
256 }
257 
258 //------------------------------------------------------------------------------
259 // Input Basis Action
260 //------------------------------------------------------------------------------
261 static inline int CeedOperatorInputBasis_Hip(CeedInt num_elem, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
262                                              CeedInt num_input_fields, const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX],
263                                              CeedOperator_Hip *impl) {
264   for (CeedInt i = 0; i < num_input_fields; i++) {
265     CeedInt             elem_size, size;
266     CeedEvalMode        eval_mode;
267     CeedElemRestriction elem_rstr;
268     CeedBasis           basis;
269 
270     // Skip active input
271     if (skip_active) {
272       CeedVector vec;
273 
274       CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
275       if (vec == CEED_VECTOR_ACTIVE) continue;
276     }
277     // Get elem_size, eval_mode, size
278     CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
279     CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
280     CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
281     CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size));
282     // Basis action
283     switch (eval_mode) {
284       case CEED_EVAL_NONE:
285         CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i]));
286         break;
287       case CEED_EVAL_INTERP:
288       case CEED_EVAL_GRAD:
289       case CEED_EVAL_DIV:
290       case CEED_EVAL_CURL:
291         CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
292         CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs[i], impl->q_vecs_in[i]));
293         break;
294       case CEED_EVAL_WEIGHT:
295         break;  // No action
296     }
297   }
298   return CEED_ERROR_SUCCESS;
299 }
300 
301 //------------------------------------------------------------------------------
302 // Restore Input Vectors
303 //------------------------------------------------------------------------------
304 static inline int CeedOperatorRestoreInputs_Hip(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
305                                                 const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], CeedOperator_Hip *impl) {
306   for (CeedInt i = 0; i < num_input_fields; i++) {
307     CeedEvalMode eval_mode;
308     CeedVector   vec;
309 
310     // Skip active input
311     if (skip_active) {
312       CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
313       if (vec == CEED_VECTOR_ACTIVE) continue;
314     }
315     CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
316     if (eval_mode == CEED_EVAL_WEIGHT) {  // Skip
317     } else {
318       if (!impl->e_vecs[i]) {  // This was a skip_restriction case
319         CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
320         CeedCallBackend(CeedVectorRestoreArrayRead(vec, (const CeedScalar **)&e_data[i]));
321       } else {
322         CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs[i], (const CeedScalar **)&e_data[i]));
323       }
324     }
325   }
326   return CEED_ERROR_SUCCESS;
327 }
328 
329 //------------------------------------------------------------------------------
330 // Apply and add to output
331 //------------------------------------------------------------------------------
332 static int CeedOperatorApplyAdd_Hip(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) {
333   CeedInt             Q, num_elem, elem_size, num_input_fields, num_output_fields, size;
334   CeedScalar         *e_data[2 * CEED_FIELD_MAX] = {NULL};
335   CeedQFunctionField *qf_input_fields, *qf_output_fields;
336   CeedQFunction       qf;
337   CeedOperatorField  *op_input_fields, *op_output_fields;
338   CeedOperator_Hip   *impl;
339 
340   CeedCallBackend(CeedOperatorGetData(op, &impl));
341   CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
342   CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q));
343   CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
344   CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
345   CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
346 
347   // Setup
348   CeedCallBackend(CeedOperatorSetup_Hip(op));
349 
350   // Input Evecs and Restriction
351   CeedCallBackend(CeedOperatorSetupInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data, impl, request));
352 
353   // Input basis apply if needed
354   CeedCallBackend(CeedOperatorInputBasis_Hip(num_elem, qf_input_fields, op_input_fields, num_input_fields, false, e_data, impl));
355 
356   // Output pointers, as necessary
357   for (CeedInt i = 0; i < num_output_fields; i++) {
358     CeedEvalMode eval_mode;
359 
360     CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
361     if (eval_mode == CEED_EVAL_NONE) {
362       // Set the output Q-Vector to use the E-Vector data directly.
363       CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs[i + impl->num_inputs], CEED_MEM_DEVICE, &e_data[i + num_input_fields]));
364       CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i + num_input_fields]));
365     }
366   }
367 
368   // Q function
369   CeedCallBackend(CeedQFunctionApply(qf, num_elem * Q, impl->q_vecs_in, impl->q_vecs_out));
370 
371   // Output basis apply if needed
372   for (CeedInt i = 0; i < num_output_fields; i++) {
373     CeedEvalMode        eval_mode;
374     CeedElemRestriction elem_rstr;
375     CeedBasis           basis;
376 
377     // Get elem_size, eval_mode, size
378     CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
379     CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
380     CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
381     CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size));
382     // Basis action
383     switch (eval_mode) {
384       case CEED_EVAL_NONE:
385         break;  // No action
386       case CEED_EVAL_INTERP:
387       case CEED_EVAL_GRAD:
388       case CEED_EVAL_DIV:
389       case CEED_EVAL_CURL:
390         CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
391         CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs[i + impl->num_inputs]));
392         break;
393       // LCOV_EXCL_START
394       case CEED_EVAL_WEIGHT: {
395         return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
396         // LCOV_EXCL_STOP
397       }
398     }
399   }
400 
401   // Output restriction
402   for (CeedInt i = 0; i < num_output_fields; i++) {
403     CeedEvalMode        eval_mode;
404     CeedVector          vec;
405     CeedElemRestriction elem_rstr;
406 
407     // Restore evec
408     CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
409     if (eval_mode == CEED_EVAL_NONE) {
410       CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs[i + impl->num_inputs], &e_data[i + num_input_fields]));
411     }
412     // Get output vector
413     CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
414     // Restrict
415     CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
416     // Active
417     if (vec == CEED_VECTOR_ACTIVE) vec = out_vec;
418 
419     CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs[i + impl->num_inputs], vec, request));
420   }
421 
422   // Restore input arrays
423   CeedCallBackend(CeedOperatorRestoreInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, false, e_data, impl));
424   return CEED_ERROR_SUCCESS;
425 }
426 
427 //------------------------------------------------------------------------------
428 // Linear QFunction Assembly Core
429 //------------------------------------------------------------------------------
430 static inline int CeedOperatorLinearAssembleQFunctionCore_Hip(CeedOperator op, bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr,
431                                                               CeedRequest *request) {
432   Ceed                ceed, ceed_parent;
433   CeedInt             num_active_in, num_active_out, Q, num_elem, num_input_fields, num_output_fields, size;
434   CeedScalar         *assembled_array, *e_data[2 * CEED_FIELD_MAX] = {NULL};
435   CeedVector         *active_inputs;
436   CeedQFunctionField *qf_input_fields, *qf_output_fields;
437   CeedQFunction       qf;
438   CeedOperatorField  *op_input_fields, *op_output_fields;
439   CeedOperator_Hip   *impl;
440 
441   CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
442   CeedCallBackend(CeedOperatorGetFallbackParentCeed(op, &ceed_parent));
443   CeedCallBackend(CeedOperatorGetData(op, &impl));
444   CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q));
445   CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
446   CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
447   CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
448   CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
449   active_inputs = impl->qf_active_in;
450   num_active_in = impl->num_active_in, num_active_out = impl->num_active_out;
451 
452   // Setup
453   CeedCallBackend(CeedOperatorSetup_Hip(op));
454 
455   // Input Evecs and Restriction
456   CeedCallBackend(CeedOperatorSetupInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request));
457 
458   // Count number of active input fields
459   if (!num_active_in) {
460     for (CeedInt i = 0; i < num_input_fields; i++) {
461       CeedScalar *q_vec_array;
462       CeedVector  vec;
463 
464       // Get input vector
465       CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
466       // Check if active input
467       if (vec == CEED_VECTOR_ACTIVE) {
468         CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size));
469         CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0));
470         CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, &q_vec_array));
471         CeedCallBackend(CeedRealloc(num_active_in + size, &active_inputs));
472         for (CeedInt field = 0; field < size; field++) {
473           CeedSize q_size = (CeedSize)Q * num_elem;
474 
475           CeedCallBackend(CeedVectorCreate(ceed, q_size, &active_inputs[num_active_in + field]));
476           CeedCallBackend(
477               CeedVectorSetArray(active_inputs[num_active_in + field], CEED_MEM_DEVICE, CEED_USE_POINTER, &q_vec_array[field * Q * num_elem]));
478         }
479         num_active_in += size;
480         CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &q_vec_array));
481       }
482     }
483     impl->num_active_in = num_active_in;
484     impl->qf_active_in  = active_inputs;
485   }
486 
487   // Count number of active output fields
488   if (!num_active_out) {
489     for (CeedInt i = 0; i < num_output_fields; i++) {
490       CeedVector vec;
491 
492       // Get output vector
493       CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
494       // Check if active output
495       if (vec == CEED_VECTOR_ACTIVE) {
496         CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size));
497         num_active_out += size;
498       }
499     }
500     impl->num_active_out = num_active_out;
501   }
502 
503   // Check sizes
504   CeedCheck(num_active_in > 0 && num_active_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs");
505 
506   // Build objects if needed
507   if (build_objects) {
508     CeedSize l_size     = (CeedSize)num_elem * Q * num_active_in * num_active_out;
509     CeedInt  strides[3] = {1, num_elem * Q, Q}; /* *NOPAD* */
510 
511     // Create output restriction
512     CeedCallBackend(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, num_active_in * num_active_out,
513                                                      num_active_in * num_active_out * num_elem * Q, strides, rstr));
514     // Create assembled vector
515     CeedCallBackend(CeedVectorCreate(ceed_parent, l_size, assembled));
516   }
517   CeedCallBackend(CeedVectorSetValue(*assembled, 0.0));
518   CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_DEVICE, &assembled_array));
519 
520   // Input basis apply
521   CeedCallBackend(CeedOperatorInputBasis_Hip(num_elem, qf_input_fields, op_input_fields, num_input_fields, true, e_data, impl));
522 
523   // Assemble QFunction
524   for (CeedInt in = 0; in < num_active_in; in++) {
525     // Set Inputs
526     CeedCallBackend(CeedVectorSetValue(active_inputs[in], 1.0));
527     if (num_active_in > 1) {
528       CeedCallBackend(CeedVectorSetValue(active_inputs[(in + num_active_in - 1) % num_active_in], 0.0));
529     }
530     // Set Outputs
531     for (CeedInt out = 0; out < num_output_fields; out++) {
532       CeedVector vec;
533 
534       // Get output vector
535       CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec));
536       // Check if active output
537       if (vec == CEED_VECTOR_ACTIVE) {
538         CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_DEVICE, CEED_USE_POINTER, assembled_array));
539         CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &size));
540         assembled_array += size * Q * num_elem;  // Advance the pointer by the size of the output
541       }
542     }
543     // Apply QFunction
544     CeedCallBackend(CeedQFunctionApply(qf, Q * num_elem, impl->q_vecs_in, impl->q_vecs_out));
545   }
546 
547   // Un-set output q_vecs to prevent accidental overwrite of Assembled
548   for (CeedInt out = 0; out < num_output_fields; out++) {
549     CeedVector vec;
550 
551     // Get output vector
552     CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec));
553     // Check if active output
554     if (vec == CEED_VECTOR_ACTIVE) {
555       CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_DEVICE, NULL));
556     }
557   }
558 
559   // Restore input arrays
560   CeedCallBackend(CeedOperatorRestoreInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl));
561 
562   // Restore output
563   CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array));
564   return CEED_ERROR_SUCCESS;
565 }
566 
567 //------------------------------------------------------------------------------
568 // Assemble Linear QFunction
569 //------------------------------------------------------------------------------
570 static int CeedOperatorLinearAssembleQFunction_Hip(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) {
571   return CeedOperatorLinearAssembleQFunctionCore_Hip(op, true, assembled, rstr, request);
572 }
573 
574 //------------------------------------------------------------------------------
575 // Update Assembled Linear QFunction
576 //------------------------------------------------------------------------------
577 static int CeedOperatorLinearAssembleQFunctionUpdate_Hip(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) {
578   return CeedOperatorLinearAssembleQFunctionCore_Hip(op, false, &assembled, &rstr, request);
579 }
580 
581 //------------------------------------------------------------------------------
582 // Assemble Diagonal Setup
583 //------------------------------------------------------------------------------
584 static inline int CeedOperatorAssembleDiagonalSetup_Hip(CeedOperator op) {
585   Ceed                ceed;
586   CeedInt             num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0;
587   CeedInt             q_comp, num_nodes, num_qpts;
588   CeedEvalMode       *eval_modes_in = NULL, *eval_modes_out = NULL;
589   CeedBasis           basis_in = NULL, basis_out = NULL;
590   CeedQFunctionField *qf_fields;
591   CeedQFunction       qf;
592   CeedOperatorField  *op_fields;
593   CeedOperator_Hip   *impl;
594 
595   CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
596   CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
597   CeedCallBackend(CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields));
598 
599   // Determine active input basis
600   CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL));
601   CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL));
602   for (CeedInt i = 0; i < num_input_fields; i++) {
603     CeedVector vec;
604 
605     CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec));
606     if (vec == CEED_VECTOR_ACTIVE) {
607       CeedBasis    basis;
608       CeedEvalMode eval_mode;
609 
610       CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
611       CeedCheck(!basis_in || basis_in == basis, ceed, CEED_ERROR_BACKEND,
612                 "Backend does not implement operator diagonal assembly with multiple active bases");
613       basis_in = basis;
614       CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
615       CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp));
616       if (eval_mode != CEED_EVAL_WEIGHT) {
617         // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF assembly
618         CeedCallBackend(CeedRealloc(num_eval_modes_in + q_comp, &eval_modes_in));
619         for (CeedInt d = 0; d < q_comp; d++) eval_modes_in[num_eval_modes_in + d] = eval_mode;
620         num_eval_modes_in += q_comp;
621       }
622     }
623   }
624 
625   // Determine active output basis
626   CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields));
627   CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields));
628   for (CeedInt i = 0; i < num_output_fields; i++) {
629     CeedVector vec;
630 
631     CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec));
632     if (vec == CEED_VECTOR_ACTIVE) {
633       CeedBasis    basis;
634       CeedEvalMode eval_mode;
635 
636       CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
637       CeedCheck(!basis_out || basis_out == basis, ceed, CEED_ERROR_BACKEND,
638                 "Backend does not implement operator diagonal assembly with multiple active bases");
639       basis_out = basis;
640       CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
641       CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp));
642       if (eval_mode != CEED_EVAL_WEIGHT) {
643         // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF assembly
644         CeedCallBackend(CeedRealloc(num_eval_modes_out + q_comp, &eval_modes_out));
645         for (CeedInt d = 0; d < q_comp; d++) eval_modes_out[num_eval_modes_out + d] = eval_mode;
646         num_eval_modes_out += q_comp;
647       }
648     }
649   }
650 
651   // Operator data struct
652   CeedCallBackend(CeedOperatorGetData(op, &impl));
653   CeedCallBackend(CeedCalloc(1, &impl->diag));
654   CeedOperatorDiag_Hip *diag = impl->diag;
655 
656   // Basis matrices
657   CeedCallBackend(CeedBasisGetNumNodes(basis_in, &num_nodes));
658   if (basis_in == CEED_BASIS_NONE) num_qpts = num_nodes;
659   else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts));
660   const CeedInt interp_bytes     = num_nodes * num_qpts * sizeof(CeedScalar);
661   const CeedInt eval_modes_bytes = sizeof(CeedEvalMode);
662   bool          has_eval_none    = false;
663 
664   // CEED_EVAL_NONE
665   for (CeedInt i = 0; i < num_eval_modes_in; i++) has_eval_none = has_eval_none || (eval_modes_in[i] == CEED_EVAL_NONE);
666   for (CeedInt i = 0; i < num_eval_modes_out; i++) has_eval_none = has_eval_none || (eval_modes_out[i] == CEED_EVAL_NONE);
667   if (has_eval_none) {
668     CeedScalar *identity = NULL;
669 
670     CeedCallBackend(CeedCalloc(num_nodes * num_qpts, &identity));
671     for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) identity[i * num_nodes + i] = 1.0;
672     CeedCallHip(ceed, hipMalloc((void **)&diag->d_identity, interp_bytes));
673     CeedCallHip(ceed, hipMemcpy(diag->d_identity, identity, interp_bytes, hipMemcpyHostToDevice));
674     CeedCallBackend(CeedFree(&identity));
675   }
676 
677   // CEED_EVAL_INTERP, CEED_EVAL_GRAD, CEED_EVAL_DIV, and CEED_EVAL_CURL
678   for (CeedInt in = 0; in < 2; in++) {
679     CeedFESpace fespace;
680     CeedBasis   basis = in ? basis_in : basis_out;
681 
682     CeedCallBackend(CeedBasisGetFESpace(basis, &fespace));
683     switch (fespace) {
684       case CEED_FE_SPACE_H1: {
685         CeedInt           q_comp_interp, q_comp_grad;
686         const CeedScalar *interp, *grad;
687         CeedScalar       *d_interp, *d_grad;
688 
689         CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp));
690         CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp_grad));
691 
692         CeedCallBackend(CeedBasisGetInterp(basis, &interp));
693         CeedCallHip(ceed, hipMalloc((void **)&d_interp, interp_bytes * q_comp_interp));
694         CeedCallHip(ceed, hipMemcpy(d_interp, interp, interp_bytes * q_comp_interp, hipMemcpyHostToDevice));
695         CeedCallBackend(CeedBasisGetGrad(basis, &grad));
696         CeedCallHip(ceed, hipMalloc((void **)&d_grad, interp_bytes * q_comp_grad));
697         CeedCallHip(ceed, hipMemcpy(d_grad, grad, interp_bytes * q_comp_grad, hipMemcpyHostToDevice));
698         if (in) {
699           diag->d_interp_in = d_interp;
700           diag->d_grad_in   = d_grad;
701         } else {
702           diag->d_interp_out = d_interp;
703           diag->d_grad_out   = d_grad;
704         }
705       } break;
706       case CEED_FE_SPACE_HDIV: {
707         CeedInt           q_comp_interp, q_comp_div;
708         const CeedScalar *interp, *div;
709         CeedScalar       *d_interp, *d_div;
710 
711         CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp));
712         CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp_div));
713 
714         CeedCallBackend(CeedBasisGetInterp(basis, &interp));
715         CeedCallHip(ceed, hipMalloc((void **)&d_interp, interp_bytes * q_comp_interp));
716         CeedCallHip(ceed, hipMemcpy(d_interp, interp, interp_bytes * q_comp_interp, hipMemcpyHostToDevice));
717         CeedCallBackend(CeedBasisGetDiv(basis, &div));
718         CeedCallHip(ceed, hipMalloc((void **)&d_div, interp_bytes * q_comp_div));
719         CeedCallHip(ceed, hipMemcpy(d_div, div, interp_bytes * q_comp_div, hipMemcpyHostToDevice));
720         if (in) {
721           diag->d_interp_in = d_interp;
722           diag->d_div_in    = d_div;
723         } else {
724           diag->d_interp_out = d_interp;
725           diag->d_div_out    = d_div;
726         }
727       } break;
728       case CEED_FE_SPACE_HCURL: {
729         CeedInt           q_comp_interp, q_comp_curl;
730         const CeedScalar *interp, *curl;
731         CeedScalar       *d_interp, *d_curl;
732 
733         CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp));
734         CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp_curl));
735 
736         CeedCallBackend(CeedBasisGetInterp(basis, &interp));
737         CeedCallHip(ceed, hipMalloc((void **)&d_interp, interp_bytes * q_comp_interp));
738         CeedCallHip(ceed, hipMemcpy(d_interp, interp, interp_bytes * q_comp_interp, hipMemcpyHostToDevice));
739         CeedCallBackend(CeedBasisGetCurl(basis, &curl));
740         CeedCallHip(ceed, hipMalloc((void **)&d_curl, interp_bytes * q_comp_curl));
741         CeedCallHip(ceed, hipMemcpy(d_curl, curl, interp_bytes * q_comp_curl, hipMemcpyHostToDevice));
742         if (in) {
743           diag->d_interp_in = d_interp;
744           diag->d_curl_in   = d_curl;
745         } else {
746           diag->d_interp_out = d_interp;
747           diag->d_curl_out   = d_curl;
748         }
749       } break;
750     }
751   }
752 
753   // Arrays of eval_modes
754   CeedCallHip(ceed, hipMalloc((void **)&diag->d_eval_modes_in, num_eval_modes_in * eval_modes_bytes));
755   CeedCallHip(ceed, hipMemcpy(diag->d_eval_modes_in, eval_modes_in, num_eval_modes_in * eval_modes_bytes, hipMemcpyHostToDevice));
756   CeedCallHip(ceed, hipMalloc((void **)&diag->d_eval_modes_out, num_eval_modes_out * eval_modes_bytes));
757   CeedCallHip(ceed, hipMemcpy(diag->d_eval_modes_out, eval_modes_out, num_eval_modes_out * eval_modes_bytes, hipMemcpyHostToDevice));
758   CeedCallBackend(CeedFree(&eval_modes_in));
759   CeedCallBackend(CeedFree(&eval_modes_out));
760   return CEED_ERROR_SUCCESS;
761 }
762 
763 //------------------------------------------------------------------------------
764 // Assemble Diagonal Setup (Compilation)
765 //------------------------------------------------------------------------------
766 static inline int CeedOperatorAssembleDiagonalSetupCompile_Hip(CeedOperator op, CeedInt use_ceedsize_idx, const bool is_point_block) {
767   Ceed                ceed;
768   char               *diagonal_kernel_source;
769   const char         *diagonal_kernel_path;
770   CeedInt             num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0;
771   CeedInt             num_comp, q_comp, num_nodes, num_qpts;
772   CeedBasis           basis_in = NULL, basis_out = NULL;
773   CeedQFunctionField *qf_fields;
774   CeedQFunction       qf;
775   CeedOperatorField  *op_fields;
776   CeedOperator_Hip   *impl;
777 
778   CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
779   CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
780   CeedCallBackend(CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields));
781 
782   // Determine active input basis
783   CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL));
784   CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL));
785   for (CeedInt i = 0; i < num_input_fields; i++) {
786     CeedVector vec;
787 
788     CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec));
789     if (vec == CEED_VECTOR_ACTIVE) {
790       CeedEvalMode eval_mode;
791 
792       CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis_in));
793       CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
794       CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp));
795       if (eval_mode != CEED_EVAL_WEIGHT) {
796         num_eval_modes_in += q_comp;
797       }
798     }
799   }
800 
801   // Determine active output basis
802   CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields));
803   CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields));
804   for (CeedInt i = 0; i < num_output_fields; i++) {
805     CeedVector vec;
806 
807     CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec));
808     if (vec == CEED_VECTOR_ACTIVE) {
809       CeedEvalMode eval_mode;
810 
811       CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis_out));
812       CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
813       CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp));
814       if (eval_mode != CEED_EVAL_WEIGHT) {
815         num_eval_modes_out += q_comp;
816       }
817     }
818   }
819 
820   // Operator data struct
821   CeedCallBackend(CeedOperatorGetData(op, &impl));
822   CeedOperatorDiag_Hip *diag = impl->diag;
823 
824   // Assemble kernel
825   hipModule_t *module          = is_point_block ? &diag->module_point_block : &diag->module;
826   CeedInt      elems_per_block = 1;
827   CeedCallBackend(CeedBasisGetNumNodes(basis_in, &num_nodes));
828   CeedCallBackend(CeedBasisGetNumComponents(basis_in, &num_comp));
829   if (basis_in == CEED_BASIS_NONE) num_qpts = num_nodes;
830   else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts));
831   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/hip/hip-ref-operator-assemble-diagonal.h", &diagonal_kernel_path));
832   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Diagonal Assembly Kernel Source -----\n");
833   CeedCallBackend(CeedLoadSourceToBuffer(ceed, diagonal_kernel_path, &diagonal_kernel_source));
834   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Diagonal Assembly Source Complete! -----\n");
835   CeedCallHip(ceed, CeedCompile_Hip(ceed, diagonal_kernel_source, module, 8, "NUM_EVAL_MODES_IN", num_eval_modes_in, "NUM_EVAL_MODES_OUT",
836                                     num_eval_modes_out, "NUM_COMP", num_comp, "NUM_NODES", num_nodes, "NUM_QPTS", num_qpts, "USE_CEEDSIZE",
837                                     use_ceedsize_idx, "USE_POINT_BLOCK", is_point_block ? 1 : 0, "BLOCK_SIZE", num_nodes * elems_per_block));
838   CeedCallHip(ceed, CeedGetKernel_Hip(ceed, *module, "LinearDiagonal", is_point_block ? &diag->LinearPointBlock : &diag->LinearDiagonal));
839   CeedCallBackend(CeedFree(&diagonal_kernel_path));
840   CeedCallBackend(CeedFree(&diagonal_kernel_source));
841   return CEED_ERROR_SUCCESS;
842 }
843 
844 //------------------------------------------------------------------------------
845 // Assemble Diagonal Core
846 //------------------------------------------------------------------------------
847 static inline int CeedOperatorAssembleDiagonalCore_Hip(CeedOperator op, CeedVector assembled, CeedRequest *request, const bool is_point_block) {
848   Ceed                ceed;
849   CeedInt             num_elem, num_nodes;
850   CeedScalar         *elem_diag_array;
851   const CeedScalar   *assembled_qf_array;
852   CeedVector          assembled_qf   = NULL, elem_diag;
853   CeedElemRestriction assembled_rstr = NULL, rstr_in, rstr_out, diag_rstr;
854   CeedOperator_Hip   *impl;
855 
856   CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
857   CeedCallBackend(CeedOperatorGetData(op, &impl));
858 
859   // Assemble QFunction
860   CeedCallBackend(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_rstr, request));
861   CeedCallBackend(CeedElemRestrictionDestroy(&assembled_rstr));
862   CeedCallBackend(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &assembled_qf_array));
863 
864   // Setup
865   if (!impl->diag) CeedCallBackend(CeedOperatorAssembleDiagonalSetup_Hip(op));
866   CeedOperatorDiag_Hip *diag = impl->diag;
867 
868   assert(diag != NULL);
869 
870   // Assemble kernel if needed
871   if ((!is_point_block && !diag->LinearDiagonal) || (is_point_block && !diag->LinearPointBlock)) {
872     CeedSize assembled_length, assembled_qf_length;
873     CeedInt  use_ceedsize_idx = 0;
874     CeedCallBackend(CeedVectorGetLength(assembled, &assembled_length));
875     CeedCallBackend(CeedVectorGetLength(assembled_qf, &assembled_qf_length));
876     if ((assembled_length > INT_MAX) || (assembled_qf_length > INT_MAX)) use_ceedsize_idx = 1;
877 
878     CeedCallBackend(CeedOperatorAssembleDiagonalSetupCompile_Hip(op, use_ceedsize_idx, is_point_block));
879   }
880 
881   // Restriction and diagonal vector
882   CeedCallBackend(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out));
883   CeedCheck(rstr_in == rstr_out, ceed, CEED_ERROR_BACKEND,
884             "Cannot assemble operator diagonal with different input and output active element restrictions");
885   if (!is_point_block && !diag->diag_rstr) {
886     CeedCallBackend(CeedElemRestrictionCreateUnsignedCopy(rstr_out, &diag->diag_rstr));
887     CeedCallBackend(CeedElemRestrictionCreateVector(diag->diag_rstr, NULL, &diag->elem_diag));
888   } else if (is_point_block && !diag->point_block_diag_rstr) {
889     CeedCallBackend(CeedOperatorCreateActivePointBlockRestriction(rstr_out, &diag->point_block_diag_rstr));
890     CeedCallBackend(CeedElemRestrictionCreateVector(diag->point_block_diag_rstr, NULL, &diag->point_block_elem_diag));
891   }
892   diag_rstr = is_point_block ? diag->point_block_diag_rstr : diag->diag_rstr;
893   elem_diag = is_point_block ? diag->point_block_elem_diag : diag->elem_diag;
894   CeedCallBackend(CeedVectorSetValue(elem_diag, 0.0));
895 
896   // Only assemble diagonal if the basis has nodes, otherwise inputs are null pointers
897   CeedCallBackend(CeedElemRestrictionGetElementSize(diag_rstr, &num_nodes));
898   if (num_nodes > 0) {
899     // Assemble element operator diagonals
900     CeedCallBackend(CeedVectorGetArray(elem_diag, CEED_MEM_DEVICE, &elem_diag_array));
901     CeedCallBackend(CeedElemRestrictionGetNumElements(diag_rstr, &num_elem));
902 
903     // Compute the diagonal of B^T D B
904     CeedInt elems_per_block = 1;
905     CeedInt grid            = CeedDivUpInt(num_elem, elems_per_block);
906     void   *args[]          = {(void *)&num_elem,      &diag->d_identity,       &diag->d_interp_in,  &diag->d_grad_in, &diag->d_div_in,
907                                &diag->d_curl_in,       &diag->d_interp_out,     &diag->d_grad_out,   &diag->d_div_out, &diag->d_curl_out,
908                                &diag->d_eval_modes_in, &diag->d_eval_modes_out, &assembled_qf_array, &elem_diag_array};
909 
910     if (is_point_block) {
911       CeedCallBackend(CeedRunKernelDim_Hip(ceed, diag->LinearPointBlock, grid, num_nodes, 1, elems_per_block, args));
912     } else {
913       CeedCallBackend(CeedRunKernelDim_Hip(ceed, diag->LinearDiagonal, grid, num_nodes, 1, elems_per_block, args));
914     }
915 
916     // Restore arrays
917     CeedCallBackend(CeedVectorRestoreArray(elem_diag, &elem_diag_array));
918     CeedCallBackend(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array));
919   }
920 
921   // Assemble local operator diagonal
922   CeedCallBackend(CeedElemRestrictionApply(diag_rstr, CEED_TRANSPOSE, elem_diag, assembled, request));
923 
924   // Cleanup
925   CeedCallBackend(CeedVectorDestroy(&assembled_qf));
926   return CEED_ERROR_SUCCESS;
927 }
928 
929 //------------------------------------------------------------------------------
930 // Assemble Linear Diagonal
931 //------------------------------------------------------------------------------
932 static int CeedOperatorLinearAssembleAddDiagonal_Hip(CeedOperator op, CeedVector assembled, CeedRequest *request) {
933   CeedCallBackend(CeedOperatorAssembleDiagonalCore_Hip(op, assembled, request, false));
934   return CEED_ERROR_SUCCESS;
935 }
936 
937 //------------------------------------------------------------------------------
938 // Assemble Linear Point Block Diagonal
939 //------------------------------------------------------------------------------
940 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip(CeedOperator op, CeedVector assembled, CeedRequest *request) {
941   CeedCallBackend(CeedOperatorAssembleDiagonalCore_Hip(op, assembled, request, true));
942   return CEED_ERROR_SUCCESS;
943 }
944 
945 //------------------------------------------------------------------------------
946 // Single Operator Assembly Setup
947 //------------------------------------------------------------------------------
948 static int CeedSingleOperatorAssembleSetup_Hip(CeedOperator op, CeedInt use_ceedsize_idx) {
949   Ceed                ceed;
950   char               *assembly_kernel_source;
951   const char         *assembly_kernel_path;
952   CeedInt             num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0;
953   CeedInt             elem_size_in, num_qpts_in = 0, num_comp_in, elem_size_out, num_qpts_out, num_comp_out, q_comp;
954   CeedEvalMode       *eval_modes_in = NULL, *eval_modes_out = NULL;
955   CeedElemRestriction rstr_in = NULL, rstr_out = NULL;
956   CeedBasis           basis_in = NULL, basis_out = NULL;
957   CeedQFunctionField *qf_fields;
958   CeedQFunction       qf;
959   CeedOperatorField  *input_fields, *output_fields;
960   CeedOperator_Hip   *impl;
961 
962   CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
963   CeedCallBackend(CeedOperatorGetData(op, &impl));
964 
965   // Get intput and output fields
966   CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields));
967 
968   // Determine active input basis eval mode
969   CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
970   CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL));
971   for (CeedInt i = 0; i < num_input_fields; i++) {
972     CeedVector vec;
973 
974     CeedCallBackend(CeedOperatorFieldGetVector(input_fields[i], &vec));
975     if (vec == CEED_VECTOR_ACTIVE) {
976       CeedBasis    basis;
977       CeedEvalMode eval_mode;
978 
979       CeedCallBackend(CeedOperatorFieldGetBasis(input_fields[i], &basis));
980       CeedCheck(!basis_in || basis_in == basis, ceed, CEED_ERROR_BACKEND, "Backend does not implement operator assembly with multiple active bases");
981       basis_in = basis;
982       CeedCallBackend(CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_in));
983       CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in));
984       if (basis_in == CEED_BASIS_NONE) num_qpts_in = elem_size_in;
985       else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts_in));
986       CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
987       CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp));
988       if (eval_mode != CEED_EVAL_WEIGHT) {
989         // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly
990         CeedCallBackend(CeedRealloc(num_eval_modes_in + q_comp, &eval_modes_in));
991         for (CeedInt d = 0; d < q_comp; d++) {
992           eval_modes_in[num_eval_modes_in + d] = eval_mode;
993         }
994         num_eval_modes_in += q_comp;
995       }
996     }
997   }
998 
999   // Determine active output basis; basis_out and rstr_out only used if same as input, TODO
1000   CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields));
1001   for (CeedInt i = 0; i < num_output_fields; i++) {
1002     CeedVector vec;
1003 
1004     CeedCallBackend(CeedOperatorFieldGetVector(output_fields[i], &vec));
1005     if (vec == CEED_VECTOR_ACTIVE) {
1006       CeedBasis    basis;
1007       CeedEvalMode eval_mode;
1008 
1009       CeedCallBackend(CeedOperatorFieldGetBasis(output_fields[i], &basis));
1010       CeedCheck(!basis_out || basis_out == basis, ceed, CEED_ERROR_BACKEND,
1011                 "Backend does not implement operator assembly with multiple active bases");
1012       basis_out = basis;
1013       CeedCallBackend(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_out));
1014       CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out));
1015       if (basis_out == CEED_BASIS_NONE) num_qpts_out = elem_size_out;
1016       else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_out, &num_qpts_out));
1017       CeedCheck(num_qpts_in == num_qpts_out, ceed, CEED_ERROR_UNSUPPORTED,
1018                 "Active input and output bases must have the same number of quadrature points");
1019       CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
1020       CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp));
1021       if (eval_mode != CEED_EVAL_WEIGHT) {
1022         // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly
1023         CeedCallBackend(CeedRealloc(num_eval_modes_out + q_comp, &eval_modes_out));
1024         for (CeedInt d = 0; d < q_comp; d++) {
1025           eval_modes_out[num_eval_modes_out + d] = eval_mode;
1026         }
1027         num_eval_modes_out += q_comp;
1028       }
1029     }
1030   }
1031   CeedCheck(num_eval_modes_in > 0 && num_eval_modes_out > 0, ceed, CEED_ERROR_UNSUPPORTED, "Cannot assemble operator without inputs/outputs");
1032 
1033   CeedCallBackend(CeedCalloc(1, &impl->asmb));
1034   CeedOperatorAssemble_Hip *asmb = impl->asmb;
1035   asmb->elems_per_block          = 1;
1036   asmb->block_size_x             = elem_size_in;
1037   asmb->block_size_y             = elem_size_out;
1038 
1039   bool fallback = asmb->block_size_x * asmb->block_size_y * asmb->elems_per_block > 1024;
1040 
1041   if (fallback) {
1042     // Use fallback kernel with 1D threadblock
1043     asmb->block_size_y = 1;
1044   }
1045 
1046   // Compile kernels
1047   CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_in, &num_comp_in));
1048   CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_out, &num_comp_out));
1049   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/hip/hip-ref-operator-assemble.h", &assembly_kernel_path));
1050   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Assembly Kernel Source -----\n");
1051   CeedCallBackend(CeedLoadSourceToBuffer(ceed, assembly_kernel_path, &assembly_kernel_source));
1052   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Assembly Source Complete! -----\n");
1053   CeedCallBackend(CeedCompile_Hip(ceed, assembly_kernel_source, &asmb->module, 10, "NUM_EVAL_MODES_IN", num_eval_modes_in, "NUM_EVAL_MODES_OUT",
1054                                   num_eval_modes_out, "NUM_COMP_IN", num_comp_in, "NUM_COMP_OUT", num_comp_out, "NUM_NODES_IN", elem_size_in,
1055                                   "NUM_NODES_OUT", elem_size_out, "NUM_QPTS", num_qpts_in, "BLOCK_SIZE",
1056                                   asmb->block_size_x * asmb->block_size_y * asmb->elems_per_block, "BLOCK_SIZE_Y", asmb->block_size_y, "USE_CEEDSIZE",
1057                                   use_ceedsize_idx));
1058   CeedCallBackend(CeedGetKernel_Hip(ceed, asmb->module, "LinearAssemble", &asmb->LinearAssemble));
1059   CeedCallBackend(CeedFree(&assembly_kernel_path));
1060   CeedCallBackend(CeedFree(&assembly_kernel_source));
1061 
1062   // Load into B_in, in order that they will be used in eval_modes_in
1063   {
1064     const CeedInt in_bytes           = elem_size_in * num_qpts_in * num_eval_modes_in * sizeof(CeedScalar);
1065     CeedInt       d_in               = 0;
1066     CeedEvalMode  eval_modes_in_prev = CEED_EVAL_NONE;
1067     bool          has_eval_none      = false;
1068     CeedScalar   *identity           = NULL;
1069 
1070     for (CeedInt i = 0; i < num_eval_modes_in; i++) {
1071       has_eval_none = has_eval_none || (eval_modes_in[i] == CEED_EVAL_NONE);
1072     }
1073     if (has_eval_none) {
1074       CeedCallBackend(CeedCalloc(elem_size_in * num_qpts_in, &identity));
1075       for (CeedInt i = 0; i < (elem_size_in < num_qpts_in ? elem_size_in : num_qpts_in); i++) identity[i * elem_size_in + i] = 1.0;
1076     }
1077 
1078     CeedCallHip(ceed, hipMalloc((void **)&asmb->d_B_in, in_bytes));
1079     for (CeedInt i = 0; i < num_eval_modes_in; i++) {
1080       const CeedScalar *h_B_in;
1081 
1082       CeedCallBackend(CeedOperatorGetBasisPointer(basis_in, eval_modes_in[i], identity, &h_B_in));
1083       CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_modes_in[i], &q_comp));
1084       if (q_comp > 1) {
1085         if (i == 0 || eval_modes_in[i] != eval_modes_in_prev) d_in = 0;
1086         else h_B_in = &h_B_in[(++d_in) * elem_size_in * num_qpts_in];
1087       }
1088       eval_modes_in_prev = eval_modes_in[i];
1089 
1090       CeedCallHip(ceed, hipMemcpy(&asmb->d_B_in[i * elem_size_in * num_qpts_in], h_B_in, elem_size_in * num_qpts_in * sizeof(CeedScalar),
1091                                   hipMemcpyHostToDevice));
1092     }
1093 
1094     if (identity) {
1095       CeedCallBackend(CeedFree(&identity));
1096     }
1097   }
1098 
1099   // Load into B_out, in order that they will be used in eval_modes_out
1100   {
1101     const CeedInt out_bytes           = elem_size_out * num_qpts_out * num_eval_modes_out * sizeof(CeedScalar);
1102     CeedInt       d_out               = 0;
1103     CeedEvalMode  eval_modes_out_prev = CEED_EVAL_NONE;
1104     bool          has_eval_none       = false;
1105     CeedScalar   *identity            = NULL;
1106 
1107     for (CeedInt i = 0; i < num_eval_modes_out; i++) {
1108       has_eval_none = has_eval_none || (eval_modes_out[i] == CEED_EVAL_NONE);
1109     }
1110     if (has_eval_none) {
1111       CeedCallBackend(CeedCalloc(elem_size_out * num_qpts_out, &identity));
1112       for (CeedInt i = 0; i < (elem_size_out < num_qpts_out ? elem_size_out : num_qpts_out); i++) identity[i * elem_size_out + i] = 1.0;
1113     }
1114 
1115     CeedCallHip(ceed, hipMalloc((void **)&asmb->d_B_out, out_bytes));
1116     for (CeedInt i = 0; i < num_eval_modes_out; i++) {
1117       const CeedScalar *h_B_out;
1118 
1119       CeedCallBackend(CeedOperatorGetBasisPointer(basis_out, eval_modes_out[i], identity, &h_B_out));
1120       CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_modes_out[i], &q_comp));
1121       if (q_comp > 1) {
1122         if (i == 0 || eval_modes_out[i] != eval_modes_out_prev) d_out = 0;
1123         else h_B_out = &h_B_out[(++d_out) * elem_size_out * num_qpts_out];
1124       }
1125       eval_modes_out_prev = eval_modes_out[i];
1126 
1127       CeedCallHip(ceed, hipMemcpy(&asmb->d_B_out[i * elem_size_out * num_qpts_out], h_B_out, elem_size_out * num_qpts_out * sizeof(CeedScalar),
1128                                   hipMemcpyHostToDevice));
1129     }
1130 
1131     if (identity) {
1132       CeedCallBackend(CeedFree(&identity));
1133     }
1134   }
1135   return CEED_ERROR_SUCCESS;
1136 }
1137 
1138 //------------------------------------------------------------------------------
1139 // Assemble matrix data for COO matrix of assembled operator.
1140 // The sparsity pattern is set by CeedOperatorLinearAssembleSymbolic.
1141 //
1142 // Note that this (and other assembly routines) currently assume only one active input restriction/basis per operator (could have multiple basis eval
1143 // modes).
1144 // TODO: allow multiple active input restrictions/basis objects
1145 //------------------------------------------------------------------------------
1146 static int CeedSingleOperatorAssemble_Hip(CeedOperator op, CeedInt offset, CeedVector values) {
1147   Ceed                ceed;
1148   CeedSize            values_length = 0, assembled_qf_length = 0;
1149   CeedInt             use_ceedsize_idx = 0, num_elem_in, num_elem_out, elem_size_in, elem_size_out;
1150   CeedScalar         *values_array;
1151   const CeedScalar   *assembled_qf_array;
1152   CeedVector          assembled_qf   = NULL;
1153   CeedElemRestriction assembled_rstr = NULL, rstr_in, rstr_out;
1154   CeedRestrictionType rstr_type_in, rstr_type_out;
1155   const bool         *orients_in = NULL, *orients_out = NULL;
1156   const CeedInt8     *curl_orients_in = NULL, *curl_orients_out = NULL;
1157   CeedOperator_Hip   *impl;
1158 
1159   CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
1160   CeedCallBackend(CeedOperatorGetData(op, &impl));
1161 
1162   // Assemble QFunction
1163   CeedCallBackend(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_rstr, CEED_REQUEST_IMMEDIATE));
1164   CeedCallBackend(CeedElemRestrictionDestroy(&assembled_rstr));
1165   CeedCallBackend(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &assembled_qf_array));
1166 
1167   CeedCallBackend(CeedVectorGetLength(values, &values_length));
1168   CeedCallBackend(CeedVectorGetLength(assembled_qf, &assembled_qf_length));
1169   if ((values_length > INT_MAX) || (assembled_qf_length > INT_MAX)) use_ceedsize_idx = 1;
1170 
1171   // Setup
1172   if (!impl->asmb) CeedCallBackend(CeedSingleOperatorAssembleSetup_Hip(op, use_ceedsize_idx));
1173   CeedOperatorAssemble_Hip *asmb = impl->asmb;
1174 
1175   assert(asmb != NULL);
1176 
1177   // Assemble element operator
1178   CeedCallBackend(CeedVectorGetArray(values, CEED_MEM_DEVICE, &values_array));
1179   values_array += offset;
1180 
1181   CeedCallBackend(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out));
1182   CeedCallBackend(CeedElemRestrictionGetNumElements(rstr_in, &num_elem_in));
1183   CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in));
1184 
1185   CeedCallBackend(CeedElemRestrictionGetType(rstr_in, &rstr_type_in));
1186   if (rstr_type_in == CEED_RESTRICTION_ORIENTED) {
1187     CeedCallBackend(CeedElemRestrictionGetOrientations(rstr_in, CEED_MEM_DEVICE, &orients_in));
1188   } else if (rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) {
1189     CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr_in, CEED_MEM_DEVICE, &curl_orients_in));
1190   }
1191 
1192   if (rstr_in != rstr_out) {
1193     CeedCallBackend(CeedElemRestrictionGetNumElements(rstr_out, &num_elem_out));
1194     CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED,
1195               "Active input and output operator restrictions must have the same number of elements");
1196     CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out));
1197 
1198     CeedCallBackend(CeedElemRestrictionGetType(rstr_out, &rstr_type_out));
1199     if (rstr_type_out == CEED_RESTRICTION_ORIENTED) {
1200       CeedCallBackend(CeedElemRestrictionGetOrientations(rstr_out, CEED_MEM_DEVICE, &orients_out));
1201     } else if (rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) {
1202       CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr_out, CEED_MEM_DEVICE, &curl_orients_out));
1203     }
1204   } else {
1205     elem_size_out    = elem_size_in;
1206     orients_out      = orients_in;
1207     curl_orients_out = curl_orients_in;
1208   }
1209 
1210   // Compute B^T D B
1211   CeedInt shared_mem =
1212       ((curl_orients_in || curl_orients_out ? elem_size_in * elem_size_out : 0) + (curl_orients_in ? elem_size_in * asmb->block_size_y : 0)) *
1213       sizeof(CeedScalar);
1214   CeedInt grid   = CeedDivUpInt(num_elem_in, asmb->elems_per_block);
1215   void   *args[] = {(void *)&num_elem_in, &asmb->d_B_in,     &asmb->d_B_out,      &orients_in,  &curl_orients_in,
1216                     &orients_out,         &curl_orients_out, &assembled_qf_array, &values_array};
1217 
1218   CeedCallBackend(
1219       CeedRunKernelDimShared_Hip(ceed, asmb->LinearAssemble, grid, asmb->block_size_x, asmb->block_size_y, asmb->elems_per_block, shared_mem, args));
1220 
1221   // Restore arrays
1222   CeedCallBackend(CeedVectorRestoreArray(values, &values_array));
1223   CeedCallBackend(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array));
1224 
1225   // Cleanup
1226   CeedCallBackend(CeedVectorDestroy(&assembled_qf));
1227   if (rstr_type_in == CEED_RESTRICTION_ORIENTED) {
1228     CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr_in, &orients_in));
1229   } else if (rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) {
1230     CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr_in, &curl_orients_in));
1231   }
1232   if (rstr_in != rstr_out) {
1233     if (rstr_type_out == CEED_RESTRICTION_ORIENTED) {
1234       CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr_out, &orients_out));
1235     } else if (rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) {
1236       CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr_out, &curl_orients_out));
1237     }
1238   }
1239   return CEED_ERROR_SUCCESS;
1240 }
1241 
1242 //------------------------------------------------------------------------------
1243 // Create operator
1244 //------------------------------------------------------------------------------
1245 int CeedOperatorCreate_Hip(CeedOperator op) {
1246   Ceed              ceed;
1247   CeedOperator_Hip *impl;
1248 
1249   CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
1250   CeedCallBackend(CeedCalloc(1, &impl));
1251   CeedCallBackend(CeedOperatorSetData(op, impl));
1252   CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Hip));
1253   CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Hip));
1254   CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonal_Hip));
1255   CeedCallBackend(
1256       CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddPointBlockDiagonal", CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip));
1257   CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleSingle", CeedSingleOperatorAssemble_Hip));
1258   CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Hip));
1259   CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Hip));
1260   return CEED_ERROR_SUCCESS;
1261 }
1262 
1263 //------------------------------------------------------------------------------
1264