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