xref: /libCEED/backends/ref/ceed-ref-operator.c (revision 74653b0dfd9ea5453a7e0d931542e0b08b318f8b) 
1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
3 // All Rights reserved. See files LICENSE and NOTICE for details.
4 //
5 // This file is part of CEED, a collection of benchmarks, miniapps, software
6 // libraries and APIs for efficient high-order finite element and spectral
7 // element discretizations for exascale applications. For more information and
8 // source code availability see http://github.com/ceed.
9 //
10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11 // a collaborative effort of two U.S. Department of Energy organizations (Office
12 // of Science and the National Nuclear Security Administration) responsible for
13 // the planning and preparation of a capable exascale ecosystem, including
14 // software, applications, hardware, advanced system engineering and early
15 // testbed platforms, in support of the nation's exascale computing imperative.
16 
17 #include <ceed/ceed.h>
18 #include <ceed/backend.h>
19 #include <math.h>
20 #include <stdbool.h>
21 #include <stddef.h>
22 #include <stdint.h>
23 #include "ceed-ref.h"
24 
25 //------------------------------------------------------------------------------
26 // Setup Input/Output Fields
27 //------------------------------------------------------------------------------
28 static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op,
29                                        bool inOrOut,
30                                        CeedVector *full_evecs, CeedVector *e_vecs,
31                                        CeedVector *q_vecs, CeedInt starte,
32                                        CeedInt num_fields, CeedInt Q) {
33   CeedInt dim, ierr, size, P;
34   Ceed ceed;
35   ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr);
36   CeedBasis basis;
37   CeedElemRestriction elem_restr;
38   CeedOperatorField *op_fields;
39   CeedQFunctionField *qf_fields;
40   if (inOrOut) {
41     ierr = CeedOperatorGetFields(op, NULL, &op_fields); CeedChkBackend(ierr);
42     ierr = CeedQFunctionGetFields(qf, NULL, &qf_fields); CeedChkBackend(ierr);
43   } else {
44     ierr = CeedOperatorGetFields(op, &op_fields, NULL); CeedChkBackend(ierr);
45     ierr = CeedQFunctionGetFields(qf, &qf_fields, NULL); CeedChkBackend(ierr);
46   }
47 
48   // Loop over fields
49   for (CeedInt i=0; i<num_fields; i++) {
50     CeedEvalMode eval_mode;
51     ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode);
52     CeedChkBackend(ierr);
53 
54     if (eval_mode != CEED_EVAL_WEIGHT) {
55       ierr = CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_restr);
56       CeedChkBackend(ierr);
57       ierr = CeedElemRestrictionCreateVector(elem_restr, NULL,
58                                              &full_evecs[i+starte]);
59       CeedChkBackend(ierr);
60     }
61 
62     switch(eval_mode) {
63     case CEED_EVAL_NONE:
64       ierr = CeedQFunctionFieldGetSize(qf_fields[i], &size); CeedChkBackend(ierr);
65       ierr = CeedVectorCreate(ceed, Q*size, &q_vecs[i]); CeedChkBackend(ierr);
66       break;
67     case CEED_EVAL_INTERP:
68       ierr = CeedQFunctionFieldGetSize(qf_fields[i], &size); CeedChkBackend(ierr);
69       ierr = CeedElemRestrictionGetElementSize(elem_restr, &P);
70       CeedChkBackend(ierr);
71       ierr = CeedVectorCreate(ceed, P*size, &e_vecs[i]); CeedChkBackend(ierr);
72       ierr = CeedVectorCreate(ceed, Q*size, &q_vecs[i]); CeedChkBackend(ierr);
73       break;
74     case CEED_EVAL_GRAD:
75       ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis); CeedChkBackend(ierr);
76       ierr = CeedQFunctionFieldGetSize(qf_fields[i], &size); CeedChkBackend(ierr);
77       ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr);
78       ierr = CeedElemRestrictionGetElementSize(elem_restr, &P);
79       CeedChkBackend(ierr);
80       ierr = CeedVectorCreate(ceed, P*size/dim, &e_vecs[i]); CeedChkBackend(ierr);
81       ierr = CeedVectorCreate(ceed, Q*size, &q_vecs[i]); CeedChkBackend(ierr);
82       break;
83     case CEED_EVAL_WEIGHT: // Only on input fields
84       ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis); CeedChkBackend(ierr);
85       ierr = CeedVectorCreate(ceed, Q, &q_vecs[i]); CeedChkBackend(ierr);
86       ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT,
87                             CEED_VECTOR_NONE, q_vecs[i]); CeedChkBackend(ierr);
88       break;
89     case CEED_EVAL_DIV:
90       break; // Not implemented
91     case CEED_EVAL_CURL:
92       break; // Not implemented
93     }
94   }
95   return CEED_ERROR_SUCCESS;
96 }
97 
98 //------------------------------------------------------------------------------
99 // Setup Operator
100 //------------------------------------------------------------------------------/*
101 static int CeedOperatorSetup_Ref(CeedOperator op) {
102   int ierr;
103   bool setup_done;
104   ierr = CeedOperatorIsSetupDone(op, &setup_done); CeedChkBackend(ierr);
105   if (setup_done) return CEED_ERROR_SUCCESS;
106   Ceed ceed;
107   ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr);
108   CeedOperator_Ref *impl;
109   ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr);
110   CeedQFunction qf;
111   ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr);
112   CeedInt Q, num_input_fields, num_output_fields;
113   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr);
114   ierr = CeedQFunctionIsIdentity(qf, &impl->identity_qf); CeedChkBackend(ierr);
115   ierr = CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields);
116   CeedChkBackend(ierr);
117   CeedOperatorField *op_input_fields, *op_output_fields;
118   ierr = CeedOperatorGetFields(op, &op_input_fields, &op_output_fields);
119   CeedChkBackend(ierr);
120   CeedQFunctionField *qf_input_fields, *qf_output_fields;
121   ierr = CeedQFunctionGetFields(qf, &qf_input_fields, &qf_output_fields);
122   CeedChkBackend(ierr);
123 
124   // Allocate
125   ierr = CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs);
126   CeedChkBackend(ierr);
127   ierr = CeedCalloc(num_input_fields + num_output_fields, &impl->e_data);
128   CeedChkBackend(ierr);
129 
130   ierr = CeedCalloc(16, &impl->input_state); CeedChkBackend(ierr);
131   ierr = CeedCalloc(16, &impl->e_vecs_in); CeedChkBackend(ierr);
132   ierr = CeedCalloc(16, &impl->e_vecs_out); CeedChkBackend(ierr);
133   ierr = CeedCalloc(16, &impl->q_vecs_in); CeedChkBackend(ierr);
134   ierr = CeedCalloc(16, &impl->q_vecs_out); CeedChkBackend(ierr);
135 
136   impl->num_e_vecs_in = num_input_fields;
137   impl->num_e_vecs_out = num_output_fields;
138 
139   // Set up infield and outfield e_vecs and q_vecs
140   // Infields
141   ierr = CeedOperatorSetupFields_Ref(qf, op, 0, impl->e_vecs,
142                                      impl->e_vecs_in, impl->q_vecs_in, 0,
143                                      num_input_fields, Q);
144   CeedChkBackend(ierr);
145   // Outfields
146   ierr = CeedOperatorSetupFields_Ref(qf, op, 1, impl->e_vecs,
147                                      impl->e_vecs_out, impl->q_vecs_out,
148                                      num_input_fields, num_output_fields, Q);
149   CeedChkBackend(ierr);
150 
151   // Identity QFunctions
152   if (impl->identity_qf) {
153     CeedEvalMode in_mode, out_mode;
154     CeedQFunctionField *in_fields, *out_fields;
155     ierr = CeedQFunctionGetFields(qf, &in_fields, &out_fields);
156     CeedChkBackend(ierr);
157 
158     for (CeedInt i=0; i<num_input_fields; i++) {
159       ierr = CeedQFunctionFieldGetEvalMode(in_fields[i], &in_mode);
160       CeedChkBackend(ierr);
161       ierr = CeedQFunctionFieldGetEvalMode(out_fields[i], &out_mode);
162       CeedChkBackend(ierr);
163 
164       ierr = CeedVectorDestroy(&impl->q_vecs_out[i]); CeedChkBackend(ierr);
165       impl->q_vecs_out[i] = impl->q_vecs_in[i];
166       ierr = CeedVectorAddReference(impl->q_vecs_in[i]); CeedChkBackend(ierr);
167     }
168   }
169 
170   ierr = CeedOperatorSetSetupDone(op); CeedChkBackend(ierr);
171 
172   return CEED_ERROR_SUCCESS;
173 }
174 
175 //------------------------------------------------------------------------------
176 // Setup Operator Inputs
177 //------------------------------------------------------------------------------
178 static inline int CeedOperatorSetupInputs_Ref(CeedInt num_input_fields,
179     CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
180     CeedVector in_vec, const bool skip_active, CeedOperator_Ref *impl,
181     CeedRequest *request) {
182   CeedInt ierr;
183   CeedEvalMode eval_mode;
184   CeedVector vec;
185   CeedElemRestriction elem_restr;
186   uint64_t state;
187 
188   for (CeedInt i=0; i<num_input_fields; i++) {
189     // Get input vector
190     ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec);
191     CeedChkBackend(ierr);
192     if (vec == CEED_VECTOR_ACTIVE) {
193       if (skip_active)
194         continue;
195       else
196         vec = in_vec;
197     }
198 
199     ierr = CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode);
200     CeedChkBackend(ierr);
201     // Restrict and Evec
202     if (eval_mode == CEED_EVAL_WEIGHT) { // Skip
203     } else {
204       // Restrict
205       ierr = CeedVectorGetState(vec, &state); CeedChkBackend(ierr);
206       // Skip restriction if input is unchanged
207       if (state != impl->input_state[i] || vec == in_vec) {
208         ierr = CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_restr);
209         CeedChkBackend(ierr);
210         ierr = CeedElemRestrictionApply(elem_restr, CEED_NOTRANSPOSE, vec,
211                                         impl->e_vecs[i], request); CeedChkBackend(ierr);
212         impl->input_state[i] = state;
213       }
214       // Get evec
215       ierr = CeedVectorGetArrayRead(impl->e_vecs[i], CEED_MEM_HOST,
216                                     (const CeedScalar **) &impl->e_data[i]);
217       CeedChkBackend(ierr);
218     }
219   }
220   return CEED_ERROR_SUCCESS;
221 }
222 
223 //------------------------------------------------------------------------------
224 // Input Basis Action
225 //------------------------------------------------------------------------------
226 static inline int CeedOperatorInputBasis_Ref(CeedInt e, CeedInt Q,
227     CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
228     CeedInt num_input_fields, const bool skip_active, CeedOperator_Ref *impl) {
229   CeedInt ierr;
230   CeedInt dim, elem_size, size;
231   CeedElemRestriction elem_restr;
232   CeedEvalMode eval_mode;
233   CeedBasis basis;
234 
235   for (CeedInt i=0; i<num_input_fields; i++) {
236     // Skip active input
237     if (skip_active) {
238       CeedVector vec;
239       ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec);
240       CeedChkBackend(ierr);
241       if (vec == CEED_VECTOR_ACTIVE)
242         continue;
243     }
244     // Get elem_size, eval_mode, size
245     ierr = CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_restr);
246     CeedChkBackend(ierr);
247     ierr = CeedElemRestrictionGetElementSize(elem_restr, &elem_size);
248     CeedChkBackend(ierr);
249     ierr = CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode);
250     CeedChkBackend(ierr);
251     ierr = CeedQFunctionFieldGetSize(qf_input_fields[i], &size);
252     CeedChkBackend(ierr);
253     // Basis action
254     switch(eval_mode) {
255     case CEED_EVAL_NONE:
256       ierr = CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST,
257                                 CEED_USE_POINTER, &impl->e_data[i][e*Q*size]);
258       CeedChkBackend(ierr);
259       break;
260     case CEED_EVAL_INTERP:
261       ierr = CeedOperatorFieldGetBasis(op_input_fields[i], &basis);
262       CeedChkBackend(ierr);
263       ierr = CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST,
264                                 CEED_USE_POINTER, &impl->e_data[i][e*elem_size*size]);
265       CeedChkBackend(ierr);
266       ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP,
267                             impl->e_vecs_in[i], impl->q_vecs_in[i]); CeedChkBackend(ierr);
268       break;
269     case CEED_EVAL_GRAD:
270       ierr = CeedOperatorFieldGetBasis(op_input_fields[i], &basis);
271       CeedChkBackend(ierr);
272       ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr);
273       ierr = CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST,
274                                 CEED_USE_POINTER, &impl->e_data[i][e*elem_size*size/dim]);
275       CeedChkBackend(ierr);
276       ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE,
277                             CEED_EVAL_GRAD, impl->e_vecs_in[i],
278                             impl->q_vecs_in[i]); CeedChkBackend(ierr);
279       break;
280     case CEED_EVAL_WEIGHT:
281       break;  // No action
282     // LCOV_EXCL_START
283     case CEED_EVAL_DIV:
284     case CEED_EVAL_CURL: {
285       ierr = CeedOperatorFieldGetBasis(op_input_fields[i], &basis);
286       CeedChkBackend(ierr);
287       Ceed ceed;
288       ierr = CeedBasisGetCeed(basis, &ceed); CeedChkBackend(ierr);
289       return CeedError(ceed, CEED_ERROR_BACKEND,
290                        "Ceed evaluation mode not implemented");
291       // LCOV_EXCL_STOP
292     }
293     }
294   }
295   return CEED_ERROR_SUCCESS;
296 }
297 
298 //------------------------------------------------------------------------------
299 // Output Basis Action
300 //------------------------------------------------------------------------------
301 static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q,
302     CeedQFunctionField *qf_output_fields, CeedOperatorField *op_output_fields,
303     CeedInt num_input_fields, CeedInt num_output_fields, CeedOperator op,
304     CeedOperator_Ref *impl) {
305   CeedInt ierr;
306   CeedInt dim, elem_size, size;
307   CeedElemRestriction elem_restr;
308   CeedEvalMode eval_mode;
309   CeedBasis basis;
310 
311   for (CeedInt i=0; i<num_output_fields; i++) {
312     // Get elem_size, eval_mode, size
313     ierr = CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_restr);
314     CeedChkBackend(ierr);
315     ierr = CeedElemRestrictionGetElementSize(elem_restr, &elem_size);
316     CeedChkBackend(ierr);
317     ierr = CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode);
318     CeedChkBackend(ierr);
319     ierr = CeedQFunctionFieldGetSize(qf_output_fields[i], &size);
320     CeedChkBackend(ierr);
321     // Basis action
322     switch(eval_mode) {
323     case CEED_EVAL_NONE:
324       break; // No action
325     case CEED_EVAL_INTERP:
326       ierr = CeedOperatorFieldGetBasis(op_output_fields[i], &basis);
327       CeedChkBackend(ierr);
328       ierr = CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST,
329                                 CEED_USE_POINTER,
330                                 &impl->e_data[i + num_input_fields][e*elem_size*size]);
331       CeedChkBackend(ierr);
332       ierr = CeedBasisApply(basis, 1, CEED_TRANSPOSE,
333                             CEED_EVAL_INTERP, impl->q_vecs_out[i],
334                             impl->e_vecs_out[i]); CeedChkBackend(ierr);
335       break;
336     case CEED_EVAL_GRAD:
337       ierr = CeedOperatorFieldGetBasis(op_output_fields[i], &basis);
338       CeedChkBackend(ierr);
339       ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr);
340       ierr = CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST,
341                                 CEED_USE_POINTER,
342                                 &impl->e_data[i + num_input_fields][e*elem_size*size/dim]);
343       CeedChkBackend(ierr);
344       ierr = CeedBasisApply(basis, 1, CEED_TRANSPOSE,
345                             CEED_EVAL_GRAD, impl->q_vecs_out[i],
346                             impl->e_vecs_out[i]); CeedChkBackend(ierr);
347       break;
348     // LCOV_EXCL_START
349     case CEED_EVAL_WEIGHT: {
350       Ceed ceed;
351       ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr);
352       return CeedError(ceed, CEED_ERROR_BACKEND,
353                        "CEED_EVAL_WEIGHT cannot be an output "
354                        "evaluation mode");
355     }
356     case CEED_EVAL_DIV:
357     case CEED_EVAL_CURL: {
358       Ceed ceed;
359       ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr);
360       return CeedError(ceed, CEED_ERROR_BACKEND,
361                        "Ceed evaluation mode not implemented");
362       // LCOV_EXCL_STOP
363     }
364     }
365   }
366   return CEED_ERROR_SUCCESS;
367 }
368 
369 //------------------------------------------------------------------------------
370 // Restore Input Vectors
371 //------------------------------------------------------------------------------
372 static inline int CeedOperatorRestoreInputs_Ref(CeedInt num_input_fields,
373     CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
374     const bool skip_active, CeedOperator_Ref *impl) {
375   CeedInt ierr;
376   CeedEvalMode eval_mode;
377 
378   for (CeedInt i=0; i<num_input_fields; i++) {
379     // Skip active inputs
380     if (skip_active) {
381       CeedVector vec;
382       ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec);
383       CeedChkBackend(ierr);
384       if (vec == CEED_VECTOR_ACTIVE)
385         continue;
386     }
387     // Restore input
388     ierr = CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode);
389     CeedChkBackend(ierr);
390     if (eval_mode == CEED_EVAL_WEIGHT) { // Skip
391     } else {
392       ierr = CeedVectorRestoreArrayRead(impl->e_vecs[i],
393                                         (const CeedScalar **) &impl->e_data[i]);
394       CeedChkBackend(ierr);
395     }
396   }
397   return CEED_ERROR_SUCCESS;
398 }
399 
400 //------------------------------------------------------------------------------
401 // Operator Apply
402 //------------------------------------------------------------------------------
403 static int CeedOperatorApplyAdd_Ref(CeedOperator op, CeedVector in_vec,
404                                     CeedVector out_vec, CeedRequest *request) {
405   int ierr;
406   CeedOperator_Ref *impl;
407   ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr);
408   CeedQFunction qf;
409   ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr);
410   CeedInt Q, num_elem, num_input_fields, num_output_fields, size;
411   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr);
412   ierr = CeedOperatorGetNumElements(op, &num_elem); CeedChkBackend(ierr);
413   ierr= CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields);
414   CeedChkBackend(ierr);
415   CeedOperatorField *op_input_fields, *op_output_fields;
416   ierr = CeedOperatorGetFields(op, &op_input_fields, &op_output_fields);
417   CeedChkBackend(ierr);
418   CeedQFunctionField *qf_input_fields, *qf_output_fields;
419   ierr = CeedQFunctionGetFields(qf, &qf_input_fields, &qf_output_fields);
420   CeedChkBackend(ierr);
421   CeedEvalMode eval_mode;
422   CeedVector vec;
423   CeedElemRestriction elem_restr;
424 
425   // Setup
426   ierr = CeedOperatorSetup_Ref(op); CeedChkBackend(ierr);
427 
428   // Input Evecs and Restriction
429   ierr = CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields,
430                                      op_input_fields, in_vec, false, impl,
431                                      request); CeedChkBackend(ierr);
432 
433   // Output Evecs
434   for (CeedInt i=0; i<num_output_fields; i++) {
435     ierr = CeedVectorGetArray(impl->e_vecs[i+impl->num_e_vecs_in], CEED_MEM_HOST,
436                               &impl->e_data[i + num_input_fields]); CeedChkBackend(ierr);
437   }
438 
439   // Loop through elements
440   for (CeedInt e=0; e<num_elem; e++) {
441     // Output pointers
442     for (CeedInt i=0; i<num_output_fields; i++) {
443       ierr = CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode);
444       CeedChkBackend(ierr);
445       if (eval_mode == CEED_EVAL_NONE) {
446         ierr = CeedQFunctionFieldGetSize(qf_output_fields[i], &size);
447         CeedChkBackend(ierr);
448         ierr = CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST,
449                                   CEED_USE_POINTER,
450                                   &impl->e_data[i + num_input_fields][e*Q*size]);
451         CeedChkBackend(ierr);
452       }
453     }
454 
455     // Input basis apply
456     ierr = CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields,
457                                       num_input_fields, false, impl);
458     CeedChkBackend(ierr);
459 
460     // Q function
461     if (!impl->identity_qf) {
462       ierr = CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out);
463       CeedChkBackend(ierr);
464     }
465 
466     // Output basis apply
467     ierr = CeedOperatorOutputBasis_Ref(e, Q, qf_output_fields, op_output_fields,
468                                        num_input_fields, num_output_fields, op, impl);
469     CeedChkBackend(ierr);
470   }
471 
472   // Output restriction
473   for (CeedInt i=0; i<num_output_fields; i++) {
474     // Restore Evec
475     ierr = CeedVectorRestoreArray(impl->e_vecs[i+impl->num_e_vecs_in],
476                                   &impl->e_data[i + num_input_fields]);
477     CeedChkBackend(ierr);
478     // Get output vector
479     ierr = CeedOperatorFieldGetVector(op_output_fields[i], &vec);
480     CeedChkBackend(ierr);
481     // Active
482     if (vec == CEED_VECTOR_ACTIVE)
483       vec = out_vec;
484     // Restrict
485     ierr = CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_restr);
486     CeedChkBackend(ierr);
487     ierr = CeedElemRestrictionApply(elem_restr, CEED_TRANSPOSE,
488                                     impl->e_vecs[i+impl->num_e_vecs_in], vec, request);
489     CeedChkBackend(ierr);
490   }
491 
492   // Restore input arrays
493   ierr = CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields,
494                                        op_input_fields, false, impl);
495   CeedChkBackend(ierr);
496 
497   return CEED_ERROR_SUCCESS;
498 }
499 
500 //------------------------------------------------------------------------------
501 // Assemble Linear QFunction
502 //------------------------------------------------------------------------------
503 static int CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op,
504     CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) {
505   int ierr;
506   CeedOperator_Ref *impl;
507   ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr);
508   CeedQFunction qf;
509   ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr);
510   CeedInt Q, num_elem, num_input_fields, num_output_fields, size;
511   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr);
512   ierr = CeedOperatorGetNumElements(op, &num_elem); CeedChkBackend(ierr);
513   ierr= CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields);
514   CeedChkBackend(ierr);
515   CeedOperatorField *op_input_fields, *op_output_fields;
516   ierr = CeedOperatorGetFields(op, &op_input_fields, &op_output_fields);
517   CeedChkBackend(ierr);
518   CeedQFunctionField *qf_input_fields, *qf_output_fields;
519   ierr = CeedQFunctionGetFields(qf, &qf_input_fields, &qf_output_fields);
520   CeedChkBackend(ierr);
521   CeedVector vec;
522   CeedInt num_active_in = 0, num_active_out = 0;
523   CeedVector *active_in = NULL;
524   CeedScalar *a, *tmp;
525   Ceed ceed, ceed_parent;
526   ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr);
527   ierr = CeedGetOperatorFallbackParentCeed(ceed, &ceed_parent);
528   CeedChkBackend(ierr);
529   ceed_parent = ceed_parent ? ceed_parent : ceed;
530 
531   // Setup
532   ierr = CeedOperatorSetup_Ref(op); CeedChkBackend(ierr);
533 
534   // Check for identity
535   if (impl->identity_qf)
536     // LCOV_EXCL_START
537     return CeedError(ceed, CEED_ERROR_BACKEND,
538                      "Assembling identity QFunctions not supported");
539   // LCOV_EXCL_STOP
540 
541   // Input Evecs and Restriction
542   ierr = CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields,
543                                      op_input_fields, NULL, true, impl, request);
544   CeedChkBackend(ierr);
545 
546   // Count number of active input fields
547   for (CeedInt i=0; i<num_input_fields; i++) {
548     // Get input vector
549     ierr = CeedOperatorFieldGetVector(op_input_fields[i], &vec);
550     CeedChkBackend(ierr);
551     // Check if active input
552     if (vec == CEED_VECTOR_ACTIVE) {
553       ierr = CeedQFunctionFieldGetSize(qf_input_fields[i], &size);
554       CeedChkBackend(ierr);
555       ierr = CeedVectorSetValue(impl->q_vecs_in[i], 0.0); CeedChkBackend(ierr);
556       ierr = CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &tmp);
557       CeedChkBackend(ierr);
558       ierr = CeedRealloc(num_active_in + size, &active_in); CeedChkBackend(ierr);
559       for (CeedInt field=0; field<size; field++) {
560         ierr = CeedVectorCreate(ceed, Q, &active_in[num_active_in+field]);
561         CeedChkBackend(ierr);
562         ierr = CeedVectorSetArray(active_in[num_active_in+field], CEED_MEM_HOST,
563                                   CEED_USE_POINTER, &tmp[field*Q]);
564         CeedChkBackend(ierr);
565       }
566       num_active_in += size;
567       ierr = CeedVectorRestoreArray(impl->q_vecs_in[i], &tmp); CeedChkBackend(ierr);
568     }
569   }
570 
571   // Count number of active output fields
572   for (CeedInt i=0; i<num_output_fields; i++) {
573     // Get output vector
574     ierr = CeedOperatorFieldGetVector(op_output_fields[i], &vec);
575     CeedChkBackend(ierr);
576     // Check if active output
577     if (vec == CEED_VECTOR_ACTIVE) {
578       ierr = CeedQFunctionFieldGetSize(qf_output_fields[i], &size);
579       CeedChkBackend(ierr);
580       num_active_out += size;
581     }
582   }
583 
584   // Check sizes
585   if (!num_active_in || !num_active_out)
586     // LCOV_EXCL_START
587     return CeedError(ceed, CEED_ERROR_BACKEND,
588                      "Cannot assemble QFunction without active inputs "
589                      "and outputs");
590   // LCOV_EXCL_STOP
591 
592   // Create output restriction
593   CeedInt strides[3] = {1, Q, num_active_in*num_active_out*Q}; /* *NOPAD* */
594   ierr = CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q,
595                                           num_active_in*num_active_out,
596                                           num_active_in*num_active_out*num_elem*Q,
597                                           strides, rstr); CeedChkBackend(ierr);
598   // Create assembled vector
599   ierr = CeedVectorCreate(ceed_parent, num_elem*Q*num_active_in*num_active_out,
600                           assembled); CeedChkBackend(ierr);
601   ierr = CeedVectorSetValue(*assembled, 0.0); CeedChkBackend(ierr);
602   ierr = CeedVectorGetArray(*assembled, CEED_MEM_HOST, &a); CeedChkBackend(ierr);
603 
604   // Loop through elements
605   for (CeedInt e=0; e<num_elem; e++) {
606     // Input basis apply
607     ierr = CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields,
608                                       num_input_fields, true, impl);
609     CeedChkBackend(ierr);
610 
611     // Assemble QFunction
612     for (CeedInt in=0; in<num_active_in; in++) {
613       // Set Inputs
614       ierr = CeedVectorSetValue(active_in[in], 1.0); CeedChkBackend(ierr);
615       if (num_active_in > 1) {
616         ierr = CeedVectorSetValue(active_in[(in+num_active_in-1)%num_active_in],
617                                   0.0); CeedChkBackend(ierr);
618       }
619       // Set Outputs
620       for (CeedInt out=0; out<num_output_fields; out++) {
621         // Get output vector
622         ierr = CeedOperatorFieldGetVector(op_output_fields[out], &vec);
623         CeedChkBackend(ierr);
624         // Check if active output
625         if (vec == CEED_VECTOR_ACTIVE) {
626           CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST,
627                              CEED_USE_POINTER, a); CeedChkBackend(ierr);
628           ierr = CeedQFunctionFieldGetSize(qf_output_fields[out], &size);
629           CeedChkBackend(ierr);
630           a += size*Q; // Advance the pointer by the size of the output
631         }
632       }
633       // Apply QFunction
634       ierr = CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out);
635       CeedChkBackend(ierr);
636     }
637   }
638 
639   // Un-set output Qvecs to prevent accidental overwrite of Assembled
640   for (CeedInt out=0; out<num_output_fields; out++) {
641     // Get output vector
642     ierr = CeedOperatorFieldGetVector(op_output_fields[out], &vec);
643     CeedChkBackend(ierr);
644     // Check if active output
645     if (vec == CEED_VECTOR_ACTIVE) {
646       CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL);
647       CeedChkBackend(ierr);
648     }
649   }
650 
651   // Restore input arrays
652   ierr = CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields,
653                                        op_input_fields, true, impl);
654   CeedChkBackend(ierr);
655 
656   // Restore output
657   ierr = CeedVectorRestoreArray(*assembled, &a); CeedChkBackend(ierr);
658 
659   // Cleanup
660   for (CeedInt i=0; i<num_active_in; i++) {
661     ierr = CeedVectorDestroy(&active_in[i]); CeedChkBackend(ierr);
662   }
663   ierr = CeedFree(&active_in); CeedChkBackend(ierr);
664 
665   return CEED_ERROR_SUCCESS;
666 }
667 
668 //------------------------------------------------------------------------------
669 // Get Basis Emode Pointer
670 //------------------------------------------------------------------------------
671 static inline void CeedOperatorGetBasisPointer_Ref(const CeedScalar **basis_ptr,
672     CeedEvalMode eval_mode, const CeedScalar *identity, const CeedScalar *interp,
673     const CeedScalar *grad) {
674   switch (eval_mode) {
675   case CEED_EVAL_NONE:
676     *basis_ptr = identity;
677     break;
678   case CEED_EVAL_INTERP:
679     *basis_ptr = interp;
680     break;
681   case CEED_EVAL_GRAD:
682     *basis_ptr = grad;
683     break;
684   case CEED_EVAL_WEIGHT:
685   case CEED_EVAL_DIV:
686   case CEED_EVAL_CURL:
687     break; // Caught by QF Assembly
688   }
689 }
690 
691 //------------------------------------------------------------------------------
692 // Create point block restriction
693 //------------------------------------------------------------------------------
694 static int CreatePBRestriction_Ref(CeedElemRestriction rstr,
695                                    CeedElemRestriction *pb_rstr) {
696   int ierr;
697   Ceed ceed;
698   ierr = CeedElemRestrictionGetCeed(rstr, &ceed); CeedChkBackend(ierr);
699   const CeedInt *offsets;
700   ierr = CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets);
701   CeedChkBackend(ierr);
702 
703   // Expand offsets
704   CeedInt num_elem, num_comp, elem_size, comp_stride, max = 1, *pbOffsets;
705   ierr = CeedElemRestrictionGetNumElements(rstr, &num_elem); CeedChkBackend(ierr);
706   ierr = CeedElemRestrictionGetNumComponents(rstr, &num_comp);
707   CeedChkBackend(ierr);
708   ierr = CeedElemRestrictionGetElementSize(rstr, &elem_size);
709   CeedChkBackend(ierr);
710   ierr = CeedElemRestrictionGetCompStride(rstr, &comp_stride);
711   CeedChkBackend(ierr);
712   CeedInt shift = num_comp;
713   if (comp_stride != 1)
714     shift *= num_comp;
715   ierr = CeedCalloc(num_elem*elem_size, &pbOffsets); CeedChkBackend(ierr);
716   for (CeedInt i = 0; i < num_elem*elem_size; i++) {
717     pbOffsets[i] = offsets[i]*shift;
718     if (pbOffsets[i] > max)
719       max = pbOffsets[i];
720   }
721 
722   // Create new restriction
723   ierr = CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp*num_comp,
724                                    1,
725                                    max + num_comp*num_comp, CEED_MEM_HOST,
726                                    CEED_OWN_POINTER, pbOffsets, pb_rstr);
727   CeedChkBackend(ierr);
728 
729   // Cleanup
730   ierr = CeedElemRestrictionRestoreOffsets(rstr, &offsets); CeedChkBackend(ierr);
731 
732   return CEED_ERROR_SUCCESS;
733 }
734 
735 //------------------------------------------------------------------------------
736 // Assemble diagonal common code
737 //------------------------------------------------------------------------------
738 static inline int CeedOperatorAssembleAddDiagonalCore_Ref(CeedOperator op,
739     CeedVector assembled, CeedRequest *request, const bool is_pointblock) {
740   int ierr;
741   Ceed ceed;
742   ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr);
743 
744   // Assemble QFunction
745   CeedQFunction qf;
746   ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr);
747   CeedInt num_input_fields, num_output_fields;
748   ierr= CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields);
749   CeedChkBackend(ierr);
750   CeedVector assembled_qf;
751   CeedElemRestriction rstr;
752   ierr = CeedOperatorLinearAssembleQFunction(op,  &assembled_qf, &rstr, request);
753   CeedChkBackend(ierr);
754   ierr = CeedElemRestrictionDestroy(&rstr); CeedChkBackend(ierr);
755   CeedScalar max_norm = 0;
756   ierr = CeedVectorNorm(assembled_qf, CEED_NORM_MAX, &max_norm);
757   CeedChkBackend(ierr);
758 
759   // Determine active input basis
760   CeedOperatorField *op_fields;
761   CeedQFunctionField *qf_fields;
762   ierr = CeedOperatorGetFields(op, &op_fields, NULL); CeedChkBackend(ierr);
763   ierr = CeedQFunctionGetFields(qf, &qf_fields, NULL); CeedChkBackend(ierr);
764   CeedInt num_eval_mode_in = 0, num_comp, dim = 1;
765   CeedEvalMode *eval_mode_in = NULL;
766   CeedBasis basis_in = NULL;
767   CeedElemRestriction rstr_in = NULL;
768   for (CeedInt i=0; i<num_input_fields; i++) {
769     CeedVector vec;
770     ierr = CeedOperatorFieldGetVector(op_fields[i], &vec); CeedChkBackend(ierr);
771     if (vec == CEED_VECTOR_ACTIVE) {
772       CeedElemRestriction rstr;
773       ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis_in); CeedChkBackend(ierr);
774       ierr = CeedBasisGetNumComponents(basis_in, &num_comp); CeedChkBackend(ierr);
775       ierr = CeedBasisGetDimension(basis_in, &dim); CeedChkBackend(ierr);
776       ierr = CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr);
777       CeedChkBackend(ierr);
778       if (rstr_in && rstr_in != rstr)
779         // LCOV_EXCL_START
780         return CeedError(ceed, CEED_ERROR_BACKEND,
781                          "Multi-field non-composite operator diagonal assembly not supported");
782       // LCOV_EXCL_STOP
783       rstr_in = rstr;
784       CeedEvalMode eval_mode;
785       ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode);
786       CeedChkBackend(ierr);
787       switch (eval_mode) {
788       case CEED_EVAL_NONE:
789       case CEED_EVAL_INTERP:
790         ierr = CeedRealloc(num_eval_mode_in + 1, &eval_mode_in); CeedChkBackend(ierr);
791         eval_mode_in[num_eval_mode_in] = eval_mode;
792         num_eval_mode_in += 1;
793         break;
794       case CEED_EVAL_GRAD:
795         ierr = CeedRealloc(num_eval_mode_in + dim, &eval_mode_in); CeedChkBackend(ierr);
796         for (CeedInt d=0; d<dim; d++)
797           eval_mode_in[num_eval_mode_in+d] = eval_mode;
798         num_eval_mode_in += dim;
799         break;
800       case CEED_EVAL_WEIGHT:
801       case CEED_EVAL_DIV:
802       case CEED_EVAL_CURL:
803         break; // Caught by QF Assembly
804       }
805     }
806   }
807 
808   // Determine active output basis
809   ierr = CeedOperatorGetFields(op, NULL, &op_fields); CeedChkBackend(ierr);
810   ierr = CeedQFunctionGetFields(qf, NULL, &qf_fields); CeedChkBackend(ierr);
811   CeedInt num_eval_mode_out = 0;
812   CeedEvalMode *eval_mode_out = NULL;
813   CeedBasis basis_out = NULL;
814   CeedElemRestriction rstr_out = NULL;
815   for (CeedInt i=0; i<num_output_fields; i++) {
816     CeedVector vec;
817     ierr = CeedOperatorFieldGetVector(op_fields[i], &vec); CeedChkBackend(ierr);
818     if (vec == CEED_VECTOR_ACTIVE) {
819       CeedElemRestriction rstr;
820       ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis_out);
821       CeedChkBackend(ierr);
822       ierr = CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr);
823       CeedChkBackend(ierr);
824       if (rstr_out && rstr_out != rstr)
825         // LCOV_EXCL_START
826         return CeedError(ceed, CEED_ERROR_BACKEND,
827                          "Multi-field non-composite operator diagonal assembly not supported");
828       // LCOV_EXCL_STOP
829       rstr_out = rstr;
830       CeedEvalMode eval_mode;
831       ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode);
832       CeedChkBackend(ierr);
833       switch (eval_mode) {
834       case CEED_EVAL_NONE:
835       case CEED_EVAL_INTERP:
836         ierr = CeedRealloc(num_eval_mode_out + 1, &eval_mode_out); CeedChkBackend(ierr);
837         eval_mode_out[num_eval_mode_out] = eval_mode;
838         num_eval_mode_out += 1;
839         break;
840       case CEED_EVAL_GRAD:
841         ierr = CeedRealloc(num_eval_mode_out + dim, &eval_mode_out);
842         CeedChkBackend(ierr);
843         for (CeedInt d=0; d<dim; d++)
844           eval_mode_out[num_eval_mode_out+d] = eval_mode;
845         num_eval_mode_out += dim;
846         break;
847       case CEED_EVAL_WEIGHT:
848       case CEED_EVAL_DIV:
849       case CEED_EVAL_CURL:
850         break; // Caught by QF Assembly
851       }
852     }
853   }
854 
855   // Assemble point-block diagonal restriction, if needed
856   CeedElemRestriction diag_rstr = rstr_out;
857   if (is_pointblock) {
858     ierr = CreatePBRestriction_Ref(rstr_out, &diag_rstr); CeedChkBackend(ierr);
859   }
860 
861   // Create diagonal vector
862   CeedVector elem_diag;
863   ierr = CeedElemRestrictionCreateVector(diag_rstr, NULL, &elem_diag);
864   CeedChkBackend(ierr);
865 
866   // Assemble element operator diagonals
867   CeedScalar *elem_diag_array, *assembled_qf_array;
868   ierr = CeedVectorSetValue(elem_diag, 0.0); CeedChkBackend(ierr);
869   ierr = CeedVectorGetArray(elem_diag, CEED_MEM_HOST, &elem_diag_array);
870   CeedChkBackend(ierr);
871   ierr = CeedVectorGetArray(assembled_qf, CEED_MEM_HOST, &assembled_qf_array);
872   CeedChkBackend(ierr);
873   CeedInt num_elem, num_nodes, num_qpts;
874   ierr = CeedElemRestrictionGetNumElements(diag_rstr, &num_elem);
875   CeedChkBackend(ierr);
876   ierr = CeedBasisGetNumNodes(basis_in, &num_nodes); CeedChkBackend(ierr);
877   ierr = CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts);
878   CeedChkBackend(ierr);
879   // Basis matrices
880   const CeedScalar *interp_in, *interp_out, *grad_in, *grad_out;
881   CeedScalar *identity = NULL;
882   bool evalNone = false;
883   for (CeedInt i=0; i<num_eval_mode_in; i++)
884     evalNone = evalNone || (eval_mode_in[i] == CEED_EVAL_NONE);
885   for (CeedInt i=0; i<num_eval_mode_out; i++)
886     evalNone = evalNone || (eval_mode_out[i] == CEED_EVAL_NONE);
887   if (evalNone) {
888     ierr = CeedCalloc(num_qpts*num_nodes, &identity); CeedChkBackend(ierr);
889     for (CeedInt i=0; i<(num_nodes<num_qpts?num_nodes:num_qpts); i++)
890       identity[i*num_nodes+i] = 1.0;
891   }
892   ierr = CeedBasisGetInterp(basis_in, &interp_in); CeedChkBackend(ierr);
893   ierr = CeedBasisGetInterp(basis_out, &interp_out); CeedChkBackend(ierr);
894   ierr = CeedBasisGetGrad(basis_in, &grad_in); CeedChkBackend(ierr);
895   ierr = CeedBasisGetGrad(basis_out, &grad_out); CeedChkBackend(ierr);
896   // Compute the diagonal of B^T D B
897   // Each element
898   const CeedScalar qf_value_bound = max_norm*1e-12;
899   for (CeedInt e=0; e<num_elem; e++) {
900     CeedInt d_out = -1;
901     // Each basis eval mode pair
902     for (CeedInt e_out=0; e_out<num_eval_mode_out; e_out++) {
903       const CeedScalar *bt = NULL;
904       if (eval_mode_out[e_out] == CEED_EVAL_GRAD)
905         d_out += 1;
906       CeedOperatorGetBasisPointer_Ref(&bt, eval_mode_out[e_out], identity, interp_out,
907                                       &grad_out[d_out*num_qpts*num_nodes]);
908       CeedInt d_in = -1;
909       for (CeedInt e_in=0; e_in<num_eval_mode_in; e_in++) {
910         const CeedScalar *b = NULL;
911         if (eval_mode_in[e_in] == CEED_EVAL_GRAD)
912           d_in += 1;
913         CeedOperatorGetBasisPointer_Ref(&b, eval_mode_in[e_in], identity, interp_in,
914                                         &grad_in[d_in*num_qpts*num_nodes]);
915         // Each component
916         for (CeedInt c_out=0; c_out<num_comp; c_out++)
917           // Each qpoint/node pair
918           for (CeedInt q=0; q<num_qpts; q++)
919             if (is_pointblock) {
920               // Point Block Diagonal
921               for (CeedInt c_in=0; c_in<num_comp; c_in++) {
922                 const CeedScalar qf_value =
923                   assembled_qf_array[((((e*num_eval_mode_in+e_in)*num_comp+c_in)*
924                                        num_eval_mode_out+e_out)*num_comp+c_out)*num_qpts+q];
925                 if (fabs(qf_value) > qf_value_bound)
926                   for (CeedInt n=0; n<num_nodes; n++)
927                     elem_diag_array[((e*num_comp+c_out)*num_comp+c_in)*num_nodes+n] +=
928                       bt[q*num_nodes+n] * qf_value * b[q*num_nodes+n];
929               }
930             } else {
931               // Diagonal Only
932               const CeedScalar qf_value =
933                 assembled_qf_array[((((e*num_eval_mode_in+e_in)*num_comp+c_out)*
934                                      num_eval_mode_out+e_out)*num_comp+c_out)*num_qpts+q];
935               if (fabs(qf_value) > qf_value_bound)
936                 for (CeedInt n=0; n<num_nodes; n++)
937                   elem_diag_array[(e*num_comp+c_out)*num_nodes+n] +=
938                     bt[q*num_nodes+n] * qf_value * b[q*num_nodes+n];
939             }
940       }
941     }
942   }
943   ierr = CeedVectorRestoreArray(elem_diag, &elem_diag_array);
944   CeedChkBackend(ierr);
945   ierr = CeedVectorRestoreArray(assembled_qf, &assembled_qf_array);
946   CeedChkBackend(ierr);
947 
948   // Assemble local operator diagonal
949   ierr = CeedElemRestrictionApply(diag_rstr, CEED_TRANSPOSE, elem_diag,
950                                   assembled, request); CeedChkBackend(ierr);
951 
952   // Cleanup
953   if (is_pointblock) {
954     ierr = CeedElemRestrictionDestroy(&diag_rstr); CeedChkBackend(ierr);
955   }
956   ierr = CeedVectorDestroy(&assembled_qf); CeedChkBackend(ierr);
957   ierr = CeedVectorDestroy(&elem_diag); CeedChkBackend(ierr);
958   ierr = CeedFree(&eval_mode_in); CeedChkBackend(ierr);
959   ierr = CeedFree(&eval_mode_out); CeedChkBackend(ierr);
960   ierr = CeedFree(&identity); CeedChkBackend(ierr);
961 
962   return CEED_ERROR_SUCCESS;
963 }
964 
965 //------------------------------------------------------------------------------
966 // Assemble composite diagonal common code
967 //------------------------------------------------------------------------------
968 static inline int CeedOperatorLinearAssembleAddDiagonalCompositeCore_Ref(
969   CeedOperator op, CeedVector assembled, CeedRequest *request,
970   const bool is_pointblock) {
971   int ierr;
972   CeedInt num_sub;
973   CeedOperator *suboperators;
974   ierr = CeedOperatorGetNumSub(op, &num_sub); CeedChkBackend(ierr);
975   ierr = CeedOperatorGetSubList(op, &suboperators); CeedChkBackend(ierr);
976   for (CeedInt i = 0; i < num_sub; i++) {
977     ierr = CeedOperatorAssembleAddDiagonalCore_Ref(suboperators[i], assembled,
978            request, is_pointblock); CeedChkBackend(ierr);
979   }
980   return CEED_ERROR_SUCCESS;
981 }
982 
983 //------------------------------------------------------------------------------
984 // Assemble Linear Diagonal
985 //------------------------------------------------------------------------------
986 static int CeedOperatorLinearAssembleAddDiagonal_Ref(CeedOperator op,
987     CeedVector assembled, CeedRequest *request) {
988   int ierr;
989   bool is_composite;
990   ierr = CeedOperatorIsComposite(op, &is_composite); CeedChkBackend(ierr);
991   if (is_composite) {
992     return CeedOperatorLinearAssembleAddDiagonalCompositeCore_Ref(op, assembled,
993            request, false);
994   } else {
995     return CeedOperatorAssembleAddDiagonalCore_Ref(op, assembled, request, false);
996   }
997 }
998 
999 //------------------------------------------------------------------------------
1000 // Assemble Linear Point Block Diagonal
1001 //------------------------------------------------------------------------------
1002 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Ref(CeedOperator op,
1003     CeedVector assembled, CeedRequest *request) {
1004   int ierr;
1005   bool is_composite;
1006   ierr = CeedOperatorIsComposite(op, &is_composite); CeedChkBackend(ierr);
1007   if (is_composite) {
1008     return CeedOperatorLinearAssembleAddDiagonalCompositeCore_Ref(op, assembled,
1009            request, true);
1010   } else {
1011     return CeedOperatorAssembleAddDiagonalCore_Ref(op, assembled, request, true);
1012   }
1013 }
1014 
1015 //------------------------------------------------------------------------------
1016 // Create FDM Element Inverse
1017 //------------------------------------------------------------------------------
1018 int CeedOperatorCreateFDMElementInverse_Ref(CeedOperator op,
1019     CeedOperator *fdm_inv, CeedRequest *request) {
1020   int ierr;
1021   Ceed ceed, ceed_parent;
1022   ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr);
1023   ierr = CeedGetOperatorFallbackParentCeed(ceed, &ceed_parent);
1024   CeedChkBackend(ierr);
1025   ceed_parent = ceed_parent ? ceed_parent : ceed;
1026   CeedQFunction qf;
1027   ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr);
1028 
1029   // Determine active input basis
1030   bool interp = false, grad = false;
1031   CeedBasis basis = NULL;
1032   CeedElemRestriction rstr = NULL;
1033   CeedOperatorField *op_fields;
1034   CeedQFunctionField *qf_fields;
1035   ierr = CeedOperatorGetFields(op, &op_fields, NULL); CeedChkBackend(ierr);
1036   ierr = CeedQFunctionGetFields(qf, &qf_fields, NULL); CeedChkBackend(ierr);
1037   CeedInt num_input_fields;
1038   ierr = CeedQFunctionGetNumArgs(qf, &num_input_fields, NULL);
1039   CeedChkBackend(ierr);
1040   for (CeedInt i=0; i<num_input_fields; i++) {
1041     CeedVector vec;
1042     ierr = CeedOperatorFieldGetVector(op_fields[i], &vec); CeedChkBackend(ierr);
1043     if (vec == CEED_VECTOR_ACTIVE) {
1044       CeedEvalMode eval_mode;
1045       ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode);
1046       CeedChkBackend(ierr);
1047       interp = interp || eval_mode == CEED_EVAL_INTERP;
1048       grad = grad || eval_mode == CEED_EVAL_GRAD;
1049       ierr = CeedOperatorFieldGetBasis(op_fields[i], &basis); CeedChkBackend(ierr);
1050       ierr = CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr);
1051       CeedChkBackend(ierr);
1052     }
1053   }
1054   if (!basis)
1055     // LCOV_EXCL_START
1056     return CeedError(ceed, CEED_ERROR_BACKEND, "No active field set");
1057   // LCOV_EXCL_STOP
1058   CeedInt P_1d, Q_1d, elem_size, num_qpts, dim, num_comp = 1, num_elem = 1,
1059                                                 l_size = 1;
1060   ierr = CeedBasisGetNumNodes1D(basis, &P_1d); CeedChkBackend(ierr);
1061   ierr = CeedBasisGetNumNodes(basis, &elem_size); CeedChkBackend(ierr);
1062   ierr = CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d); CeedChkBackend(ierr);
1063   ierr = CeedBasisGetNumQuadraturePoints(basis, &num_qpts); CeedChkBackend(ierr);
1064   ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr);
1065   ierr = CeedBasisGetNumComponents(basis, &num_comp); CeedChkBackend(ierr);
1066   ierr = CeedElemRestrictionGetNumElements(rstr, &num_elem); CeedChkBackend(ierr);
1067   ierr = CeedElemRestrictionGetLVectorSize(rstr, &l_size); CeedChkBackend(ierr);
1068 
1069   // Build and diagonalize 1D Mass and Laplacian
1070   bool tensor_basis;
1071   ierr = CeedBasisIsTensor(basis, &tensor_basis); CeedChkBackend(ierr);
1072   if (!tensor_basis)
1073     // LCOV_EXCL_START
1074     return CeedError(ceed, CEED_ERROR_BACKEND,
1075                      "FDMElementInverse only supported for tensor "
1076                      "bases");
1077   // LCOV_EXCL_STOP
1078   CeedScalar *work, *mass, *laplace, *x, *x2, *lambda;
1079   ierr = CeedMalloc(Q_1d*P_1d, &work); CeedChkBackend(ierr);
1080   ierr = CeedMalloc(P_1d*P_1d, &mass); CeedChkBackend(ierr);
1081   ierr = CeedMalloc(P_1d*P_1d, &laplace); CeedChkBackend(ierr);
1082   ierr = CeedMalloc(P_1d*P_1d, &x); CeedChkBackend(ierr);
1083   ierr = CeedMalloc(P_1d*P_1d, &x2); CeedChkBackend(ierr);
1084   ierr = CeedMalloc(P_1d, &lambda); CeedChkBackend(ierr);
1085   // -- Mass
1086   const CeedScalar *interp_1d, *grad_1d, *q_weight_1d;
1087   ierr = CeedBasisGetInterp1D(basis, &interp_1d); CeedChkBackend(ierr);
1088   ierr = CeedBasisGetGrad1D(basis, &grad_1d); CeedChkBackend(ierr);
1089   ierr = CeedBasisGetQWeights(basis, &q_weight_1d); CeedChkBackend(ierr);
1090   for (CeedInt i=0; i<Q_1d; i++)
1091     for (CeedInt j=0; j<P_1d; j++)
1092       work[i+j*Q_1d] = interp_1d[i*P_1d+j]*q_weight_1d[i];
1093   ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)work,
1094                             (const CeedScalar *)interp_1d, mass, P_1d, P_1d, Q_1d);
1095   CeedChkBackend(ierr);
1096   // -- Laplacian
1097   for (CeedInt i=0; i<Q_1d; i++)
1098     for (CeedInt j=0; j<P_1d; j++)
1099       work[i+j*Q_1d] = grad_1d[i*P_1d+j]*q_weight_1d[i];
1100   ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)work,
1101                             (const CeedScalar *)grad_1d, laplace, P_1d, P_1d, Q_1d);
1102   CeedChkBackend(ierr);
1103   // -- Diagonalize
1104   ierr = CeedSimultaneousDiagonalization(ceed, laplace, mass, x, lambda, P_1d);
1105   CeedChkBackend(ierr);
1106   ierr = CeedFree(&work); CeedChkBackend(ierr);
1107   ierr = CeedFree(&mass); CeedChkBackend(ierr);
1108   ierr = CeedFree(&laplace); CeedChkBackend(ierr);
1109   for (CeedInt i=0; i<P_1d; i++)
1110     for (CeedInt j=0; j<P_1d; j++)
1111       x2[i+j*P_1d] = x[j+i*P_1d];
1112   ierr = CeedFree(&x); CeedChkBackend(ierr);
1113 
1114   // Assemble QFunction
1115   CeedVector assembled;
1116   CeedElemRestriction rstr_qf;
1117   ierr =  CeedOperatorLinearAssembleQFunction(op, &assembled, &rstr_qf,
1118           request); CeedChkBackend(ierr);
1119   ierr = CeedElemRestrictionDestroy(&rstr_qf); CeedChkBackend(ierr);
1120   CeedScalar max_norm = 0;
1121   ierr = CeedVectorNorm(assembled, CEED_NORM_MAX, &max_norm);
1122   CeedChkBackend(ierr);
1123 
1124   // Calculate element averages
1125   CeedInt num_fields = ((interp?1:0) + (grad?dim:0))*((interp?1:0) +
1126                        (grad?dim:0));
1127   CeedScalar *elem_avg;
1128   const CeedScalar *assembledarray, *q_weight_array;
1129   CeedVector q_weight;
1130   ierr = CeedVectorCreate(ceed_parent, num_qpts, &q_weight); CeedChkBackend(ierr);
1131   ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT,
1132                         CEED_VECTOR_NONE, q_weight); CeedChkBackend(ierr);
1133   ierr = CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembledarray);
1134   CeedChkBackend(ierr);
1135   ierr = CeedVectorGetArrayRead(q_weight, CEED_MEM_HOST, &q_weight_array);
1136   CeedChkBackend(ierr);
1137   ierr = CeedCalloc(num_elem, &elem_avg); CeedChkBackend(ierr);
1138   for (CeedInt e=0; e<num_elem; e++) {
1139     CeedInt count = 0;
1140     for (CeedInt q=0; q<num_qpts; q++)
1141       for (CeedInt i=0; i<num_comp*num_comp*num_fields; i++)
1142         if (fabs(assembledarray[e*num_elem*num_qpts*num_comp*num_comp*num_fields +
1143                                                                                  i*num_qpts + q]) > max_norm*1e-12) {
1144           elem_avg[e] += assembledarray[e*num_elem*num_qpts*num_comp*num_comp*num_fields +
1145                                         i*num_qpts + q] / q_weight_array[q];
1146           count++;
1147         }
1148     if (count)
1149       elem_avg[e] /= count;
1150   }
1151   ierr = CeedVectorRestoreArrayRead(assembled, &assembledarray);
1152   CeedChkBackend(ierr);
1153   ierr = CeedVectorDestroy(&assembled); CeedChkBackend(ierr);
1154   ierr = CeedVectorRestoreArrayRead(q_weight, &q_weight_array);
1155   CeedChkBackend(ierr);
1156   ierr = CeedVectorDestroy(&q_weight); CeedChkBackend(ierr);
1157 
1158   // Build FDM diagonal
1159   CeedVector q_data;
1160   CeedScalar *q_data_array;
1161   ierr = CeedVectorCreate(ceed_parent, num_elem*num_comp*l_size, &q_data);
1162   CeedChkBackend(ierr);
1163   ierr = CeedVectorSetArray(q_data, CEED_MEM_HOST, CEED_COPY_VALUES, NULL);
1164   CeedChkBackend(ierr);
1165   ierr = CeedVectorGetArray(q_data, CEED_MEM_HOST, &q_data_array);
1166   CeedChkBackend(ierr);
1167   for (CeedInt e=0; e<num_elem; e++)
1168     for (CeedInt c=0; c<num_comp; c++)
1169       for (CeedInt n=0; n<l_size; n++) {
1170         if (interp)
1171           q_data_array[(e*num_comp+c)*l_size+n] = 1;
1172         if (grad)
1173           for (CeedInt d=0; d<dim; d++) {
1174             CeedInt i = (n / CeedIntPow(P_1d, d)) % P_1d;
1175             q_data_array[(e*num_comp+c)*l_size+n] += lambda[i];
1176           }
1177         q_data_array[(e*num_comp+c)*l_size+n] = 1 / (elem_avg[e] *
1178                                                 q_data_array[(e*num_comp+c)*l_size+n]);
1179       }
1180   ierr = CeedFree(&elem_avg); CeedChkBackend(ierr);
1181   ierr = CeedVectorRestoreArray(q_data, &q_data_array); CeedChkBackend(ierr);
1182 
1183   // Setup FDM operator
1184   // -- Basis
1185   CeedBasis fdm_basis;
1186   CeedScalar *grad_dummy, *q_ref_dummy, *q_weight_dummy;
1187   ierr = CeedCalloc(P_1d*P_1d, &grad_dummy); CeedChkBackend(ierr);
1188   ierr = CeedCalloc(P_1d, &q_ref_dummy); CeedChkBackend(ierr);
1189   ierr = CeedCalloc(P_1d, &q_weight_dummy); CeedChkBackend(ierr);
1190   ierr = CeedBasisCreateTensorH1(ceed_parent, dim, num_comp, P_1d, P_1d, x2,
1191                                  grad_dummy, q_ref_dummy, q_weight_dummy,
1192                                  &fdm_basis); CeedChkBackend(ierr);
1193   ierr = CeedFree(&grad_dummy); CeedChkBackend(ierr);
1194   ierr = CeedFree(&q_ref_dummy); CeedChkBackend(ierr);
1195   ierr = CeedFree(&q_weight_dummy); CeedChkBackend(ierr);
1196   ierr = CeedFree(&x2); CeedChkBackend(ierr);
1197   ierr = CeedFree(&lambda); CeedChkBackend(ierr);
1198 
1199   // -- Restriction
1200   CeedElemRestriction rstr_i;
1201   CeedInt strides[3] = {1, l_size, l_size*num_comp};
1202   ierr = CeedElemRestrictionCreateStrided(ceed_parent, num_elem, l_size, num_comp,
1203                                           l_size*num_elem*num_comp, strides, &rstr_i);
1204   CeedChkBackend(ierr);
1205   // -- QFunction
1206   CeedQFunction mass_qf;
1207   ierr = CeedQFunctionCreateInteriorByName(ceed_parent, "MassApply", &mass_qf);
1208   CeedChkBackend(ierr);
1209   // -- Operator
1210   ierr = CeedOperatorCreate(ceed_parent, mass_qf, NULL, NULL, fdm_inv);
1211   CeedChkBackend(ierr);
1212   CeedOperatorSetField(*fdm_inv, "u", rstr_i, fdm_basis, CEED_VECTOR_ACTIVE);
1213   CeedChkBackend(ierr);
1214   CeedOperatorSetField(*fdm_inv, "qdata", rstr_i, CEED_BASIS_COLLOCATED, q_data);
1215   CeedChkBackend(ierr);
1216   CeedOperatorSetField(*fdm_inv, "v", rstr_i, fdm_basis, CEED_VECTOR_ACTIVE);
1217   CeedChkBackend(ierr);
1218 
1219   // Cleanup
1220   ierr = CeedVectorDestroy(&q_data); CeedChkBackend(ierr);
1221   ierr = CeedBasisDestroy(&fdm_basis); CeedChkBackend(ierr);
1222   ierr = CeedElemRestrictionDestroy(&rstr_i); CeedChkBackend(ierr);
1223   ierr = CeedQFunctionDestroy(&mass_qf); CeedChkBackend(ierr);
1224 
1225   return CEED_ERROR_SUCCESS;
1226 }
1227 
1228 //------------------------------------------------------------------------------
1229 // Operator Destroy
1230 //------------------------------------------------------------------------------
1231 static int CeedOperatorDestroy_Ref(CeedOperator op) {
1232   int ierr;
1233   CeedOperator_Ref *impl;
1234   ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr);
1235 
1236   for (CeedInt i=0; i<impl->num_e_vecs_in+impl->num_e_vecs_out; i++) {
1237     ierr = CeedVectorDestroy(&impl->e_vecs[i]); CeedChkBackend(ierr);
1238   }
1239   ierr = CeedFree(&impl->e_vecs); CeedChkBackend(ierr);
1240   ierr = CeedFree(&impl->e_data); CeedChkBackend(ierr);
1241   ierr = CeedFree(&impl->input_state); CeedChkBackend(ierr);
1242 
1243   for (CeedInt i=0; i<impl->num_e_vecs_in; i++) {
1244     ierr = CeedVectorDestroy(&impl->e_vecs_in[i]); CeedChkBackend(ierr);
1245     ierr = CeedVectorDestroy(&impl->q_vecs_in[i]); CeedChkBackend(ierr);
1246   }
1247   ierr = CeedFree(&impl->e_vecs_in); CeedChkBackend(ierr);
1248   ierr = CeedFree(&impl->q_vecs_in); CeedChkBackend(ierr);
1249 
1250   for (CeedInt i=0; i<impl->num_e_vecs_out; i++) {
1251     ierr = CeedVectorDestroy(&impl->e_vecs_out[i]); CeedChkBackend(ierr);
1252     ierr = CeedVectorDestroy(&impl->q_vecs_out[i]); CeedChkBackend(ierr);
1253   }
1254   ierr = CeedFree(&impl->e_vecs_out); CeedChkBackend(ierr);
1255   ierr = CeedFree(&impl->q_vecs_out); CeedChkBackend(ierr);
1256 
1257   ierr = CeedFree(&impl); CeedChkBackend(ierr);
1258   return CEED_ERROR_SUCCESS;
1259 }
1260 
1261 //------------------------------------------------------------------------------
1262 // Operator Create
1263 //------------------------------------------------------------------------------
1264 int CeedOperatorCreate_Ref(CeedOperator op) {
1265   int ierr;
1266   Ceed ceed;
1267   ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr);
1268   CeedOperator_Ref *impl;
1269 
1270   ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr);
1271   ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr);
1272 
1273   ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction",
1274                                 CeedOperatorLinearAssembleQFunction_Ref);
1275   CeedChkBackend(ierr);
1276   ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal",
1277                                 CeedOperatorLinearAssembleAddDiagonal_Ref);
1278   CeedChkBackend(ierr);
1279   ierr = CeedSetBackendFunction(ceed, "Operator", op,
1280                                 "LinearAssembleAddPointBlockDiagonal",
1281                                 CeedOperatorLinearAssembleAddPointBlockDiagonal_Ref);
1282   CeedChkBackend(ierr);
1283   ierr = CeedSetBackendFunction(ceed, "Operator", op, "CreateFDMElementInverse",
1284                                 CeedOperatorCreateFDMElementInverse_Ref);
1285   CeedChkBackend(ierr);
1286   ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd",
1287                                 CeedOperatorApplyAdd_Ref); CeedChkBackend(ierr);
1288   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy",
1289                                 CeedOperatorDestroy_Ref); CeedChkBackend(ierr);
1290   return CEED_ERROR_SUCCESS;
1291 }
1292 
1293 //------------------------------------------------------------------------------
1294 // Composite Operator Create
1295 //------------------------------------------------------------------------------
1296 int CeedCompositeOperatorCreate_Ref(CeedOperator op) {
1297   int ierr;
1298   Ceed ceed;
1299   ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr);
1300   ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal",
1301                                 CeedOperatorLinearAssembleAddDiagonal_Ref);
1302   CeedChkBackend(ierr);
1303   ierr = CeedSetBackendFunction(ceed, "Operator", op,
1304                                 "LinearAssembleAddPointBlockDiagonal",
1305                                 CeedOperatorLinearAssembleAddPointBlockDiagonal_Ref);
1306   CeedChkBackend(ierr);
1307   return CEED_ERROR_SUCCESS;
1308 }
1309 //------------------------------------------------------------------------------
1310