xref: /libCEED/rust/libceed-sys/c-src/backends/sycl-shared/ceed-sycl-shared-basis.sycl.cpp (revision 7f96cfa0bc215ede2af32822e69abcbe2cbc7f79)
1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED:  http://github.com/ceed
7 
8 #include <ceed/backend.h>
9 #include <ceed/ceed.h>
10 #include <ceed/jit-tools.h>
11 
12 #include <map>
13 #include <string_view>
14 #include <sycl/sycl.hpp>
15 
16 #include "../sycl/ceed-sycl-compile.hpp"
17 #include "ceed-sycl-shared.hpp"
18 
19 //------------------------------------------------------------------------------
20 // Compute the local range of for basis kernels
21 //------------------------------------------------------------------------------
22 static int ComputeLocalRange(Ceed ceed, CeedInt dim, CeedInt thread_1d, CeedInt *local_range, CeedInt max_group_size = 256) {
23   local_range[0]               = thread_1d;
24   local_range[1]               = (dim > 1) ? thread_1d : 1;
25   const CeedInt min_group_size = local_range[0] * local_range[1];
26 
27   CeedCheck(min_group_size <= max_group_size, ceed, CEED_ERROR_BACKEND, "Requested group size is smaller than the required minimum.");
28 
29   local_range[2] = max_group_size / min_group_size;  // elements per group
30   return CEED_ERROR_SUCCESS;
31 }
32 
33 //------------------------------------------------------------------------------
34 // Apply basis
35 //------------------------------------------------------------------------------
36 int CeedBasisApplyTensor_Sycl_shared(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u,
37                                      CeedVector v) {
38   Ceed                   ceed;
39   Ceed_Sycl             *ceed_Sycl;
40   const CeedScalar      *d_u;
41   CeedScalar            *d_v;
42   CeedBasis_Sycl_shared *impl;
43 
44   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
45   CeedCallBackend(CeedGetData(ceed, &ceed_Sycl));
46   CeedCallBackend(CeedBasisGetData(basis, &impl));
47 
48   // Get read/write access to u, v
49   if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u));
50   else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode");
51   CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v));
52 
53   // Apply basis operation
54   switch (eval_mode) {
55     case CEED_EVAL_INTERP: {
56       CeedInt       *lrange         = impl->interp_local_range;
57       const CeedInt &elem_per_group = lrange[2];
58       const CeedInt  group_count    = (num_elem / elem_per_group) + !!(num_elem % elem_per_group);
59       //-----------
60       sycl::range<3>    local_range(lrange[2], lrange[1], lrange[0]);
61       sycl::range<3>    global_range(group_count * lrange[2], lrange[1], lrange[0]);
62       sycl::nd_range<3> kernel_range(global_range, local_range);
63       //-----------
64       sycl::kernel *interp_kernel = (t_mode == CEED_TRANSPOSE) ? impl->interp_transpose_kernel : impl->interp_kernel;
65 
66       // Order queue
67       sycl::event e = ceed_Sycl->sycl_queue.ext_oneapi_submit_barrier();
68 
69       ceed_Sycl->sycl_queue.submit([&](sycl::handler &cgh) {
70         cgh.depends_on(e);
71         cgh.set_args(num_elem, impl->d_interp_1d, d_u, d_v);
72         cgh.parallel_for(kernel_range, *interp_kernel);
73       });
74 
75     } break;
76     case CEED_EVAL_GRAD: {
77       CeedInt       *lrange         = impl->grad_local_range;
78       const CeedInt &elem_per_group = lrange[2];
79       const CeedInt  group_count    = (num_elem / elem_per_group) + !!(num_elem % elem_per_group);
80       //-----------
81       sycl::range<3>    local_range(lrange[2], lrange[1], lrange[0]);
82       sycl::range<3>    global_range(group_count * lrange[2], lrange[1], lrange[0]);
83       sycl::nd_range<3> kernel_range(global_range, local_range);
84       //-----------
85       sycl::kernel     *grad_kernel = (t_mode == CEED_TRANSPOSE) ? impl->grad_transpose_kernel : impl->grad_kernel;
86       const CeedScalar *d_grad_1d   = (impl->d_collo_grad_1d) ? impl->d_collo_grad_1d : impl->d_grad_1d;
87       // Order queue
88       sycl::event e = ceed_Sycl->sycl_queue.ext_oneapi_submit_barrier();
89 
90       ceed_Sycl->sycl_queue.submit([&](sycl::handler &cgh) {
91         cgh.depends_on(e);
92         cgh.set_args(num_elem, impl->d_interp_1d, d_grad_1d, d_u, d_v);
93         cgh.parallel_for(kernel_range, *grad_kernel);
94       });
95     } break;
96     case CEED_EVAL_WEIGHT: {
97       CeedInt       *lrange         = impl->weight_local_range;
98       const CeedInt &elem_per_group = lrange[2];
99       const CeedInt  group_count    = (num_elem / elem_per_group) + !!(num_elem % elem_per_group);
100       //-----------
101       sycl::range<3>    local_range(lrange[2], lrange[1], lrange[0]);
102       sycl::range<3>    global_range(group_count * lrange[2], lrange[1], lrange[0]);
103       sycl::nd_range<3> kernel_range(global_range, local_range);
104       //-----------
105       // Order queue
106       sycl::event e = ceed_Sycl->sycl_queue.ext_oneapi_submit_barrier();
107 
108       ceed_Sycl->sycl_queue.submit([&](sycl::handler &cgh) {
109         cgh.depends_on(e);
110         cgh.set_args(num_elem, impl->d_q_weight_1d, d_v);
111         cgh.parallel_for(kernel_range, *(impl->weight_kernel));
112       });
113     } break;
114     case CEED_EVAL_NONE: /* handled separately below */
115       break;
116     // LCOV_EXCL_START
117     case CEED_EVAL_DIV:
118     case CEED_EVAL_CURL:
119       return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]);
120       // LCOV_EXCL_STOP
121   }
122 
123   // Restore vectors, cover CEED_EVAL_NONE
124   CeedCallBackend(CeedVectorRestoreArray(v, &d_v));
125   if (eval_mode == CEED_EVAL_NONE) CeedCallBackend(CeedVectorSetArray(v, CEED_MEM_DEVICE, CEED_COPY_VALUES, (CeedScalar *)d_u));
126   if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u));
127 
128   return CEED_ERROR_SUCCESS;
129 }
130 
131 //------------------------------------------------------------------------------
132 // Destroy basis
133 //------------------------------------------------------------------------------
134 static int CeedBasisDestroy_Sycl_shared(CeedBasis basis) {
135   Ceed                   ceed;
136   Ceed_Sycl             *data;
137   CeedBasis_Sycl_shared *impl;
138 
139   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
140   CeedCallBackend(CeedBasisGetData(basis, &impl));
141   CeedCallBackend(CeedGetData(ceed, &data));
142   CeedCallSycl(ceed, data->sycl_queue.wait_and_throw());
143   CeedCallSycl(ceed, sycl::free(impl->d_q_weight_1d, data->sycl_context));
144   CeedCallSycl(ceed, sycl::free(impl->d_interp_1d, data->sycl_context));
145   CeedCallSycl(ceed, sycl::free(impl->d_grad_1d, data->sycl_context));
146   CeedCallSycl(ceed, sycl::free(impl->d_collo_grad_1d, data->sycl_context));
147 
148   delete impl->interp_kernel;
149   delete impl->interp_transpose_kernel;
150   delete impl->grad_kernel;
151   delete impl->grad_transpose_kernel;
152   delete impl->weight_kernel;
153   delete impl->sycl_module;
154 
155   CeedCallBackend(CeedFree(&impl));
156   return CEED_ERROR_SUCCESS;
157 }
158 
159 //------------------------------------------------------------------------------
160 // Create tensor basis
161 // TODO: Refactor
162 //------------------------------------------------------------------------------
163 int CeedBasisCreateTensorH1_Sycl_shared(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d,
164                                         const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) {
165   Ceed                   ceed;
166   Ceed_Sycl             *data;
167   char                  *basis_kernel_source;
168   const char            *basis_kernel_path;
169   CeedInt                num_comp;
170   CeedBasis_Sycl_shared *impl;
171 
172   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
173   CeedCallBackend(CeedCalloc(1, &impl));
174   CeedCallBackend(CeedGetData(ceed, &data));
175   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
176 
177   const CeedInt thread_1d = CeedIntMax(Q_1d, P_1d);
178   const CeedInt num_nodes = CeedIntPow(P_1d, dim);
179   const CeedInt num_qpts  = CeedIntPow(Q_1d, dim);
180 
181   CeedInt *interp_lrange = impl->interp_local_range;
182 
183   CeedCallBackend(ComputeLocalRange(ceed, dim, thread_1d, interp_lrange));
184   const CeedInt interp_group_size = interp_lrange[0] * interp_lrange[1] * interp_lrange[2];
185 
186   CeedInt *grad_lrange = impl->grad_local_range;
187 
188   CeedCallBackend(ComputeLocalRange(ceed, dim, thread_1d, grad_lrange));
189   const CeedInt grad_group_size = grad_lrange[0] * grad_lrange[1] * grad_lrange[2];
190 
191   CeedCallBackend(ComputeLocalRange(ceed, dim, Q_1d, impl->weight_local_range));
192 
193   // Copy basis data to GPU
194   CeedCallSycl(ceed, impl->d_q_weight_1d = sycl::malloc_device<CeedScalar>(Q_1d, data->sycl_device, data->sycl_context));
195   sycl::event copy_weight = data->sycl_queue.copy<CeedScalar>(q_weight_1d, impl->d_q_weight_1d, Q_1d);
196 
197   const CeedInt interp_length = Q_1d * P_1d;
198   CeedCallSycl(ceed, impl->d_interp_1d = sycl::malloc_device<CeedScalar>(interp_length, data->sycl_device, data->sycl_context));
199   sycl::event copy_interp = data->sycl_queue.copy<CeedScalar>(interp_1d, impl->d_interp_1d, interp_length);
200 
201   CeedCallSycl(ceed, impl->d_grad_1d = sycl::malloc_device<CeedScalar>(interp_length, data->sycl_device, data->sycl_context));
202   sycl::event copy_grad = data->sycl_queue.copy<CeedScalar>(grad_1d, impl->d_grad_1d, interp_length);
203 
204   CeedCallSycl(ceed, sycl::event::wait_and_throw({copy_weight, copy_interp, copy_grad}));
205 
206   // Compute collocated gradient and copy to GPU
207   impl->d_collo_grad_1d          = NULL;
208   const bool has_collocated_grad = (dim == 3) && (Q_1d >= P_1d);
209 
210   if (has_collocated_grad) {
211     CeedScalar   *collo_grad_1d;
212     const CeedInt cgrad_length = Q_1d * Q_1d;
213 
214     CeedCallBackend(CeedMalloc(Q_1d * Q_1d, &collo_grad_1d));
215     CeedCallBackend(CeedBasisGetCollocatedGrad(basis, collo_grad_1d));
216     CeedCallSycl(ceed, impl->d_collo_grad_1d = sycl::malloc_device<CeedScalar>(cgrad_length, data->sycl_device, data->sycl_context));
217     CeedCallSycl(ceed, data->sycl_queue.copy<CeedScalar>(collo_grad_1d, impl->d_collo_grad_1d, cgrad_length).wait_and_throw());
218     CeedCallBackend(CeedFree(&collo_grad_1d));
219   }
220 
221   // ---[Refactor into separate function]------>
222   // Define compile-time constants
223   std::map<std::string, CeedInt> jit_constants;
224   jit_constants["BASIS_DIM"]                 = dim;
225   jit_constants["BASIS_Q_1D"]                = Q_1d;
226   jit_constants["BASIS_P_1D"]                = P_1d;
227   jit_constants["T_1D"]                      = thread_1d;
228   jit_constants["BASIS_NUM_COMP"]            = num_comp;
229   jit_constants["BASIS_NUM_NODES"]           = num_nodes;
230   jit_constants["BASIS_NUM_QPTS"]            = num_qpts;
231   jit_constants["BASIS_HAS_COLLOCATED_GRAD"] = has_collocated_grad;
232   jit_constants["BASIS_INTERP_SCRATCH_SIZE"] = interp_group_size;
233   jit_constants["BASIS_GRAD_SCRATCH_SIZE"]   = grad_group_size;
234 
235   // Load kernel source
236   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/sycl/sycl-shared-basis-tensor.h", &basis_kernel_path));
237   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n");
238   {
239     char *source;
240 
241     CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &source));
242     basis_kernel_source = source;
243   }
244   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete -----\n");
245 
246   // Compile kernels into a kernel bundle
247   CeedCallBackend(CeedBuildModule_Sycl(ceed, basis_kernel_source, &impl->sycl_module, jit_constants));
248 
249   // Load kernel functions
250   CeedCallBackend(CeedGetKernel_Sycl(ceed, impl->sycl_module, "Interp", &impl->interp_kernel));
251   CeedCallBackend(CeedGetKernel_Sycl(ceed, impl->sycl_module, "InterpTranspose", &impl->interp_transpose_kernel));
252   CeedCallBackend(CeedGetKernel_Sycl(ceed, impl->sycl_module, "Grad", &impl->grad_kernel));
253   CeedCallBackend(CeedGetKernel_Sycl(ceed, impl->sycl_module, "GradTranspose", &impl->grad_transpose_kernel));
254   CeedCallBackend(CeedGetKernel_Sycl(ceed, impl->sycl_module, "Weight", &impl->weight_kernel));
255 
256   // Clean-up
257   CeedCallBackend(CeedFree(&basis_kernel_path));
258   CeedCallBackend(CeedFree(&basis_kernel_source));
259   // <---[Refactor into separate function]------
260 
261   CeedCallBackend(CeedBasisSetData(basis, impl));
262 
263   // Register backend functions
264   CeedCallBackend(CeedSetBackendFunctionCpp(ceed, "Basis", basis, "Apply", CeedBasisApplyTensor_Sycl_shared));
265   CeedCallBackend(CeedSetBackendFunctionCpp(ceed, "Basis", basis, "Destroy", CeedBasisDestroy_Sycl_shared));
266   return CEED_ERROR_SUCCESS;
267 }
268 
269 //------------------------------------------------------------------------------
270