xref: /libCEED/backends/ref/ceed-ref-operator.c (revision 16911fdad6ca4b1dd37f9d3206958ee664667dbe) 
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-ref.h"
18 
19 static int CeedOperatorDestroy_Ref(CeedOperator op) {
20   int ierr;
21   CeedOperator_Ref *impl;
22   ierr = CeedOperatorGetData(op, (void *)&impl); CeedChk(ierr);
23 
24   for (CeedInt i=0; i<impl->numein+impl->numeout; i++) {
25     ierr = CeedVectorDestroy(&impl->evecs[i]); CeedChk(ierr);
26   }
27   ierr = CeedFree(&impl->evecs); CeedChk(ierr);
28   ierr = CeedFree(&impl->edata); CeedChk(ierr);
29   ierr = CeedFree(&impl->inputstate); CeedChk(ierr);
30 
31   for (CeedInt i=0; i<impl->numein; i++) {
32     ierr = CeedVectorDestroy(&impl->evecsin[i]); CeedChk(ierr);
33     ierr = CeedVectorDestroy(&impl->qvecsin[i]); CeedChk(ierr);
34   }
35   ierr = CeedFree(&impl->evecsin); CeedChk(ierr);
36   ierr = CeedFree(&impl->qvecsin); CeedChk(ierr);
37 
38   for (CeedInt i=0; i<impl->numeout; i++) {
39     ierr = CeedVectorDestroy(&impl->evecsout[i]); CeedChk(ierr);
40     if (!impl->identityqf) {
41       ierr = CeedVectorDestroy(&impl->qvecsout[i]); CeedChk(ierr);
42     }
43   }
44   ierr = CeedFree(&impl->evecsout); CeedChk(ierr);
45   ierr = CeedFree(&impl->qvecsout); CeedChk(ierr);
46 
47   ierr = CeedFree(&impl); CeedChk(ierr);
48   return 0;
49 }
50 
51 /*
52   Setup infields or outfields
53  */
54 static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op,
55                                        bool inOrOut,
56                                        CeedVector *fullevecs, CeedVector *evecs,
57                                        CeedVector *qvecs, CeedInt starte,
58                                        CeedInt numfields, CeedInt Q) {
59   CeedInt dim, ierr, size, P;
60   Ceed ceed;
61   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
62   CeedBasis basis;
63   CeedElemRestriction Erestrict;
64   CeedOperatorField *opfields;
65   CeedQFunctionField *qffields;
66   if (inOrOut) {
67     ierr = CeedOperatorGetFields(op, NULL, &opfields); CeedChk(ierr);
68     ierr = CeedQFunctionGetFields(qf, NULL, &qffields); CeedChk(ierr);
69   } else {
70     ierr = CeedOperatorGetFields(op, &opfields, NULL); CeedChk(ierr);
71     ierr = CeedQFunctionGetFields(qf, &qffields, NULL); CeedChk(ierr);
72   }
73 
74   // Loop over fields
75   for (CeedInt i=0; i<numfields; i++) {
76     CeedEvalMode emode;
77     ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); CeedChk(ierr);
78 
79     if (emode != CEED_EVAL_WEIGHT) {
80       ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &Erestrict);
81       CeedChk(ierr);
82       ierr = CeedElemRestrictionCreateVector(Erestrict, NULL,
83                                              &fullevecs[i+starte]);
84       CeedChk(ierr);
85     }
86 
87     switch(emode) {
88     case CEED_EVAL_NONE:
89       ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChk(ierr);
90       ierr = CeedVectorCreate(ceed, Q*size, &qvecs[i]); CeedChk(ierr);
91       break;
92     case CEED_EVAL_INTERP:
93       ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChk(ierr);
94       ierr = CeedElemRestrictionGetElementSize(Erestrict, &P);
95       CeedChk(ierr);
96       ierr = CeedVectorCreate(ceed, P*size, &evecs[i]); CeedChk(ierr);
97       ierr = CeedVectorCreate(ceed, Q*size, &qvecs[i]); CeedChk(ierr);
98       break;
99     case CEED_EVAL_GRAD:
100       ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChk(ierr);
101       ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChk(ierr);
102       ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr);
103       ierr = CeedElemRestrictionGetElementSize(Erestrict, &P);
104       CeedChk(ierr);
105       ierr = CeedVectorCreate(ceed, P*size/dim, &evecs[i]); CeedChk(ierr);
106       ierr = CeedVectorCreate(ceed, Q*size, &qvecs[i]); CeedChk(ierr);
107       break;
108     case CEED_EVAL_WEIGHT: // Only on input fields
109       ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChk(ierr);
110       ierr = CeedVectorCreate(ceed, Q, &qvecs[i]); CeedChk(ierr);
111       ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT,
112                             NULL, qvecs[i]); CeedChk(ierr);
113       break;
114     case CEED_EVAL_DIV:
115       break; // Not implemented
116     case CEED_EVAL_CURL:
117       break; // Not implemented
118     }
119   }
120   return 0;
121 }
122 
123 /*
124   CeedOperator needs to connect all the named fields (be they active or passive)
125   to the named inputs and outputs of its CeedQFunction.
126  */
127 static int CeedOperatorSetup_Ref(CeedOperator op) {
128   int ierr;
129   bool setupdone;
130   ierr = CeedOperatorGetSetupStatus(op, &setupdone); CeedChk(ierr);
131   if (setupdone) return 0;
132   Ceed ceed;
133   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
134   CeedOperator_Ref *impl;
135   ierr = CeedOperatorGetData(op, (void *)&impl); CeedChk(ierr);
136   CeedQFunction qf;
137   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
138   CeedInt Q, numinputfields, numoutputfields;
139   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChk(ierr);
140   ierr = CeedQFunctionGetIdentityStatus(qf, &impl->identityqf); CeedChk(ierr);
141   ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
142   CeedChk(ierr);
143   CeedOperatorField *opinputfields, *opoutputfields;
144   ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields);
145   CeedChk(ierr);
146   CeedQFunctionField *qfinputfields, *qfoutputfields;
147   ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields);
148   CeedChk(ierr);
149 
150   // Allocate
151   ierr = CeedCalloc(numinputfields + numoutputfields, &impl->evecs);
152   CeedChk(ierr);
153   ierr = CeedCalloc(numinputfields + numoutputfields, &impl->edata);
154   CeedChk(ierr);
155 
156   ierr = CeedCalloc(16, &impl->inputstate); CeedChk(ierr);
157   ierr = CeedCalloc(16, &impl->evecsin); CeedChk(ierr);
158   ierr = CeedCalloc(16, &impl->evecsout); CeedChk(ierr);
159   ierr = CeedCalloc(16, &impl->qvecsin); CeedChk(ierr);
160   ierr = CeedCalloc(16, &impl->qvecsout); CeedChk(ierr);
161 
162   impl->numein = numinputfields; impl->numeout = numoutputfields;
163 
164   // Set up infield and outfield evecs and qvecs
165   // Infields
166   ierr = CeedOperatorSetupFields_Ref(qf, op, 0, impl->evecs,
167                                      impl->evecsin, impl->qvecsin, 0,
168                                      numinputfields, Q);
169   CeedChk(ierr);
170   // Outfields
171   ierr = CeedOperatorSetupFields_Ref(qf, op, 1, impl->evecs,
172                                      impl->evecsout, impl->qvecsout,
173                                      numinputfields, numoutputfields, Q);
174   CeedChk(ierr);
175 
176   // Identity QFunctions
177   if (impl->identityqf) {
178     CeedEvalMode inmode, outmode;
179     CeedQFunctionField *infields, *outfields;
180     ierr = CeedQFunctionGetFields(qf, &infields, &outfields); CeedChk(ierr);
181 
182     for (CeedInt i=0; i<numinputfields; i++) {
183       ierr = CeedQFunctionFieldGetEvalMode(infields[i], &inmode);
184       CeedChk(ierr);
185       ierr = CeedQFunctionFieldGetEvalMode(outfields[i], &outmode);
186       CeedChk(ierr);
187 
188       if (inmode == CEED_EVAL_NONE && outmode == CEED_EVAL_NONE)
189         // LCOV_EXCL_START
190         return CeedError(ceed, 1, "CEED_EVAL_NONE for a matching input and "
191                          "output does not make sense with identity QFunction");
192       // LCOV_EXCL_STOP
193 
194       ierr = CeedVectorDestroy(&impl->qvecsout[i]); CeedChk(ierr);
195       impl->qvecsout[i] = impl->qvecsin[i];
196     }
197   }
198 
199   ierr = CeedOperatorSetSetupDone(op); CeedChk(ierr);
200 
201   return 0;
202 }
203 
204 // Setup Input fields
205 static inline int CeedOperatorSetupInputs_Ref(CeedInt numinputfields,
206     CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields,
207     CeedVector invec, const bool skipactive, CeedOperator_Ref *impl,
208     CeedRequest *request) {
209   CeedInt ierr;
210   CeedEvalMode emode;
211   CeedVector vec;
212   CeedElemRestriction Erestrict;
213   CeedTransposeMode lmode;
214   uint64_t state;
215 
216   for (CeedInt i=0; i<numinputfields; i++) {
217     // Get input vector
218     ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
219     if (vec == CEED_VECTOR_ACTIVE) {
220       if (skipactive)
221         continue;
222       else
223         vec = invec;
224     }
225 
226     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
227     CeedChk(ierr);
228     // Restrict and Evec
229     if (emode == CEED_EVAL_WEIGHT) { // Skip
230     } else {
231       // Restrict
232       ierr = CeedVectorGetState(vec, &state); CeedChk(ierr);
233       // Skip restriction if input is unchanged
234       if (state != impl->inputstate[i] || vec == invec) {
235         ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict);
236         CeedChk(ierr);
237         ierr = CeedOperatorFieldGetLMode(opinputfields[i], &lmode);
238         CeedChk(ierr);
239         ierr = CeedElemRestrictionApply(Erestrict, CEED_NOTRANSPOSE,
240                                         lmode, vec, impl->evecs[i],
241                                         request); CeedChk(ierr);
242         impl->inputstate[i] = state;
243       }
244       // Get evec
245       ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST,
246                                     (const CeedScalar **) &impl->edata[i]);
247       CeedChk(ierr);
248     }
249   }
250   return 0;
251 }
252 
253 // Input basis action
254 static inline int CeedOperatorInputBasis_Ref(CeedInt e, CeedInt Q,
255     CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields,
256     CeedInt numinputfields, const bool skipactive, CeedOperator_Ref *impl) {
257   CeedInt ierr;
258   CeedInt dim, elemsize, size;
259   CeedElemRestriction Erestrict;
260   CeedEvalMode emode;
261   CeedBasis basis;
262 
263   for (CeedInt i=0; i<numinputfields; i++) {
264     // Skip active input
265     if (skipactive) {
266       CeedVector vec;
267       ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
268       if (vec == CEED_VECTOR_ACTIVE)
269         continue;
270     }
271     // Get elemsize, emode, size
272     ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict);
273     CeedChk(ierr);
274     ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize);
275     CeedChk(ierr);
276     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
277     CeedChk(ierr);
278     ierr = CeedQFunctionFieldGetSize(qfinputfields[i], &size); CeedChk(ierr);
279     // Basis action
280     switch(emode) {
281     case CEED_EVAL_NONE:
282       ierr = CeedVectorSetArray(impl->qvecsin[i], CEED_MEM_HOST,
283                                 CEED_USE_POINTER,
284                                 &impl->edata[i][e*Q*size]); CeedChk(ierr);
285       break;
286     case CEED_EVAL_INTERP:
287       ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); CeedChk(ierr);
288       ierr = CeedVectorSetArray(impl->evecsin[i], CEED_MEM_HOST,
289                                 CEED_USE_POINTER,
290                                 &impl->edata[i][e*elemsize*size]);
291       CeedChk(ierr);
292       ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE,
293                             CEED_EVAL_INTERP, impl->evecsin[i],
294                             impl->qvecsin[i]); CeedChk(ierr);
295       break;
296     case CEED_EVAL_GRAD:
297       ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); CeedChk(ierr);
298       ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr);
299       ierr = CeedVectorSetArray(impl->evecsin[i], CEED_MEM_HOST,
300                                 CEED_USE_POINTER,
301                                 &impl->edata[i][e*elemsize*size/dim]);
302       CeedChk(ierr);
303       ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE,
304                             CEED_EVAL_GRAD, impl->evecsin[i],
305                             impl->qvecsin[i]); CeedChk(ierr);
306       break;
307     case CEED_EVAL_WEIGHT:
308       break;  // No action
309     case CEED_EVAL_DIV:
310     case CEED_EVAL_CURL: {
311       // LCOV_EXCL_START
312       ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis);
313       CeedChk(ierr);
314       Ceed ceed;
315       ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr);
316       return CeedError(ceed, 1, "Ceed evaluation mode not implemented");
317       // LCOV_EXCL_STOP
318       break; // Not implemented
319     }
320     }
321   }
322   return 0;
323 }
324 
325 // Output basis action
326 static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q,
327     CeedQFunctionField *qfoutputfields, CeedOperatorField *opoutputfields,
328     CeedInt numinputfields, CeedInt numoutputfields, CeedOperator op,
329     CeedOperator_Ref *impl) {
330   CeedInt ierr;
331   CeedInt dim, elemsize, size;
332   CeedElemRestriction Erestrict;
333   CeedEvalMode emode;
334   CeedBasis basis;
335 
336   for (CeedInt i=0; i<numoutputfields; i++) {
337     // Get elemsize, emode, size
338     ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict);
339     CeedChk(ierr);
340     ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize);
341     CeedChk(ierr);
342     ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode);
343     CeedChk(ierr);
344     ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size); CeedChk(ierr);
345     // Basis action
346     switch(emode) {
347     case CEED_EVAL_NONE:
348       break; // No action
349     case CEED_EVAL_INTERP:
350       ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis);
351       CeedChk(ierr);
352       ierr = CeedVectorSetArray(impl->evecsout[i], CEED_MEM_HOST,
353                                 CEED_USE_POINTER,
354                                 &impl->edata[i + numinputfields][e*elemsize*size]);
355       CeedChk(ierr);
356       ierr = CeedBasisApply(basis, 1, CEED_TRANSPOSE,
357                             CEED_EVAL_INTERP, impl->qvecsout[i],
358                             impl->evecsout[i]); CeedChk(ierr);
359       break;
360     case CEED_EVAL_GRAD:
361       ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis);
362       CeedChk(ierr);
363       ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr);
364       ierr = CeedVectorSetArray(impl->evecsout[i], CEED_MEM_HOST,
365                                 CEED_USE_POINTER,
366                                 &impl->edata[i + numinputfields][e*elemsize*size/dim]);
367       CeedChk(ierr);
368       ierr = CeedBasisApply(basis, 1, CEED_TRANSPOSE,
369                             CEED_EVAL_GRAD, impl->qvecsout[i],
370                             impl->evecsout[i]); CeedChk(ierr);
371       break;
372     case CEED_EVAL_WEIGHT: {
373       // LCOV_EXCL_START
374       Ceed ceed;
375       ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
376       return CeedError(ceed, 1, "CEED_EVAL_WEIGHT cannot be an output "
377                        "evaluation mode");
378       // LCOV_EXCL_STOP
379       break; // Should not occur
380     }
381     case CEED_EVAL_DIV:
382     case CEED_EVAL_CURL: {
383       // LCOV_EXCL_START
384       Ceed ceed;
385       ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
386       return CeedError(ceed, 1, "Ceed evaluation mode not implemented");
387       // LCOV_EXCL_STOP
388       break; // Not implemented
389     }
390     }
391   }
392   return 0;
393 }
394 
395 // Restore Inputs
396 static inline int CeedOperatorRestoreInputs_Ref(CeedInt numinputfields,
397     CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields,
398     const bool skipactive, CeedOperator_Ref *impl) {
399   CeedInt ierr;
400   CeedEvalMode emode;
401 
402   for (CeedInt i=0; i<numinputfields; i++) {
403     // Skip active inputs
404     if (skipactive) {
405       CeedVector vec;
406       ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
407       if (vec == CEED_VECTOR_ACTIVE)
408         continue;
409     }
410     // Restore input
411     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
412     CeedChk(ierr);
413     if (emode == CEED_EVAL_WEIGHT) { // Skip
414     } else {
415       ierr = CeedVectorRestoreArrayRead(impl->evecs[i],
416                                         (const CeedScalar **) &impl->edata[i]);
417       CeedChk(ierr);
418     }
419   }
420   return 0;
421 }
422 
423 // Apply Operator
424 static int CeedOperatorApply_Ref(CeedOperator op, CeedVector invec,
425                                  CeedVector outvec, CeedRequest *request) {
426   int ierr;
427   CeedOperator_Ref *impl;
428   ierr = CeedOperatorGetData(op, (void *)&impl); CeedChk(ierr);
429   CeedQFunction qf;
430   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
431   CeedInt Q, numelements, numinputfields, numoutputfields, size;
432   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChk(ierr);
433   ierr = CeedOperatorGetNumElements(op, &numelements); CeedChk(ierr);
434   ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
435   CeedChk(ierr);
436   CeedTransposeMode lmode;
437   CeedOperatorField *opinputfields, *opoutputfields;
438   ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields);
439   CeedChk(ierr);
440   CeedQFunctionField *qfinputfields, *qfoutputfields;
441   ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields);
442   CeedChk(ierr);
443   CeedEvalMode emode;
444   CeedVector vec;
445   CeedElemRestriction Erestrict;
446 
447   // Setup
448   ierr = CeedOperatorSetup_Ref(op); CeedChk(ierr);
449 
450   // Input Evecs and Restriction
451   ierr = CeedOperatorSetupInputs_Ref(numinputfields, qfinputfields,
452                                      opinputfields, invec, false, impl,
453                                      request); CeedChk(ierr);
454 
455   // Output Evecs
456   for (CeedInt i=0; i<numoutputfields; i++) {
457     ierr = CeedVectorGetArray(impl->evecs[i+impl->numein], CEED_MEM_HOST,
458                               &impl->edata[i + numinputfields]); CeedChk(ierr);
459   }
460 
461   // Loop through elements
462   for (CeedInt e=0; e<numelements; e++) {
463     // Output pointers
464     for (CeedInt i=0; i<numoutputfields; i++) {
465       ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode);
466       CeedChk(ierr);
467       if (emode == CEED_EVAL_NONE) {
468         ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size);
469         CeedChk(ierr);
470         ierr = CeedVectorSetArray(impl->qvecsout[i], CEED_MEM_HOST,
471                                   CEED_USE_POINTER,
472                                   &impl->edata[i + numinputfields][e*Q*size]);
473         CeedChk(ierr);
474       }
475     }
476 
477     // Input basis apply
478     ierr = CeedOperatorInputBasis_Ref(e, Q, qfinputfields, opinputfields,
479                                       numinputfields, false, impl);
480     CeedChk(ierr);
481 
482     // Q function
483     if (!impl->identityqf) {
484       ierr = CeedQFunctionApply(qf, Q, impl->qvecsin, impl->qvecsout);
485       CeedChk(ierr);
486     }
487 
488     // Output basis apply
489     ierr = CeedOperatorOutputBasis_Ref(e, Q, qfoutputfields, opoutputfields,
490                                        numinputfields, numoutputfields, op, impl);
491     CeedChk(ierr);
492   }
493 
494   // Zero lvecs
495   for (CeedInt i=0; i<numoutputfields; i++) {
496     ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
497     if (vec == CEED_VECTOR_ACTIVE) {
498       if (!impl->add) {
499         vec = outvec;
500         ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr);
501       }
502     } else {
503       ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr);
504     }
505   }
506   impl->add = false;
507 
508   // Output restriction
509   for (CeedInt i=0; i<numoutputfields; i++) {
510     // Restore evec
511     ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein],
512                                   &impl->edata[i + numinputfields]);
513     CeedChk(ierr);
514     // Get output vector
515     ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
516     // Active
517     if (vec == CEED_VECTOR_ACTIVE)
518       vec = outvec;
519     // Restrict
520     ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict);
521     CeedChk(ierr);
522     ierr = CeedOperatorFieldGetLMode(opoutputfields[i], &lmode); CeedChk(ierr);
523     ierr = CeedElemRestrictionApply(Erestrict, CEED_TRANSPOSE,
524                                     lmode, impl->evecs[i+impl->numein], vec,
525                                     request); CeedChk(ierr);
526   }
527 
528   // Restore input arrays
529   ierr = CeedOperatorRestoreInputs_Ref(numinputfields, qfinputfields,
530                                        opinputfields, false, impl);
531   CeedChk(ierr);
532 
533   return 0;
534 }
535 
536 // Composite operators
537 static int CeedCompositeOperatorApply_Ref(CeedOperator op, CeedVector invec,
538     CeedVector outvec,
539     CeedRequest *request) {
540   int ierr;
541   CeedInt numsub;
542   CeedOperator_Ref *impl;
543   CeedOperator *suboperators;
544   ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr);
545   ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr);
546 
547   // Overwrite outvec with first output
548   ierr = CeedOperatorApply(suboperators[0], invec, outvec, request);
549   CeedChk(ierr);
550   // Add to outvec with subsequent outputs
551   for (CeedInt i=1; i<numsub; i++) {
552     ierr = CeedOperatorGetData(suboperators[i], (void *)&impl); CeedChk(ierr);
553     impl->add = true;
554     ierr = CeedOperatorApply(suboperators[i], invec, outvec, request);
555     CeedChk(ierr);
556   }
557 
558   return 0;
559 }
560 
561 // Assemble Linear QFunction
562 static int CeedOperatorAssembleLinearQFunction_Ref(CeedOperator op,
563     CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) {
564   int ierr;
565   CeedOperator_Ref *impl;
566   ierr = CeedOperatorGetData(op, (void *)&impl); CeedChk(ierr);
567   CeedQFunction qf;
568   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
569   CeedInt Q, numelements, numinputfields, numoutputfields, size;
570   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChk(ierr);
571   ierr = CeedOperatorGetNumElements(op, &numelements); CeedChk(ierr);
572   ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
573   CeedChk(ierr);
574   CeedOperatorField *opinputfields, *opoutputfields;
575   ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields);
576   CeedChk(ierr);
577   CeedQFunctionField *qfinputfields, *qfoutputfields;
578   ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields);
579   CeedChk(ierr);
580   CeedVector vec;
581   CeedInt numactivein = 0, numactiveout = 0;
582   CeedVector *activein = NULL;
583   CeedScalar *a, *tmp;
584   Ceed ceed;
585   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
586 
587   // Setup
588   ierr = CeedOperatorSetup_Ref(op); CeedChk(ierr);
589 
590   // Check for identity
591   if (impl->identityqf)
592    // LCOV_EXCL_START
593     return CeedError(ceed, 1, "Assembling identity qfunction does not make sense");
594   // LCOV_EXCL_STOP
595 
596   // Input Evecs and Restriction
597   ierr = CeedOperatorSetupInputs_Ref(numinputfields, qfinputfields,
598                                      opinputfields, NULL, true, impl, request);
599   CeedChk(ierr);
600 
601   // Count number of active input fields
602   for (CeedInt i=0; i<numinputfields; i++) {
603     // Get input vector
604     ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
605     // Check if active input
606     if (vec == CEED_VECTOR_ACTIVE) {
607       ierr = CeedQFunctionFieldGetSize(qfinputfields[i], &size); CeedChk(ierr);
608       ierr = CeedVectorSetValue(impl->qvecsin[i], 0.0); CeedChk(ierr);
609       ierr = CeedVectorGetArray(impl->qvecsin[i], CEED_MEM_HOST, &tmp);
610       CeedChk(ierr);
611       ierr = CeedRealloc(numactivein + size, &activein); CeedChk(ierr);
612       for (CeedInt field=0; field<size; field++) {
613         ierr = CeedVectorCreate(ceed, Q, &activein[numactivein+field]);
614         CeedChk(ierr);
615         ierr = CeedVectorSetArray(activein[numactivein+field], CEED_MEM_HOST,
616                                   CEED_USE_POINTER, &tmp[field*Q]);
617         CeedChk(ierr);
618       }
619       numactivein += size;
620       ierr = CeedVectorRestoreArray(impl->qvecsin[i], &tmp); CeedChk(ierr);
621     }
622   }
623 
624   // Count number of active output fields
625   for (CeedInt i=0; i<numoutputfields; i++) {
626     // Get output vector
627     ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
628     // Check if active output
629     if (vec == CEED_VECTOR_ACTIVE) {
630       ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size); CeedChk(ierr);
631       numactiveout += size;
632     }
633   }
634 
635   // Check sizes
636   if (!numactivein || !numactiveout)
637     // LCOV_EXCL_START
638     return CeedError(ceed, 1, "Cannot assemble QFunction without active inputs "
639                      "and outputs");
640   // LCOV_EXCL_STOP
641 
642   // Create output restriction
643   ierr = CeedElemRestrictionCreateIdentity(ceed, numelements, Q,
644          numelements*Q, numactivein*numactiveout, rstr); CeedChk(ierr);
645   // Create assembled vector
646   ierr = CeedVectorCreate(ceed, numelements*Q*numactivein*numactiveout,
647                           assembled); CeedChk(ierr);
648   ierr = CeedVectorSetValue(*assembled, 0.0); CeedChk(ierr);
649   ierr = CeedVectorGetArray(*assembled, CEED_MEM_HOST, &a); CeedChk(ierr);
650 
651   // Loop through elements
652   for (CeedInt e=0; e<numelements; e++) {
653     // Input basis apply
654     ierr = CeedOperatorInputBasis_Ref(e, Q, qfinputfields, opinputfields,
655                                       numinputfields, true, impl);
656     CeedChk(ierr);
657 
658     // Assemble QFunction
659     for (CeedInt in=0; in<numactivein; in++) {
660       // Set Inputs
661       ierr = CeedVectorSetValue(activein[in], 1.0); CeedChk(ierr);
662       if (numactivein > 1) {
663         ierr = CeedVectorSetValue(activein[(in+numactivein-1)%numactivein],
664                                   0.0); CeedChk(ierr);
665       }
666       // Set Outputs
667       for (CeedInt out=0; out<numoutputfields; out++) {
668         // Get output vector
669         ierr = CeedOperatorFieldGetVector(opoutputfields[out], &vec);
670         CeedChk(ierr);
671         // Check if active output
672         if (vec == CEED_VECTOR_ACTIVE) {
673           CeedVectorSetArray(impl->qvecsout[out], CEED_MEM_HOST,
674                              CEED_USE_POINTER, a); CeedChk(ierr);
675           ierr = CeedQFunctionFieldGetSize(qfoutputfields[out], &size);
676           CeedChk(ierr);
677           a += size*Q; // Advance the pointer by the size of the output
678         }
679       }
680       // Apply QFunction
681       ierr = CeedQFunctionApply(qf, Q, impl->qvecsin, impl->qvecsout);
682       CeedChk(ierr);
683     }
684   }
685 
686   // Un-set output Qvecs to prevent accidental overwrite of Assembled
687   for (CeedInt out=0; out<numoutputfields; out++) {
688     // Get output vector
689     ierr = CeedOperatorFieldGetVector(opoutputfields[out], &vec);
690     CeedChk(ierr);
691     // Check if active output
692     if (vec == CEED_VECTOR_ACTIVE) {
693       CeedVectorSetArray(impl->qvecsout[out], CEED_MEM_HOST, CEED_COPY_VALUES,
694                          NULL); CeedChk(ierr);
695     }
696   }
697 
698   // Restore input arrays
699   ierr = CeedOperatorRestoreInputs_Ref(numinputfields, qfinputfields,
700                                        opinputfields, true, impl);
701   CeedChk(ierr);
702 
703   // Restore output
704   ierr = CeedVectorRestoreArray(*assembled, &a); CeedChk(ierr);
705 
706   // Cleanup
707   for (CeedInt i=0; i<numactivein; i++) {
708     ierr = CeedVectorDestroy(&activein[i]); CeedChk(ierr);
709   }
710   ierr = CeedFree(&activein); CeedChk(ierr);
711 
712   return 0;
713 }
714 
715 int CeedOperatorCreate_Ref(CeedOperator op) {
716   int ierr;
717   Ceed ceed;
718   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
719   CeedOperator_Ref *impl;
720 
721   ierr = CeedCalloc(1, &impl); CeedChk(ierr);
722   impl->add = false;
723   ierr = CeedOperatorSetData(op, (void *)&impl); CeedChk(ierr);
724 
725   ierr = CeedSetBackendFunction(ceed, "Operator", op, "AssembleLinearQFunction",
726                                 CeedOperatorAssembleLinearQFunction_Ref);
727   CeedChk(ierr);
728   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Apply",
729                                 CeedOperatorApply_Ref); CeedChk(ierr);
730   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy",
731                                 CeedOperatorDestroy_Ref); CeedChk(ierr);
732   return 0;
733 }
734 
735 int CeedCompositeOperatorCreate_Ref(CeedOperator op) {
736   int ierr;
737   Ceed ceed;
738   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
739 
740   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Apply",
741                                 CeedCompositeOperatorApply_Ref); CeedChk(ierr);
742   return 0;
743 }
744