xref: /libCEED/rust/libceed-sys/c-src/backends/cuda-gen/ceed-cuda-gen-operator.c (revision 3d576824e8d990e1f48c6609089904bee9170514)
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 <stddef.h>
20 #include "ceed-cuda-gen.h"
21 #include "ceed-cuda-gen-operator-build.h"
22 #include "../cuda/ceed-cuda.h"
23 
24 //------------------------------------------------------------------------------
25 // Destroy operator
26 //------------------------------------------------------------------------------
27 static int CeedOperatorDestroy_Cuda_gen(CeedOperator op) {
28   int ierr;
29   CeedOperator_Cuda_gen *impl;
30   ierr = CeedOperatorGetData(op, &impl); CeedChk(ierr);
31   ierr = CeedFree(&impl); CeedChk(ierr);
32   return 0;
33 }
34 
35 //------------------------------------------------------------------------------
36 // Apply and add to output
37 //------------------------------------------------------------------------------
38 static int CeedOperatorApplyAdd_Cuda_gen(CeedOperator op, CeedVector invec,
39     CeedVector outvec, CeedRequest *request) {
40   int ierr;
41   Ceed ceed;
42   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
43   CeedOperator_Cuda_gen *data;
44   ierr = CeedOperatorGetData(op, &data); CeedChk(ierr);
45   CeedQFunction qf;
46   CeedQFunction_Cuda_gen *qf_data;
47   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
48   ierr = CeedQFunctionGetData(qf, &qf_data); CeedChk(ierr);
49   CeedInt nelem, numinputfields, numoutputfields;
50   ierr = CeedOperatorGetNumElements(op, &nelem); CeedChk(ierr);
51   ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
52   CeedChk(ierr);
53   CeedOperatorField *opinputfields, *opoutputfields;
54   ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields);
55   CeedChk(ierr);
56   CeedQFunctionField *qfinputfields, *qfoutputfields;
57   ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields);
58   CeedChk(ierr);
59   CeedEvalMode emode;
60   CeedVector vec, outvecs[16] = {};
61 
62   //Creation of the operator
63   ierr = CeedCudaGenOperatorBuild(op); CeedChk(ierr);
64 
65   // Input vectors
66   for (CeedInt i = 0; i < numinputfields; i++) {
67     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
68     CeedChk(ierr);
69     if (emode == CEED_EVAL_WEIGHT) { // Skip
70       data->fields.in[i] = NULL;
71     } else {
72       // Get input vector
73       ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
74       if (vec == CEED_VECTOR_ACTIVE) vec = invec;
75       ierr = CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, &data->fields.in[i]);
76       CeedChk(ierr);
77     }
78   }
79 
80   // Output vectors
81   for (CeedInt i = 0; i < numoutputfields; i++) {
82     ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode);
83     CeedChk(ierr);
84     if (emode == CEED_EVAL_WEIGHT) { // Skip
85       data->fields.out[i] = NULL;
86     } else {
87       // Get output vector
88       ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
89       if (vec == CEED_VECTOR_ACTIVE) vec = outvec;
90       outvecs[i] = vec;
91       // Check for multiple output modes
92       CeedInt index = -1;
93       for (CeedInt j = 0; j < i; j++) {
94         if (vec == outvecs[j]) {
95           index = j;
96           break;
97         }
98       }
99       if (index == -1) {
100         ierr = CeedVectorGetArray(vec, CEED_MEM_DEVICE, &data->fields.out[i]);
101         CeedChk(ierr);
102       } else {
103         data->fields.out[i] = data->fields.out[index];
104       }
105     }
106   }
107 
108   // Get context data
109   CeedQFunctionContext ctx;
110   ierr = CeedQFunctionGetInnerContext(qf, &ctx); CeedChk(ierr);
111   if (ctx) {
112     ierr = CeedQFunctionContextGetData(ctx, CEED_MEM_DEVICE, &qf_data->d_c);
113     CeedChk(ierr);
114   }
115 
116   // Apply operator
117   void *opargs[] = {(void *) &nelem, &qf_data->d_c, &data->indices,
118                     &data->fields, &data->B, &data->G, &data->W
119                    };
120   const CeedInt dim = data->dim;
121   const CeedInt Q1d = data->Q1d;
122   const CeedInt P1d = data->maxP1d;
123   const CeedInt thread1d = CeedIntMax(Q1d, P1d);
124   if (dim==1) {
125     const CeedInt elemsPerBlock = 32;
126     CeedInt grid = nelem/elemsPerBlock + ( (nelem/elemsPerBlock*elemsPerBlock<nelem)
127                                            ? 1 : 0 );
128     CeedInt sharedMem = elemsPerBlock*thread1d*sizeof(CeedScalar);
129     ierr = CeedRunKernelDimSharedCuda(ceed, data->op, grid, thread1d, 1,
130                                       elemsPerBlock, sharedMem, opargs);
131   } else if (dim==2) {
132     const CeedInt elemsPerBlock = thread1d<4? 16 : 2;
133     CeedInt grid = nelem/elemsPerBlock + ( (nelem/elemsPerBlock*elemsPerBlock<nelem)
134                                            ? 1 : 0 );
135     CeedInt sharedMem = elemsPerBlock*thread1d*thread1d*sizeof(CeedScalar);
136     ierr = CeedRunKernelDimSharedCuda(ceed, data->op, grid, thread1d, thread1d,
137                                       elemsPerBlock, sharedMem, opargs);
138   } else if (dim==3) {
139     const CeedInt elemsPerBlock = thread1d<6? 4 : (thread1d<8? 2 : 1);
140     CeedInt grid = nelem/elemsPerBlock + ( (nelem/elemsPerBlock*elemsPerBlock<nelem)
141                                            ? 1 : 0 );
142     CeedInt sharedMem = elemsPerBlock*thread1d*thread1d*sizeof(CeedScalar);
143     ierr = CeedRunKernelDimSharedCuda(ceed, data->op, grid, thread1d, thread1d,
144                                       elemsPerBlock, sharedMem, opargs);
145   }
146   CeedChk(ierr);
147 
148   // Restore input arrays
149   for (CeedInt i = 0; i < numinputfields; i++) {
150     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
151     CeedChk(ierr);
152     if (emode == CEED_EVAL_WEIGHT) { // Skip
153     } else {
154       ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
155       if (vec == CEED_VECTOR_ACTIVE) vec = invec;
156       ierr = CeedVectorRestoreArrayRead(vec, &data->fields.in[i]);
157       CeedChk(ierr);
158     }
159   }
160 
161   // Restore output arrays
162   for (CeedInt i = 0; i < numoutputfields; i++) {
163     ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode);
164     CeedChk(ierr);
165     if (emode == CEED_EVAL_WEIGHT) { // Skip
166     } else {
167       ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
168       if (vec == CEED_VECTOR_ACTIVE) vec = outvec;
169       // Check for multiple output modes
170       CeedInt index = -1;
171       for (CeedInt j = 0; j < i; j++) {
172         if (vec == outvecs[j]) {
173           index = j;
174           break;
175         }
176       }
177       if (index == -1) {
178         ierr = CeedVectorRestoreArray(vec, &data->fields.out[i]);
179         CeedChk(ierr);
180       }
181     }
182   }
183 
184   // Restore context data
185   if (ctx) {
186     ierr = CeedQFunctionContextRestoreData(ctx, &qf_data->d_c);
187     CeedChk(ierr);
188   }
189   return 0;
190 }
191 
192 //------------------------------------------------------------------------------
193 // Create FDM element inverse not supported
194 //------------------------------------------------------------------------------
195 static int CeedOperatorCreateFDMElementInverse_Cuda(CeedOperator op) {
196   // LCOV_EXCL_START
197   int ierr;
198   Ceed ceed;
199   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
200   return CeedError(ceed, 1, "Backend does not implement FDM inverse creation");
201   // LCOV_EXCL_STOP
202 }
203 
204 //------------------------------------------------------------------------------
205 // Create operator
206 //------------------------------------------------------------------------------
207 int CeedOperatorCreate_Cuda_gen(CeedOperator op) {
208   int ierr;
209   Ceed ceed;
210   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
211   CeedOperator_Cuda_gen *impl;
212 
213   ierr = CeedCalloc(1, &impl); CeedChk(ierr);
214   ierr = CeedOperatorSetData(op, impl); CeedChk(ierr);
215 
216   ierr = CeedSetBackendFunction(ceed, "Operator", op, "CreateFDMElementInverse",
217                                 CeedOperatorCreateFDMElementInverse_Cuda);
218   CeedChk(ierr);
219   ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd",
220                                 CeedOperatorApplyAdd_Cuda_gen); CeedChk(ierr);
221   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy",
222                                 CeedOperatorDestroy_Cuda_gen); CeedChk(ierr);
223   return 0;
224 }
225 //------------------------------------------------------------------------------
226