xref: /libCEED/rust/libceed-sys/c-src/interface/ceed-fortran.c (revision 12bf2c1193fc7a0678cab875df72cd7865cfc60b)
1 // Fortran interface
2 #include <ceed.h>
3 #include <ceed-impl.h>
4 #include <ceed-fortran-name.h>
5 
6 #include <stdlib.h>
7 #include <string.h>
8 
9 #define FORTRAN_REQUEST_IMMEDIATE -1
10 #define FORTRAN_REQUEST_ORDERED -2
11 #define FORTRAN_NULL -3
12 #define FORTRAN_BASIS_COLOCATED -1
13 #define FORTRAN_QDATA_NONE -1
14 #define FORTRAN_VECTOR_ACTIVE -1
15 #define FORTRAN_VECTOR_NONE -2
16 
17 static Ceed *Ceed_dict = NULL;
18 static int Ceed_count = 0;
19 static int Ceed_n = 0;
20 static int Ceed_count_max = 0;
21 
22 #define fCeedInit FORTRAN_NAME(ceedinit,CEEDINIT)
23 void fCeedInit(const char* resource, int *ceed, int *err) {
24   if (Ceed_count == Ceed_count_max) {
25     Ceed_count_max += Ceed_count_max/2 + 1;
26     CeedRealloc(Ceed_count_max, &Ceed_dict);
27   }
28 
29   Ceed *ceed_ = &Ceed_dict[Ceed_count];
30   *err = CeedInit(resource, ceed_);
31 
32   if (*err == 0) {
33     *ceed = Ceed_count++;
34     Ceed_n++;
35   }
36 }
37 
38 #define fCeedDestroy FORTRAN_NAME(ceeddestroy,CEEDDESTROY)
39 void fCeedDestroy(int *ceed, int *err) {
40   *err = CeedDestroy(&Ceed_dict[*ceed]);
41 
42   if (*err == 0) {
43     Ceed_n--;
44     if (Ceed_n == 0) {
45       CeedFree(&Ceed_dict);
46       Ceed_count = 0;
47       Ceed_count_max = 0;
48     }
49   }
50 }
51 
52 static CeedVector *CeedVector_dict = NULL;
53 static int CeedVector_count = 0;
54 static int CeedVector_n = 0;
55 static int CeedVector_count_max = 0;
56 
57 #define fCeedVectorCreate FORTRAN_NAME(ceedvectorcreate,CEEDVECTORCREATE)
58 void fCeedVectorCreate(int *ceed, int *length, int *vec, int *err) {
59   if (CeedVector_count == CeedVector_count_max) {
60     CeedVector_count_max += CeedVector_count_max/2 + 1;
61     CeedRealloc(CeedVector_count_max, &CeedVector_dict);
62   }
63 
64   CeedVector* vec_ = &CeedVector_dict[CeedVector_count];
65   *err = CeedVectorCreate(Ceed_dict[*ceed], *length, vec_);
66 
67   if (*err == 0) {
68     *vec = CeedVector_count++;
69     CeedVector_n++;
70   }
71 }
72 
73 #define fCeedVectorSetArray FORTRAN_NAME(ceedvectorsetarray,CEEDVECTORSETARRAY)
74 void fCeedVectorSetArray(int *vec, int *memtype, int *copymode,
75                          CeedScalar *array, int *err) {
76   *err = CeedVectorSetArray(CeedVector_dict[*vec], *memtype, *copymode, array);
77 }
78 
79 #define fCeedVectorSetValue FORTRAN_NAME(ceedvectorsetvalue,CEEDVECTORSETVALUE)
80 void fCeedVectorSetValue(int *vec, CeedScalar *value, int *err) {
81   *err = CeedVectorSetValue(CeedVector_dict[*vec], *value);
82 }
83 
84 #define fCeedVectorGetArray FORTRAN_NAME(ceedvectorgetarray,CEEDVECTORGETARRAY)
85 void fCeedVectorGetArray(int *vec, int *memtype, CeedScalar *array, int64_t *offset,
86                          int *err) {
87   CeedScalar *b;
88   CeedVector vec_ = CeedVector_dict[*vec];
89   *err = CeedVectorGetArray(vec_, *memtype, &b);
90   *offset = b - array;
91 }
92 
93 #define fCeedVectorGetArrayRead \
94     FORTRAN_NAME(ceedvectorgetarrayread,CEEDVECTORGETARRAYREAD)
95 void fCeedVectorGetArrayRead(int *vec, int *memtype, CeedScalar *array,
96                              int64_t *offset, int *err) {
97   const CeedScalar *b;
98   CeedVector vec_ = CeedVector_dict[*vec];
99   *err = CeedVectorGetArrayRead(vec_, *memtype, &b);
100   *offset = b - array;
101 }
102 
103 #define fCeedVectorRestoreArray \
104     FORTRAN_NAME(ceedvectorrestorearray,CEEDVECTORRESTOREARRAY)
105 void fCeedVectorRestoreArray(int *vec, CeedScalar *array,
106                              int64_t *offset, int *err) {
107   *err = CeedVectorRestoreArray(CeedVector_dict[*vec], &array);
108   *offset = 0;
109 }
110 
111 #define fCeedVectorRestoreArrayRead \
112     FORTRAN_NAME(ceedvectorrestorearrayread,CEEDVECTORRESTOREARRAYREAD)
113 void fCeedVectorRestoreArrayRead(int *vec, const CeedScalar *array,
114                                  int64_t *offset, int *err) {
115   *err = CeedVectorRestoreArrayRead(CeedVector_dict[*vec], &array);
116   *offset = 0;
117 }
118 
119 #define fCeedVectorView FORTRAN_NAME(ceedvectorview,CEEDVECTORVIEW)
120 void fCeedVectorView(int *vec, int *err) {
121   *err = CeedVectorView(CeedVector_dict[*vec], "%12.8f", stdout);
122 }
123 
124 #define fCeedVectorDestroy FORTRAN_NAME(ceedvectordestroy,CEEDVECTORDESTROY)
125 void fCeedVectorDestroy(int *vec, int *err) {
126   *err = CeedVectorDestroy(&CeedVector_dict[*vec]);
127 
128   if (*err == 0) {
129     CeedVector_n--;
130     if (CeedVector_n == 0) {
131       CeedFree(&CeedVector_dict);
132       CeedVector_count = 0;
133       CeedVector_count_max = 0;
134     }
135   }
136 }
137 
138 static CeedElemRestriction *CeedElemRestriction_dict = NULL;
139 static int CeedElemRestriction_count = 0;
140 static int CeedElemRestriction_n = 0;
141 static int CeedElemRestriction_count_max = 0;
142 
143 #define fCeedElemRestrictionCreate \
144     FORTRAN_NAME(ceedelemrestrictioncreate, CEEDELEMRESTRICTIONCREATE)
145 void fCeedElemRestrictionCreate(int *ceed, int *nelements,
146                                 int *esize, int *ndof, int *ncomp, int *memtype, int *copymode,
147                                 const int *indices, int *elemrestriction, int *err) {
148   if (CeedElemRestriction_count == CeedElemRestriction_count_max) {
149     CeedElemRestriction_count_max += CeedElemRestriction_count_max/2 + 1;
150     CeedRealloc(CeedElemRestriction_count_max, &CeedElemRestriction_dict);
151   }
152 
153   const int *indices_ = indices;
154 
155   CeedElemRestriction *elemrestriction_ =
156     &CeedElemRestriction_dict[CeedElemRestriction_count];
157   *err = CeedElemRestrictionCreate(Ceed_dict[*ceed], *nelements, *esize, *ndof,
158                                    *ncomp,
159                                    *memtype, *copymode, indices_, elemrestriction_);
160 
161   if (*err == 0) {
162     *elemrestriction = CeedElemRestriction_count++;
163     CeedElemRestriction_n++;
164   }
165 }
166 
167 #define fCeedElemRestrictionCreateIdentity \
168     FORTRAN_NAME(ceedelemrestrictioncreateidentity, CEEDELEMRESTRICTIONCREATEIDENTITY)
169 void fCeedElemRestrictionCreateIdentity(int *ceed, int *nelements,
170                                         int *esize, int *ndof, int *ncomp,
171                                         int *elemrestriction, int *err) {
172   if (CeedElemRestriction_count == CeedElemRestriction_count_max) {
173     CeedElemRestriction_count_max += CeedElemRestriction_count_max/2 + 1;
174     CeedRealloc(CeedElemRestriction_count_max, &CeedElemRestriction_dict);
175   }
176 
177   CeedElemRestriction *elemrestriction_ =
178     &CeedElemRestriction_dict[CeedElemRestriction_count];
179   *err = CeedElemRestrictionCreateIdentity(Ceed_dict[*ceed], *nelements, *esize,
180          *ndof,
181          *ncomp, elemrestriction_);
182 
183   if (*err == 0) {
184     *elemrestriction = CeedElemRestriction_count++;
185     CeedElemRestriction_n++;
186   }
187 }
188 
189 #define fCeedElemRestrictionCreateBlocked \
190     FORTRAN_NAME(ceedelemrestrictioncreateblocked,CEEDELEMRESTRICTIONCREATEBLOCKED)
191 void fCeedElemRestrictionCreateBlocked(int *ceed, int *nelements,
192                                        int *esize, int *blocksize, int *ndof, int *ncomp,
193                                        int *mtype, int *cmode,
194                                        int *blkindices, int *elemrestriction, int *err) {
195 
196   if (CeedElemRestriction_count == CeedElemRestriction_count_max) {
197     CeedElemRestriction_count_max += CeedElemRestriction_count_max/2 + 1;
198     CeedRealloc(CeedElemRestriction_count_max, &CeedElemRestriction_dict);
199   }
200 
201   CeedElemRestriction *elemrestriction_ =
202     &CeedElemRestriction_dict[CeedElemRestriction_count];
203   *err = CeedElemRestrictionCreateBlocked(Ceed_dict[*ceed], *nelements, *esize,
204                                           *blocksize, *ndof, *ncomp, *mtype, *cmode, blkindices,
205                                           elemrestriction_);
206 
207   if (*err == 0) {
208     *elemrestriction = CeedElemRestriction_count++;
209     CeedElemRestriction_n++;
210   }
211 }
212 
213 static CeedRequest *CeedRequest_dict = NULL;
214 static int CeedRequest_count = 0;
215 static int CeedRequest_n = 0;
216 static int CeedRequest_count_max = 0;
217 
218 #define fCeedElemRestrictionApply \
219     FORTRAN_NAME(ceedelemrestrictionapply,CEEDELEMRESTRICTIONAPPLY)
220 void fCeedElemRestrictionApply(int *elemr, int *tmode, int *lmode,
221                                int *uvec, int *ruvec, int *rqst, int *err) {
222   int createRequest = 1;
223   // Check if input is CEED_REQUEST_ORDERED(-2) or CEED_REQUEST_IMMEDIATE(-1)
224   if (*rqst == FORTRAN_REQUEST_IMMEDIATE || *rqst == FORTRAN_REQUEST_ORDERED)
225     createRequest = 0;
226 
227   if (createRequest && CeedRequest_count == CeedRequest_count_max) {
228     CeedRequest_count_max += CeedRequest_count_max/2 + 1;
229     CeedRealloc(CeedRequest_count_max, &CeedRequest_dict);
230   }
231 
232   CeedRequest *rqst_;
233   if      (*rqst == FORTRAN_REQUEST_IMMEDIATE) rqst_ = CEED_REQUEST_IMMEDIATE;
234   else if (*rqst == FORTRAN_REQUEST_ORDERED  ) rqst_ = CEED_REQUEST_ORDERED;
235   else rqst_ = &CeedRequest_dict[CeedRequest_count];
236 
237   *err = CeedElemRestrictionApply(CeedElemRestriction_dict[*elemr], *tmode,
238                                   *lmode, CeedVector_dict[*uvec], CeedVector_dict[*ruvec], rqst_);
239 
240   if (*err == 0 && createRequest) {
241     *rqst = CeedRequest_count++;
242     CeedRequest_n++;
243   }
244 }
245 
246 #define fCeedRequestWait FORTRAN_NAME(ceedrequestwait, CEEDREQUESTWAIT)
247 void fCeedRequestWait(int *rqst, int *err) {
248   // TODO Uncomment this once CeedRequestWait is implemented
249   //*err = CeedRequestWait(&CeedRequest_dict[*rqst]);
250 
251   if (*err == 0) {
252     CeedRequest_n--;
253     if (CeedRequest_n == 0) {
254       CeedFree(&CeedRequest_dict);
255       CeedRequest_count = 0;
256       CeedRequest_count_max = 0;
257     }
258   }
259 }
260 
261 #define fCeedElemRestrictionDestroy \
262     FORTRAN_NAME(ceedelemrestrictiondestroy,CEEDELEMRESTRICTIONDESTROY)
263 void fCeedElemRestrictionDestroy(int *elem, int *err) {
264   *err = CeedElemRestrictionDestroy(&CeedElemRestriction_dict[*elem]);
265 
266   if (*err == 0) {
267     CeedElemRestriction_n--;
268     if (CeedElemRestriction_n == 0) {
269       CeedFree(&CeedElemRestriction_dict);
270       CeedElemRestriction_count = 0;
271       CeedElemRestriction_count_max = 0;
272     }
273   }
274 }
275 
276 static CeedBasis *CeedBasis_dict = NULL;
277 static int CeedBasis_count = 0;
278 static int CeedBasis_n = 0;
279 static int CeedBasis_count_max = 0;
280 
281 #define fCeedBasisCreateTensorH1Lagrange \
282     FORTRAN_NAME(ceedbasiscreatetensorh1lagrange, CEEDBASISCREATETENSORH1LAGRANGE)
283 void fCeedBasisCreateTensorH1Lagrange(int *ceed, int *dim,
284                                       int *ndof, int *P, int *Q, int *quadmode, int *basis,
285                                       int *err) {
286   if (CeedBasis_count == CeedBasis_count_max) {
287     CeedBasis_count_max += CeedBasis_count_max/2 + 1;
288     CeedRealloc(CeedBasis_count_max, &CeedBasis_dict);
289   }
290 
291   *err = CeedBasisCreateTensorH1Lagrange(Ceed_dict[*ceed], *dim, *ndof, *P, *Q,
292                                          *quadmode, &CeedBasis_dict[CeedBasis_count]);
293 
294   if (*err == 0) {
295     *basis = CeedBasis_count++;
296     CeedBasis_n++;
297   }
298 }
299 
300 #define fCeedBasisCreateTensorH1 \
301     FORTRAN_NAME(ceedbasiscreatetensorh1, CEEDBASISCREATETENSORH1)
302 void fCeedBasisCreateTensorH1(int *ceed, int *dim, int *ndof, int *P1d,
303                               int *Q1d, const CeedScalar *interp1d, const CeedScalar *grad1d,
304                               const CeedScalar *qref1d, const CeedScalar *qweight1d, int *basis, int *err) {
305   if (CeedBasis_count == CeedBasis_count_max) {
306     CeedBasis_count_max += CeedBasis_count_max/2 + 1;
307     CeedRealloc(CeedBasis_count_max, &CeedBasis_dict);
308   }
309 
310   *err = CeedBasisCreateTensorH1(Ceed_dict[*ceed], *dim, *ndof, *P1d, *Q1d,
311                                  interp1d, grad1d,
312                                  qref1d, qweight1d, &CeedBasis_dict[CeedBasis_count]);
313 
314   if (*err == 0) {
315     *basis = CeedBasis_count++;
316     CeedBasis_n++;
317   }
318 }
319 
320 #define fCeedBasisView FORTRAN_NAME(ceedbasisview, CEEDBASISVIEW)
321 void fCeedBasisView(int *basis, int *err) {
322   *err = CeedBasisView(CeedBasis_dict[*basis], stdout);
323 }
324 
325 #define fCeedQRFactorization \
326     FORTRAN_NAME(ceedqrfactorization, CEEDQRFACTORIZATION)
327 void fCeedQRFactorization(CeedScalar *mat, CeedScalar *tau, int *m, int *n,
328                           int *err) {
329   *err = CeedQRFactorization(mat, tau, *m, *n);
330 }
331 
332 #define fCeedBasisGetColocatedGrad \
333     FORTRAN_NAME(ceedbasisgetcolocatedgrad, CEEDBASISGETCOLOCATEDGRAD)
334 void fCeedBasisGetColocatedGrad(int *basis, CeedScalar *colograd1d,
335                                 int *err) {
336   *err = CeedBasisGetColocatedGrad(CeedBasis_dict[*basis], colograd1d);
337 }
338 
339 #define fCeedBasisApply FORTRAN_NAME(ceedbasisapply, CEEDBASISAPPLY)
340 void fCeedBasisApply(int *basis, int *nelem, int *tmode, int *emode,
341                      const CeedScalar *u, CeedScalar *v, int *err) {
342   *err = CeedBasisApply(CeedBasis_dict[*basis], *nelem, *tmode, *emode, u, v);
343 }
344 
345 #define fCeedBasisGetNumNodes \
346     FORTRAN_NAME(ceedbasisgetnumnodes, CEEDBASISGETNUMNODES)
347 void fCeedBasisGetNumNodes(int *basis, int *P, int *err) {
348   *err = CeedBasisGetNumNodes(CeedBasis_dict[*basis], P);
349 }
350 
351 #define fCeedBasisGetNumQuadraturePoints \
352     FORTRAN_NAME(ceedbasisgetnumquadraturepoints, CEEDBASISGETNUMQUADRATUREPOINTS)
353 void fCeedBasisGetNumQuadraturePoints(int *basis, int *Q, int *err) {
354   *err = CeedBasisGetNumQuadraturePoints(CeedBasis_dict[*basis], Q);
355 }
356 
357 #define fCeedBasisDestroy FORTRAN_NAME(ceedbasisdestroy,CEEDBASISDESTROY)
358 void fCeedBasisDestroy(int *basis, int *err) {
359   *err = CeedBasisDestroy(&CeedBasis_dict[*basis]);
360 
361   if (*err == 0) {
362     CeedBasis_n--;
363     if (CeedBasis_n == 0) {
364       CeedFree(&CeedBasis_dict);
365       CeedBasis_count = 0;
366       CeedBasis_count_max = 0;
367     }
368   }
369 }
370 
371 #define fCeedGaussQuadrature FORTRAN_NAME(ceedgaussquadrature, CEEDGAUSSQUADRATURE)
372 void fCeedGaussQuadrature(int *Q, CeedScalar *qref1d, CeedScalar *qweight1d,
373                           int *err) {
374   *err = CeedGaussQuadrature(*Q, qref1d, qweight1d);
375 }
376 
377 #define fCeedLobattoQuadrature \
378     FORTRAN_NAME(ceedlobattoquadrature, CEEDLOBATTOQUADRATURE)
379 void fCeedLobattoQuadrature(int *Q, CeedScalar *qref1d, CeedScalar *qweight1d,
380                             int *err) {
381   *err = CeedLobattoQuadrature(*Q, qref1d, qweight1d);
382 }
383 
384 static CeedQFunction *CeedQFunction_dict = NULL;
385 static int CeedQFunction_count = 0;
386 static int CeedQFunction_n = 0;
387 static int CeedQFunction_count_max = 0;
388 
389 struct fContext {
390   void (*f)(void *ctx, int *nq,
391             const CeedScalar *u,const CeedScalar *u1,const CeedScalar *u2,
392             const CeedScalar *u3,
393             const CeedScalar *u4,const CeedScalar *u5,const CeedScalar *u6,
394             const CeedScalar *u7,
395             const CeedScalar *u8,const CeedScalar *u9,const CeedScalar *u10,
396             const CeedScalar *u11,
397             const CeedScalar *u12,const CeedScalar *u13,const CeedScalar *u14,
398             const CeedScalar *u15,
399             CeedScalar *v,CeedScalar *v1, CeedScalar *v2,CeedScalar *v3,
400             CeedScalar *v4,CeedScalar *v5, CeedScalar *v6,CeedScalar *v7,
401             CeedScalar *v8,CeedScalar *v9, CeedScalar *v10,CeedScalar *v11,
402             CeedScalar *v12,CeedScalar *v13, CeedScalar *v14,CeedScalar *v15, int *err);
403   void *innerctx;
404 };
405 
406 static int CeedQFunctionFortranStub(void *ctx, int nq,
407                                     const CeedScalar *const *u, CeedScalar *const *v) {
408   struct fContext *fctx = ctx;
409   int ierr;
410 
411   CeedScalar ctx_=1.0;
412   fctx->f((void*)&ctx_,&nq,u[0],u[1],u[2],u[3],u[4],u[5],u[6],u[7],
413           u[8],u[9],u[10],u[11],u[12],u[13],u[14],u[15],
414           v[0],v[1],v[2],v[3],v[4],v[5],v[6],v[7],
415           v[8],v[9],v[10],v[11],v[12],v[13],v[14],v[15],&ierr);
416   return ierr;
417 }
418 
419 #define fCeedQFunctionCreateInterior \
420     FORTRAN_NAME(ceedqfunctioncreateinterior, CEEDQFUNCTIONCREATEINTERIOR)
421 void fCeedQFunctionCreateInterior(int* ceed, int* vlength,
422                                   void (*f)(void *ctx, int *nq,
423                                       const CeedScalar *u,const CeedScalar *u1,const CeedScalar *u2,
424                                       const CeedScalar *u3,
425                                       const CeedScalar *u4,const CeedScalar *u5,const CeedScalar *u6,
426                                       const CeedScalar *u7,
427                                       const CeedScalar *u8,const CeedScalar *u9,const CeedScalar *u10,
428                                       const CeedScalar *u11,
429                                       const CeedScalar *u12,const CeedScalar *u13,const CeedScalar *u14,
430                                       const CeedScalar *u15,
431                                       CeedScalar *v,CeedScalar *v1, CeedScalar *v2,CeedScalar *v3,
432                                       CeedScalar *v4,CeedScalar *v5, CeedScalar *v6,CeedScalar *v7,
433                                       CeedScalar *v8,CeedScalar *v9, CeedScalar *v10,CeedScalar *v11,
434                                       CeedScalar *v12,CeedScalar *v13, CeedScalar *v14,CeedScalar *v15,
435                                       int *err),
436                                   const char *focca, int *qf, int *err) {
437   if (CeedQFunction_count == CeedQFunction_count_max) {
438     CeedQFunction_count_max += CeedQFunction_count_max/2 + 1;
439     CeedRealloc(CeedQFunction_count_max, &CeedQFunction_dict);
440   }
441 
442   CeedQFunction *qf_ = &CeedQFunction_dict[CeedQFunction_count];
443   *err = CeedQFunctionCreateInterior(Ceed_dict[*ceed], *vlength,
444                                      CeedQFunctionFortranStub,focca, qf_);
445 
446   if (*err == 0) {
447     *qf = CeedQFunction_count++;
448     CeedQFunction_n++;
449   }
450 
451   struct fContext *fctx;
452   *err = CeedMalloc(1, &fctx);
453   if (*err) return;
454   fctx->f = f; fctx->innerctx = NULL;
455 
456   *err = CeedQFunctionSetContext(*qf_, fctx, sizeof(struct fContext));
457 
458 }
459 
460 #define fCeedQFunctionAddInput \
461     FORTRAN_NAME(ceedqfunctionaddinput,CEEDQFUNCTIONADDINPUT)
462 void fCeedQFunctionAddInput(int *qf, const char *fieldname,
463                             CeedInt *ncomp, CeedEvalMode *emode, int *err) {
464   CeedQFunction qf_ = CeedQFunction_dict[*qf];
465 
466   *err = CeedQFunctionAddInput(qf_, fieldname, *ncomp, *emode);
467 }
468 
469 #define fCeedQFunctionAddOutput \
470     FORTRAN_NAME(ceedqfunctionaddoutput,CEEDQFUNCTIONADDOUTPUT)
471 void fCeedQFunctionAddOutput(int *qf, const char *fieldname,
472                              CeedInt *ncomp, CeedEvalMode *emode, int *err) {
473   CeedQFunction qf_ = CeedQFunction_dict[*qf];
474 
475   *err = CeedQFunctionAddOutput(qf_, fieldname, *ncomp, *emode);
476 }
477 
478 #define fCeedQFunctionApply \
479     FORTRAN_NAME(ceedqfunctionapply,CEEDQFUNCTIONAPPLY)
480 //TODO Need Fixing, double pointer
481 void fCeedQFunctionApply(int *qf, int *Q,
482                          const CeedScalar *u,const CeedScalar *u1,const CeedScalar *u2,
483                          const CeedScalar *u3,
484                          const CeedScalar *u4,const CeedScalar *u5,const CeedScalar *u6,
485                          const CeedScalar *u7,
486                          const CeedScalar *u8,const CeedScalar *u9,const CeedScalar *u10,
487                          const CeedScalar *u11,
488                          const CeedScalar *u12,const CeedScalar *u13,const CeedScalar *u14,
489                          const CeedScalar *u15,
490                          CeedScalar *v,CeedScalar *v1, CeedScalar *v2,CeedScalar *v3,
491                          CeedScalar *v4,CeedScalar *v5, CeedScalar *v6,CeedScalar *v7,
492                          CeedScalar *v8,CeedScalar *v9, CeedScalar *v10,CeedScalar *v11,
493                          CeedScalar *v12,CeedScalar *v13, CeedScalar *v14,CeedScalar *v15, int *err) {
494   CeedQFunction qf_ = CeedQFunction_dict[*qf];
495   const CeedScalar **in;
496   *err = CeedCalloc(16, &in);
497   if (*err) return;
498   in[0] = u;
499   in[1] = u1;
500   in[2] = u2;
501   in[3] = u3;
502   in[4] = u4;
503   in[5] = u5;
504   in[6] = u6;
505   in[7] = u7;
506   in[8] = u8;
507   in[9] = u9;
508   in[10] = u10;
509   in[11] = u11;
510   in[12] = u12;
511   in[13] = u13;
512   in[14] = u14;
513   in[15] = u15;
514   CeedScalar **out;
515   *err = CeedCalloc(16, &out);
516   if (*err) return;
517   out[0] = v;
518   out[1] = v1;
519   out[2] = v2;
520   out[3] = v3;
521   out[4] = v4;
522   out[5] = v5;
523   out[6] = v6;
524   out[7] = v7;
525   out[8] = v8;
526   out[9] = v9;
527   out[10] = v10;
528   out[11] = v11;
529   out[12] = v12;
530   out[13] = v13;
531   out[14] = v14;
532   out[15] = v15;
533   *err = CeedQFunctionApply(qf_, *Q, (const CeedScalar * const*)in, out);
534   if (*err) return;
535 
536   *err = CeedFree(&in);
537   if (*err) return;
538   *err = CeedFree(&out);
539 }
540 
541 #define fCeedQFunctionDestroy \
542     FORTRAN_NAME(ceedqfunctiondestroy,CEEDQFUNCTIONDESTROY)
543 void fCeedQFunctionDestroy(int *qf, int *err) {
544   CeedFree(&CeedQFunction_dict[*qf]->ctx);
545   *err = CeedQFunctionDestroy(&CeedQFunction_dict[*qf]);
546 
547   if (*err) return;
548   CeedQFunction_n--;
549   if (CeedQFunction_n == 0) {
550     *err = CeedFree(&CeedQFunction_dict);
551     CeedQFunction_count = 0;
552     CeedQFunction_count_max = 0;
553   }
554 }
555 
556 static CeedOperator *CeedOperator_dict = NULL;
557 static int CeedOperator_count = 0;
558 static int CeedOperator_n = 0;
559 static int CeedOperator_count_max = 0;
560 
561 #define fCeedOperatorCreate \
562     FORTRAN_NAME(ceedoperatorcreate, CEEDOPERATORCREATE)
563 void fCeedOperatorCreate(int* ceed,
564                          int* qf, int* dqf, int* dqfT, int *op, int *err) {
565   if (CeedOperator_count == CeedOperator_count_max)
566     CeedOperator_count_max += CeedOperator_count_max/2 + 1,
567                               CeedOperator_dict =
568                                 realloc(CeedOperator_dict, sizeof(CeedOperator)*CeedOperator_count_max);
569 
570   CeedOperator *op_ = &CeedOperator_dict[CeedOperator_count];
571 
572   CeedQFunction dqf_  = NULL, dqfT_ = NULL;
573   if (*dqf  != FORTRAN_NULL) dqf_  = CeedQFunction_dict[*dqf ];
574   if (*dqfT != FORTRAN_NULL) dqfT_ = CeedQFunction_dict[*dqfT];
575 
576   *err = CeedOperatorCreate(Ceed_dict[*ceed], CeedQFunction_dict[*qf], dqf_,
577                             dqfT_, op_);
578   if (*err) return;
579   *op = CeedOperator_count++;
580   CeedOperator_n++;
581 }
582 
583 #define fCeedOperatorSetField \
584     FORTRAN_NAME(ceedoperatorsetfield,CEEDOPERATORSETFIELD)
585 void fCeedOperatorSetField(int *op, const char *fieldname,
586                            int *r, int *b, int *v, int *err) {
587   CeedElemRestriction r_;
588   CeedBasis b_;
589   CeedVector v_;
590 
591   CeedOperator op_ = CeedOperator_dict[*op];
592 
593   if (*r == FORTRAN_NULL) {
594     r_ = NULL;
595   } else {
596     r_ = CeedElemRestriction_dict[*r];
597   }
598   if (*b == FORTRAN_NULL) {
599     b_ = NULL;
600   } else if (*b == FORTRAN_BASIS_COLOCATED) {
601     b_ = CEED_BASIS_COLOCATED;
602   } else {
603     b_ = CeedBasis_dict[*b];
604   }
605   if (*v == FORTRAN_NULL) {
606     v_ = NULL;
607   } else if (*v == FORTRAN_VECTOR_ACTIVE) {
608     v_ = CEED_VECTOR_ACTIVE;
609   } else if (*v == FORTRAN_VECTOR_NONE) {
610     v_ = CEED_VECTOR_NONE;
611   } else {
612     v_ = CeedVector_dict[*v];
613   }
614 
615   *err = CeedOperatorSetField(op_, fieldname, r_, b_, v_);
616 }
617 
618 #define fCeedOperatorApply FORTRAN_NAME(ceedoperatorapply, CEEDOPERATORAPPLY)
619 void fCeedOperatorApply(int *op, int *ustatevec,
620                         int *resvec, int *rqst, int *err) {
621   CeedVector ustatevec_ = *ustatevec == FORTRAN_NULL
622     ? NULL : CeedVector_dict[*ustatevec];
623   CeedVector resvec_ = *resvec == FORTRAN_NULL
624     ? NULL : CeedVector_dict[*resvec];
625 
626   int createRequest = 1;
627   // Check if input is CEED_REQUEST_ORDERED(-2) or CEED_REQUEST_IMMEDIATE(-1)
628   if (*rqst == -1 || *rqst == -2) {
629     createRequest = 0;
630   }
631 
632   if (createRequest && CeedRequest_count == CeedRequest_count_max) {
633     CeedRequest_count_max += CeedRequest_count_max/2 + 1;
634     CeedRealloc(CeedRequest_count_max, &CeedRequest_dict);
635   }
636 
637   CeedRequest *rqst_;
638   if (*rqst == -1) rqst_ = CEED_REQUEST_IMMEDIATE;
639   else if (*rqst == -2) rqst_ = CEED_REQUEST_ORDERED;
640   else rqst_ = &CeedRequest_dict[CeedRequest_count];
641 
642   *err = CeedOperatorApply(CeedOperator_dict[*op],
643                            ustatevec_, resvec_, rqst_);
644   if (*err) return;
645   if (createRequest) {
646     *rqst = CeedRequest_count++;
647     CeedRequest_n++;
648   }
649 }
650 
651 #define fCeedOperatorApplyJacobian \
652     FORTRAN_NAME(ceedoperatorapplyjacobian, CEEDOPERATORAPPLYJACOBIAN)
653 void fCeedOperatorApplyJacobian(int *op, int *qdatavec, int *ustatevec,
654                                 int *dustatevec, int *dresvec, int *rqst, int *err) {
655 // TODO Uncomment this when CeedOperatorApplyJacobian is implemented
656 //  *err = CeedOperatorApplyJacobian(CeedOperator_dict[*op], CeedVector_dict[*qdatavec],
657 //             CeedVector_dict[*ustatevec], CeedVector_dict[*dustatevec],
658 //             CeedVector_dict[*dresvec], &CeedRequest_dict[*rqst]);
659 }
660 
661 #define fCeedOperatorDestroy \
662     FORTRAN_NAME(ceedoperatordestroy, CEEDOPERATORDESTROY)
663 void fCeedOperatorDestroy(int *op, int *err) {
664   *err = CeedOperatorDestroy(&CeedOperator_dict[*op]);
665   if (*err) return;
666   CeedOperator_n--;
667   if (CeedOperator_n == 0) {
668     *err = CeedFree(&CeedOperator_dict);
669     CeedOperator_count = 0;
670     CeedOperator_count_max = 0;
671   }
672 }
673