xref: /libCEED/backends/ref/ceed-ref-operator.c (revision 53d298d0b283eed7d6d59fec0c162c97e87aea98) 
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       ierr = CeedVectorDestroy(&impl->qvecsout[i]); CeedChk(ierr);
189       impl->qvecsout[i] = impl->qvecsin[i];
190     }
191   }
192 
193   ierr = CeedOperatorSetSetupDone(op); CeedChk(ierr);
194 
195   return 0;
196 }
197 
198 // Setup Input fields
199 static inline int CeedOperatorSetupInputs_Ref(CeedInt numinputfields,
200     CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields,
201     CeedVector invec, const bool skipactive, CeedOperator_Ref *impl,
202     CeedRequest *request) {
203   CeedInt ierr;
204   CeedEvalMode emode;
205   CeedVector vec;
206   CeedElemRestriction Erestrict;
207   CeedTransposeMode lmode;
208   uint64_t state;
209 
210   for (CeedInt i=0; i<numinputfields; i++) {
211     // Get input vector
212     ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
213     if (vec == CEED_VECTOR_ACTIVE) {
214       if (skipactive)
215         continue;
216       else
217         vec = invec;
218     }
219 
220     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
221     CeedChk(ierr);
222     // Restrict and Evec
223     if (emode == CEED_EVAL_WEIGHT) { // Skip
224     } else {
225       // Restrict
226       ierr = CeedVectorGetState(vec, &state); CeedChk(ierr);
227       // Skip restriction if input is unchanged
228       if (state != impl->inputstate[i] || vec == invec) {
229         ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict);
230         CeedChk(ierr);
231         ierr = CeedOperatorFieldGetLMode(opinputfields[i], &lmode);
232         CeedChk(ierr);
233         ierr = CeedElemRestrictionApply(Erestrict, CEED_NOTRANSPOSE,
234                                         lmode, vec, impl->evecs[i],
235                                         request); CeedChk(ierr);
236         impl->inputstate[i] = state;
237       }
238       // Get evec
239       ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST,
240                                     (const CeedScalar **) &impl->edata[i]);
241       CeedChk(ierr);
242     }
243   }
244   return 0;
245 }
246 
247 // Input basis action
248 static inline int CeedOperatorInputBasis_Ref(CeedInt e, CeedInt Q,
249     CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields,
250     CeedInt numinputfields, const bool skipactive, CeedOperator_Ref *impl) {
251   CeedInt ierr;
252   CeedInt dim, elemsize, size;
253   CeedElemRestriction Erestrict;
254   CeedEvalMode emode;
255   CeedBasis basis;
256 
257   for (CeedInt i=0; i<numinputfields; i++) {
258     // Skip active input
259     if (skipactive) {
260       CeedVector vec;
261       ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
262       if (vec == CEED_VECTOR_ACTIVE)
263         continue;
264     }
265     // Get elemsize, emode, size
266     ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict);
267     CeedChk(ierr);
268     ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize);
269     CeedChk(ierr);
270     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
271     CeedChk(ierr);
272     ierr = CeedQFunctionFieldGetSize(qfinputfields[i], &size); CeedChk(ierr);
273     // Basis action
274     switch(emode) {
275     case CEED_EVAL_NONE:
276       ierr = CeedVectorSetArray(impl->qvecsin[i], CEED_MEM_HOST,
277                                 CEED_USE_POINTER,
278                                 &impl->edata[i][e*Q*size]); CeedChk(ierr);
279       break;
280     case CEED_EVAL_INTERP:
281       ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); CeedChk(ierr);
282       ierr = CeedVectorSetArray(impl->evecsin[i], CEED_MEM_HOST,
283                                 CEED_USE_POINTER,
284                                 &impl->edata[i][e*elemsize*size]);
285       CeedChk(ierr);
286       ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE,
287                             CEED_EVAL_INTERP, impl->evecsin[i],
288                             impl->qvecsin[i]); CeedChk(ierr);
289       break;
290     case CEED_EVAL_GRAD:
291       ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); CeedChk(ierr);
292       ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr);
293       ierr = CeedVectorSetArray(impl->evecsin[i], CEED_MEM_HOST,
294                                 CEED_USE_POINTER,
295                                 &impl->edata[i][e*elemsize*size/dim]);
296       CeedChk(ierr);
297       ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE,
298                             CEED_EVAL_GRAD, impl->evecsin[i],
299                             impl->qvecsin[i]); CeedChk(ierr);
300       break;
301     case CEED_EVAL_WEIGHT:
302       break;  // No action
303     case CEED_EVAL_DIV:
304     case CEED_EVAL_CURL: {
305       // LCOV_EXCL_START
306       ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis);
307       CeedChk(ierr);
308       Ceed ceed;
309       ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr);
310       return CeedError(ceed, 1, "Ceed evaluation mode not implemented");
311       // LCOV_EXCL_STOP
312       break; // Not implemented
313     }
314     }
315   }
316   return 0;
317 }
318 
319 // Output basis action
320 static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q,
321     CeedQFunctionField *qfoutputfields, CeedOperatorField *opoutputfields,
322     CeedInt numinputfields, CeedInt numoutputfields, CeedOperator op,
323     CeedOperator_Ref *impl) {
324   CeedInt ierr;
325   CeedInt dim, elemsize, size;
326   CeedElemRestriction Erestrict;
327   CeedEvalMode emode;
328   CeedBasis basis;
329 
330   for (CeedInt i=0; i<numoutputfields; i++) {
331     // Get elemsize, emode, size
332     ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict);
333     CeedChk(ierr);
334     ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize);
335     CeedChk(ierr);
336     ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode);
337     CeedChk(ierr);
338     ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size); CeedChk(ierr);
339     // Basis action
340     switch(emode) {
341     case CEED_EVAL_NONE:
342       break; // No action
343     case CEED_EVAL_INTERP:
344       ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis);
345       CeedChk(ierr);
346       ierr = CeedVectorSetArray(impl->evecsout[i], CEED_MEM_HOST,
347                                 CEED_USE_POINTER,
348                                 &impl->edata[i + numinputfields][e*elemsize*size]);
349       CeedChk(ierr);
350       ierr = CeedBasisApply(basis, 1, CEED_TRANSPOSE,
351                             CEED_EVAL_INTERP, impl->qvecsout[i],
352                             impl->evecsout[i]); CeedChk(ierr);
353       break;
354     case CEED_EVAL_GRAD:
355       ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis);
356       CeedChk(ierr);
357       ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr);
358       ierr = CeedVectorSetArray(impl->evecsout[i], CEED_MEM_HOST,
359                                 CEED_USE_POINTER,
360                                 &impl->edata[i + numinputfields][e*elemsize*size/dim]);
361       CeedChk(ierr);
362       ierr = CeedBasisApply(basis, 1, CEED_TRANSPOSE,
363                             CEED_EVAL_GRAD, impl->qvecsout[i],
364                             impl->evecsout[i]); CeedChk(ierr);
365       break;
366     case CEED_EVAL_WEIGHT: {
367       // LCOV_EXCL_START
368       Ceed ceed;
369       ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
370       return CeedError(ceed, 1, "CEED_EVAL_WEIGHT cannot be an output "
371                        "evaluation mode");
372       // LCOV_EXCL_STOP
373       break; // Should not occur
374     }
375     case CEED_EVAL_DIV:
376     case CEED_EVAL_CURL: {
377       // LCOV_EXCL_START
378       Ceed ceed;
379       ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
380       return CeedError(ceed, 1, "Ceed evaluation mode not implemented");
381       // LCOV_EXCL_STOP
382       break; // Not implemented
383     }
384     }
385   }
386   return 0;
387 }
388 
389 // Restore Inputs
390 static inline int CeedOperatorRestoreInputs_Ref(CeedInt numinputfields,
391     CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields,
392     const bool skipactive, CeedOperator_Ref *impl) {
393   CeedInt ierr;
394   CeedEvalMode emode;
395 
396   for (CeedInt i=0; i<numinputfields; i++) {
397     // Skip active inputs
398     if (skipactive) {
399       CeedVector vec;
400       ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
401       if (vec == CEED_VECTOR_ACTIVE)
402         continue;
403     }
404     // Restore input
405     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
406     CeedChk(ierr);
407     if (emode == CEED_EVAL_WEIGHT) { // Skip
408     } else {
409       ierr = CeedVectorRestoreArrayRead(impl->evecs[i],
410                                         (const CeedScalar **) &impl->edata[i]);
411       CeedChk(ierr);
412     }
413   }
414   return 0;
415 }
416 
417 // Apply Operator
418 static int CeedOperatorApply_Ref(CeedOperator op, CeedVector invec,
419                                  CeedVector outvec, CeedRequest *request) {
420   int ierr;
421   CeedOperator_Ref *impl;
422   ierr = CeedOperatorGetData(op, (void *)&impl); CeedChk(ierr);
423   CeedQFunction qf;
424   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
425   CeedInt Q, numelements, numinputfields, numoutputfields, size;
426   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChk(ierr);
427   ierr = CeedOperatorGetNumElements(op, &numelements); CeedChk(ierr);
428   ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
429   CeedChk(ierr);
430   CeedTransposeMode lmode;
431   CeedOperatorField *opinputfields, *opoutputfields;
432   ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields);
433   CeedChk(ierr);
434   CeedQFunctionField *qfinputfields, *qfoutputfields;
435   ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields);
436   CeedChk(ierr);
437   CeedEvalMode emode;
438   CeedVector vec;
439   CeedElemRestriction Erestrict;
440 
441   // Setup
442   ierr = CeedOperatorSetup_Ref(op); CeedChk(ierr);
443 
444   // Input Evecs and Restriction
445   ierr = CeedOperatorSetupInputs_Ref(numinputfields, qfinputfields,
446                                      opinputfields, invec, false, impl,
447                                      request); CeedChk(ierr);
448 
449   // Output Evecs
450   for (CeedInt i=0; i<numoutputfields; i++) {
451     ierr = CeedVectorGetArray(impl->evecs[i+impl->numein], CEED_MEM_HOST,
452                               &impl->edata[i + numinputfields]); CeedChk(ierr);
453   }
454 
455   // Loop through elements
456   for (CeedInt e=0; e<numelements; e++) {
457     // Output pointers
458     for (CeedInt i=0; i<numoutputfields; i++) {
459       ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode);
460       CeedChk(ierr);
461       if (emode == CEED_EVAL_NONE) {
462         ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size);
463         CeedChk(ierr);
464         ierr = CeedVectorSetArray(impl->qvecsout[i], CEED_MEM_HOST,
465                                   CEED_USE_POINTER,
466                                   &impl->edata[i + numinputfields][e*Q*size]);
467         CeedChk(ierr);
468       }
469     }
470 
471     // Input basis apply
472     ierr = CeedOperatorInputBasis_Ref(e, Q, qfinputfields, opinputfields,
473                                       numinputfields, false, impl);
474     CeedChk(ierr);
475 
476     // Q function
477     if (!impl->identityqf) {
478       ierr = CeedQFunctionApply(qf, Q, impl->qvecsin, impl->qvecsout);
479       CeedChk(ierr);
480     }
481 
482     // Output basis apply
483     ierr = CeedOperatorOutputBasis_Ref(e, Q, qfoutputfields, opoutputfields,
484                                        numinputfields, numoutputfields, op, impl);
485     CeedChk(ierr);
486   }
487 
488   // Zero lvecs
489   for (CeedInt i=0; i<numoutputfields; i++) {
490     ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
491     if (vec == CEED_VECTOR_ACTIVE) {
492       if (!impl->add) {
493         vec = outvec;
494         ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr);
495       }
496     } else {
497       ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr);
498     }
499   }
500   impl->add = false;
501 
502   // Output restriction
503   for (CeedInt i=0; i<numoutputfields; i++) {
504     // Restore evec
505     ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein],
506                                   &impl->edata[i + numinputfields]);
507     CeedChk(ierr);
508     // Get output vector
509     ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
510     // Active
511     if (vec == CEED_VECTOR_ACTIVE)
512       vec = outvec;
513     // Restrict
514     ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict);
515     CeedChk(ierr);
516     ierr = CeedOperatorFieldGetLMode(opoutputfields[i], &lmode); CeedChk(ierr);
517     ierr = CeedElemRestrictionApply(Erestrict, CEED_TRANSPOSE,
518                                     lmode, impl->evecs[i+impl->numein], vec,
519                                     request); CeedChk(ierr);
520   }
521 
522   // Restore input arrays
523   ierr = CeedOperatorRestoreInputs_Ref(numinputfields, qfinputfields,
524                                        opinputfields, false, impl);
525   CeedChk(ierr);
526 
527   return 0;
528 }
529 
530 // Composite operators
531 static int CeedCompositeOperatorApply_Ref(CeedOperator op, CeedVector invec,
532     CeedVector outvec,
533     CeedRequest *request) {
534   int ierr;
535   CeedInt numsub;
536   CeedOperator_Ref *impl;
537   CeedOperator *suboperators;
538   ierr = CeedOperatorGetNumSub(op, &numsub); CeedChk(ierr);
539   ierr = CeedOperatorGetSubList(op, &suboperators); CeedChk(ierr);
540 
541   // Overwrite outvec with first output
542   ierr = CeedOperatorApply(suboperators[0], invec, outvec, request);
543   CeedChk(ierr);
544   // Add to outvec with subsequent outputs
545   for (CeedInt i=1; i<numsub; i++) {
546     ierr = CeedOperatorGetData(suboperators[i], (void *)&impl); CeedChk(ierr);
547     impl->add = true;
548     ierr = CeedOperatorApply(suboperators[i], invec, outvec, request);
549     CeedChk(ierr);
550   }
551 
552   return 0;
553 }
554 
555 // Assemble Linear QFunction
556 static int CeedOperatorAssembleLinearQFunction_Ref(CeedOperator op,
557     CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) {
558   int ierr;
559   CeedOperator_Ref *impl;
560   ierr = CeedOperatorGetData(op, (void *)&impl); CeedChk(ierr);
561   CeedQFunction qf;
562   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
563   CeedInt Q, numelements, numinputfields, numoutputfields, size;
564   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChk(ierr);
565   ierr = CeedOperatorGetNumElements(op, &numelements); CeedChk(ierr);
566   ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
567   CeedChk(ierr);
568   CeedOperatorField *opinputfields, *opoutputfields;
569   ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields);
570   CeedChk(ierr);
571   CeedQFunctionField *qfinputfields, *qfoutputfields;
572   ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields);
573   CeedChk(ierr);
574   CeedVector vec;
575   CeedInt numactivein = 0, numactiveout = 0;
576   CeedVector *activein = NULL;
577   CeedScalar *a, *tmp;
578   Ceed ceed;
579   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
580 
581   // Setup
582   ierr = CeedOperatorSetup_Ref(op); CeedChk(ierr);
583 
584   // Check for identity
585   if (impl->identityqf)
586     // LCOV_EXCL_START
587     return CeedError(ceed, 1, "Assembling identity qfunction does not make sense");
588   // LCOV_EXCL_STOP
589 
590   // Input Evecs and Restriction
591   ierr = CeedOperatorSetupInputs_Ref(numinputfields, qfinputfields,
592                                      opinputfields, NULL, true, impl, request);
593   CeedChk(ierr);
594 
595   // Count number of active input fields
596   for (CeedInt i=0; i<numinputfields; i++) {
597     // Get input vector
598     ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
599     // Check if active input
600     if (vec == CEED_VECTOR_ACTIVE) {
601       ierr = CeedQFunctionFieldGetSize(qfinputfields[i], &size); CeedChk(ierr);
602       ierr = CeedVectorSetValue(impl->qvecsin[i], 0.0); CeedChk(ierr);
603       ierr = CeedVectorGetArray(impl->qvecsin[i], CEED_MEM_HOST, &tmp);
604       CeedChk(ierr);
605       ierr = CeedRealloc(numactivein + size, &activein); CeedChk(ierr);
606       for (CeedInt field=0; field<size; field++) {
607         ierr = CeedVectorCreate(ceed, Q, &activein[numactivein+field]);
608         CeedChk(ierr);
609         ierr = CeedVectorSetArray(activein[numactivein+field], CEED_MEM_HOST,
610                                   CEED_USE_POINTER, &tmp[field*Q]);
611         CeedChk(ierr);
612       }
613       numactivein += size;
614       ierr = CeedVectorRestoreArray(impl->qvecsin[i], &tmp); CeedChk(ierr);
615     }
616   }
617 
618   // Count number of active output fields
619   for (CeedInt i=0; i<numoutputfields; i++) {
620     // Get output vector
621     ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
622     // Check if active output
623     if (vec == CEED_VECTOR_ACTIVE) {
624       ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size); CeedChk(ierr);
625       numactiveout += size;
626     }
627   }
628 
629   // Check sizes
630   if (!numactivein || !numactiveout)
631     // LCOV_EXCL_START
632     return CeedError(ceed, 1, "Cannot assemble QFunction without active inputs "
633                      "and outputs");
634   // LCOV_EXCL_STOP
635 
636   // Create output restriction
637   ierr = CeedElemRestrictionCreateIdentity(ceed, numelements, Q,
638          numelements*Q, numactivein*numactiveout, rstr); CeedChk(ierr);
639   // Create assembled vector
640   ierr = CeedVectorCreate(ceed, numelements*Q*numactivein*numactiveout,
641                           assembled); CeedChk(ierr);
642   ierr = CeedVectorSetValue(*assembled, 0.0); CeedChk(ierr);
643   ierr = CeedVectorGetArray(*assembled, CEED_MEM_HOST, &a); CeedChk(ierr);
644 
645   // Loop through elements
646   for (CeedInt e=0; e<numelements; e++) {
647     // Input basis apply
648     ierr = CeedOperatorInputBasis_Ref(e, Q, qfinputfields, opinputfields,
649                                       numinputfields, true, impl);
650     CeedChk(ierr);
651 
652     // Assemble QFunction
653     for (CeedInt in=0; in<numactivein; in++) {
654       // Set Inputs
655       ierr = CeedVectorSetValue(activein[in], 1.0); CeedChk(ierr);
656       if (numactivein > 1) {
657         ierr = CeedVectorSetValue(activein[(in+numactivein-1)%numactivein],
658                                   0.0); CeedChk(ierr);
659       }
660       // Set Outputs
661       for (CeedInt out=0; out<numoutputfields; out++) {
662         // Get output vector
663         ierr = CeedOperatorFieldGetVector(opoutputfields[out], &vec);
664         CeedChk(ierr);
665         // Check if active output
666         if (vec == CEED_VECTOR_ACTIVE) {
667           CeedVectorSetArray(impl->qvecsout[out], CEED_MEM_HOST,
668                              CEED_USE_POINTER, a); CeedChk(ierr);
669           ierr = CeedQFunctionFieldGetSize(qfoutputfields[out], &size);
670           CeedChk(ierr);
671           a += size*Q; // Advance the pointer by the size of the output
672         }
673       }
674       // Apply QFunction
675       ierr = CeedQFunctionApply(qf, Q, impl->qvecsin, impl->qvecsout);
676       CeedChk(ierr);
677     }
678   }
679 
680   // Un-set output Qvecs to prevent accidental overwrite of Assembled
681   for (CeedInt out=0; out<numoutputfields; out++) {
682     // Get output vector
683     ierr = CeedOperatorFieldGetVector(opoutputfields[out], &vec);
684     CeedChk(ierr);
685     // Check if active output
686     if (vec == CEED_VECTOR_ACTIVE) {
687       CeedVectorSetArray(impl->qvecsout[out], CEED_MEM_HOST, CEED_COPY_VALUES,
688                          NULL); CeedChk(ierr);
689     }
690   }
691 
692   // Restore input arrays
693   ierr = CeedOperatorRestoreInputs_Ref(numinputfields, qfinputfields,
694                                        opinputfields, true, impl);
695   CeedChk(ierr);
696 
697   // Restore output
698   ierr = CeedVectorRestoreArray(*assembled, &a); CeedChk(ierr);
699 
700   // Cleanup
701   for (CeedInt i=0; i<numactivein; i++) {
702     ierr = CeedVectorDestroy(&activein[i]); CeedChk(ierr);
703   }
704   ierr = CeedFree(&activein); CeedChk(ierr);
705 
706   return 0;
707 }
708 
709 int CeedOperatorCreate_Ref(CeedOperator op) {
710   int ierr;
711   Ceed ceed;
712   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
713   CeedOperator_Ref *impl;
714 
715   ierr = CeedCalloc(1, &impl); CeedChk(ierr);
716   impl->add = false;
717   ierr = CeedOperatorSetData(op, (void *)&impl); CeedChk(ierr);
718 
719   ierr = CeedSetBackendFunction(ceed, "Operator", op, "AssembleLinearQFunction",
720                                 CeedOperatorAssembleLinearQFunction_Ref);
721   CeedChk(ierr);
722   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Apply",
723                                 CeedOperatorApply_Ref); CeedChk(ierr);
724   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy",
725                                 CeedOperatorDestroy_Ref); CeedChk(ierr);
726   return 0;
727 }
728 
729 int CeedCompositeOperatorCreate_Ref(CeedOperator op) {
730   int ierr;
731   Ceed ceed;
732   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
733 
734   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Apply",
735                                 CeedCompositeOperatorApply_Ref); CeedChk(ierr);
736   return 0;
737 }
738