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