xref: /libCEED/rust/libceed-sys/c-src/backends/hip-shared/ceed-hip-shared-basis.c (revision 1dda9c1ab34a6c4075a3e9d5ebc86d8c1194a9e1)
15aed82e4SJeremy L Thompson // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors.
23d8e8822SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
37d8d0e25Snbeams //
43d8e8822SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause
57d8d0e25Snbeams //
63d8e8822SJeremy L Thompson // This file is part of CEED:  http://github.com/ceed
77d8d0e25Snbeams 
849aac155SJeremy L Thompson #include <ceed.h>
9ec3da8bcSJed Brown #include <ceed/backend.h>
10437930d1SJeremy L Thompson #include <ceed/jit-tools.h>
1149aac155SJeremy L Thompson #include <stdbool.h>
123d576824SJeremy L Thompson #include <stddef.h>
13c85e8640SSebastian Grimberg #include <hip/hip_runtime.h>
142b730f8bSJeremy L Thompson 
157fcac036SJeremy L Thompson #include "../hip/ceed-hip-common.h"
167d8d0e25Snbeams #include "../hip/ceed-hip-compile.h"
172b730f8bSJeremy L Thompson #include "ceed-hip-shared.h"
187d8d0e25Snbeams 
197d8d0e25Snbeams //------------------------------------------------------------------------------
209e31c45bSnbeams // Compute a block size based on required minimum threads
219e31c45bSnbeams //------------------------------------------------------------------------------
229e31c45bSnbeams static CeedInt ComputeBlockSizeFromRequirement(const CeedInt required) {
239e31c45bSnbeams   CeedInt maxSize     = 1024;  // Max total threads per block
249e31c45bSnbeams   CeedInt currentSize = 64;    // Start with one group
259e31c45bSnbeams 
269e31c45bSnbeams   while (currentSize < maxSize) {
272b730f8bSJeremy L Thompson     if (currentSize > required) break;
282b730f8bSJeremy L Thompson     else currentSize = currentSize * 2;
299e31c45bSnbeams   }
309e31c45bSnbeams   return currentSize;
319e31c45bSnbeams }
329e31c45bSnbeams 
339e31c45bSnbeams //------------------------------------------------------------------------------
349e31c45bSnbeams // Compute required thread block sizes for basis kernels given P, Q, dim, and
359e201c85SYohann // num_comp (num_comp not currently used, but may be again in other basis
369e201c85SYohann // parallelization options)
379e31c45bSnbeams //------------------------------------------------------------------------------
382b730f8bSJeremy L Thompson static int ComputeBasisThreadBlockSizes(const CeedInt dim, const CeedInt P_1d, const CeedInt Q_1d, const CeedInt num_comp, CeedInt *block_sizes) {
399e31c45bSnbeams   // Note that this will use the same block sizes for all dimensions when compiling,
409e31c45bSnbeams   // but as each basis object is defined for a particular dimension, we will never
419e31c45bSnbeams   // call any kernels except the ones for the dimension for which we have computed the
429e31c45bSnbeams   // block sizes.
43437930d1SJeremy L Thompson   const CeedInt thread_1d = CeedIntMax(P_1d, Q_1d);
44b7453713SJeremy L Thompson 
459e31c45bSnbeams   switch (dim) {
469e31c45bSnbeams     case 1: {
479e31c45bSnbeams       // Interp kernels:
48437930d1SJeremy L Thompson       block_sizes[0] = 256;
499e31c45bSnbeams 
509e31c45bSnbeams       // Grad kernels:
51437930d1SJeremy L Thompson       block_sizes[1] = 256;
529e31c45bSnbeams 
539e31c45bSnbeams       // Weight kernels:
54437930d1SJeremy L Thompson       block_sizes[2] = 256;
559e31c45bSnbeams     } break;
569e31c45bSnbeams     case 2: {
579e31c45bSnbeams       // Interp kernels:
589e201c85SYohann       CeedInt required = thread_1d * thread_1d;
59b7453713SJeremy L Thompson 
609e201c85SYohann       block_sizes[0] = CeedIntMax(256, ComputeBlockSizeFromRequirement(required));
619e31c45bSnbeams 
629e31c45bSnbeams       // Grad kernels: currently use same required minimum threads
639e201c85SYohann       block_sizes[1] = CeedIntMax(256, ComputeBlockSizeFromRequirement(required));
649e31c45bSnbeams 
659e31c45bSnbeams       // Weight kernels:
66437930d1SJeremy L Thompson       required       = CeedIntMax(64, Q_1d * Q_1d);
679e201c85SYohann       block_sizes[2] = CeedIntMax(256, ComputeBlockSizeFromRequirement(required));
689e31c45bSnbeams 
699e31c45bSnbeams     } break;
709e31c45bSnbeams     case 3: {
719e31c45bSnbeams       // Interp kernels:
729e201c85SYohann       CeedInt required = thread_1d * thread_1d;
73b7453713SJeremy L Thompson 
749e201c85SYohann       block_sizes[0] = CeedIntMax(256, ComputeBlockSizeFromRequirement(required));
759e31c45bSnbeams 
769e31c45bSnbeams       // Grad kernels: currently use same required minimum threads
779e201c85SYohann       block_sizes[1] = CeedIntMax(256, ComputeBlockSizeFromRequirement(required));
789e31c45bSnbeams 
799e31c45bSnbeams       // Weight kernels:
80437930d1SJeremy L Thompson       required       = Q_1d * Q_1d * Q_1d;
819e201c85SYohann       block_sizes[2] = CeedIntMax(256, ComputeBlockSizeFromRequirement(required));
829e31c45bSnbeams     }
839e31c45bSnbeams   }
84e15f9bd0SJeremy L Thompson   return CEED_ERROR_SUCCESS;
859e31c45bSnbeams }
869e31c45bSnbeams 
879e31c45bSnbeams //------------------------------------------------------------------------------
887d8d0e25Snbeams // Apply basis
897d8d0e25Snbeams //------------------------------------------------------------------------------
902b730f8bSJeremy L Thompson int CeedBasisApplyTensor_Hip_shared(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u,
917d8d0e25Snbeams                                     CeedVector v) {
927d8d0e25Snbeams   Ceed                  ceed;
936dbfb411Snbeams   Ceed_Hip             *ceed_Hip;
94437930d1SJeremy L Thompson   CeedInt               dim, num_comp;
95b7453713SJeremy L Thompson   const CeedScalar     *d_u;
96b7453713SJeremy L Thompson   CeedScalar           *d_v;
97b7453713SJeremy L Thompson   CeedBasis_Hip_shared *data;
98b7453713SJeremy L Thompson 
99b7453713SJeremy L Thompson   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
100b7453713SJeremy L Thompson   CeedCallBackend(CeedGetData(ceed, &ceed_Hip));
101b7453713SJeremy L Thompson   CeedCallBackend(CeedBasisGetData(basis, &data));
1022b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisGetDimension(basis, &dim));
1032b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
1047d8d0e25Snbeams 
1059ea2cfd9SJeremy L Thompson   // Get read/write access to u, v
1066574a04fSJeremy L Thompson   if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u));
1076574a04fSJeremy L Thompson   else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode");
1082b730f8bSJeremy L Thompson   CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v));
1097d8d0e25Snbeams 
1107d8d0e25Snbeams   // Apply basis operation
111437930d1SJeremy L Thompson   switch (eval_mode) {
1127d8d0e25Snbeams     case CEED_EVAL_INTERP: {
113437930d1SJeremy L Thompson       CeedInt P_1d, Q_1d;
114437930d1SJeremy L Thompson       CeedInt block_size = data->block_sizes[0];
115b7453713SJeremy L Thompson 
1162b730f8bSJeremy L Thompson       CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
1172b730f8bSJeremy L Thompson       CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
118437930d1SJeremy L Thompson       CeedInt thread_1d     = CeedIntMax(Q_1d, P_1d);
1192b730f8bSJeremy L Thompson       void   *interp_args[] = {(void *)&num_elem, &data->d_interp_1d, &d_u, &d_v};
120b7453713SJeremy L Thompson 
1217d8d0e25Snbeams       if (dim == 1) {
122437930d1SJeremy L Thompson         CeedInt elems_per_block = 64 * thread_1d > 256 ? 256 / thread_1d : 64;
123437930d1SJeremy L Thompson         elems_per_block         = elems_per_block > 0 ? elems_per_block : 1;
1242b730f8bSJeremy L Thompson         CeedInt grid            = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0);
125437930d1SJeremy L Thompson         CeedInt shared_mem      = elems_per_block * thread_1d * sizeof(CeedScalar);
126b2165e7aSSebastian Grimberg 
1279e201c85SYohann         if (t_mode == CEED_TRANSPOSE) {
128eb7e6cafSJeremy L Thompson           CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->InterpTranspose, grid, thread_1d, 1, elems_per_block, shared_mem, interp_args));
1299e201c85SYohann         } else {
130eb7e6cafSJeremy L Thompson           CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->Interp, grid, thread_1d, 1, elems_per_block, shared_mem, interp_args));
1319e201c85SYohann         }
1327d8d0e25Snbeams       } else if (dim == 2) {
1339e31c45bSnbeams         // Check if required threads is small enough to do multiple elems
1342b730f8bSJeremy L Thompson         const CeedInt elems_per_block = CeedIntMax(block_size / (thread_1d * thread_1d), 1);
1352b730f8bSJeremy L Thompson         CeedInt       grid            = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0);
1362b730f8bSJeremy L Thompson         CeedInt       shared_mem      = elems_per_block * thread_1d * thread_1d * sizeof(CeedScalar);
137b2165e7aSSebastian Grimberg 
1389e201c85SYohann         if (t_mode == CEED_TRANSPOSE) {
1392b730f8bSJeremy L Thompson           CeedCallBackend(
140eb7e6cafSJeremy L Thompson               CeedRunKernelDimShared_Hip(ceed, data->InterpTranspose, grid, thread_1d, thread_1d, elems_per_block, shared_mem, interp_args));
1419e201c85SYohann         } else {
142eb7e6cafSJeremy L Thompson           CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->Interp, grid, thread_1d, thread_1d, elems_per_block, shared_mem, interp_args));
1439e201c85SYohann         }
1447d8d0e25Snbeams       } else if (dim == 3) {
1452b730f8bSJeremy L Thompson         const CeedInt elems_per_block = CeedIntMax(block_size / (thread_1d * thread_1d), 1);
1462b730f8bSJeremy L Thompson         CeedInt       grid            = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0);
1472b730f8bSJeremy L Thompson         CeedInt       shared_mem      = elems_per_block * thread_1d * thread_1d * sizeof(CeedScalar);
148b2165e7aSSebastian Grimberg 
1499e201c85SYohann         if (t_mode == CEED_TRANSPOSE) {
1502b730f8bSJeremy L Thompson           CeedCallBackend(
151eb7e6cafSJeremy L Thompson               CeedRunKernelDimShared_Hip(ceed, data->InterpTranspose, grid, thread_1d, thread_1d, elems_per_block, shared_mem, interp_args));
1529e201c85SYohann         } else {
153eb7e6cafSJeremy L Thompson           CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->Interp, grid, thread_1d, thread_1d, elems_per_block, shared_mem, interp_args));
1549e201c85SYohann         }
1557d8d0e25Snbeams       }
1567d8d0e25Snbeams     } break;
1577d8d0e25Snbeams     case CEED_EVAL_GRAD: {
158437930d1SJeremy L Thompson       CeedInt P_1d, Q_1d;
159437930d1SJeremy L Thompson       CeedInt block_size = data->block_sizes[1];
160b7453713SJeremy L Thompson 
1612b730f8bSJeremy L Thompson       CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
1622b730f8bSJeremy L Thompson       CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
163437930d1SJeremy L Thompson       CeedInt     thread_1d = CeedIntMax(Q_1d, P_1d);
1649e201c85SYohann       CeedScalar *d_grad_1d = data->d_grad_1d;
165b7453713SJeremy L Thompson 
1669e201c85SYohann       if (data->d_collo_grad_1d) {
1679e201c85SYohann         d_grad_1d = data->d_collo_grad_1d;
1689e201c85SYohann       }
1692b730f8bSJeremy L Thompson       void *grad_args[] = {(void *)&num_elem, &data->d_interp_1d, &d_grad_1d, &d_u, &d_v};
1707d8d0e25Snbeams       if (dim == 1) {
171437930d1SJeremy L Thompson         CeedInt elems_per_block = 64 * thread_1d > 256 ? 256 / thread_1d : 64;
172437930d1SJeremy L Thompson         elems_per_block         = elems_per_block > 0 ? elems_per_block : 1;
1732b730f8bSJeremy L Thompson         CeedInt grid            = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0);
174437930d1SJeremy L Thompson         CeedInt shared_mem      = elems_per_block * thread_1d * sizeof(CeedScalar);
175b2165e7aSSebastian Grimberg 
1769e201c85SYohann         if (t_mode == CEED_TRANSPOSE) {
177eb7e6cafSJeremy L Thompson           CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->GradTranspose, grid, thread_1d, 1, elems_per_block, shared_mem, grad_args));
1789e201c85SYohann         } else {
179eb7e6cafSJeremy L Thompson           CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->Grad, grid, thread_1d, 1, elems_per_block, shared_mem, grad_args));
1809e201c85SYohann         }
1817d8d0e25Snbeams       } else if (dim == 2) {
1829e31c45bSnbeams         // Check if required threads is small enough to do multiple elems
1832b730f8bSJeremy L Thompson         const CeedInt elems_per_block = CeedIntMax(block_size / (thread_1d * thread_1d), 1);
1842b730f8bSJeremy L Thompson         CeedInt       grid            = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0);
1852b730f8bSJeremy L Thompson         CeedInt       shared_mem      = elems_per_block * thread_1d * thread_1d * sizeof(CeedScalar);
186b2165e7aSSebastian Grimberg 
1879e201c85SYohann         if (t_mode == CEED_TRANSPOSE) {
188eb7e6cafSJeremy L Thompson           CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->GradTranspose, grid, thread_1d, thread_1d, elems_per_block, shared_mem, grad_args));
1899e201c85SYohann         } else {
190eb7e6cafSJeremy L Thompson           CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->Grad, grid, thread_1d, thread_1d, elems_per_block, shared_mem, grad_args));
1919e201c85SYohann         }
1927d8d0e25Snbeams       } else if (dim == 3) {
1932b730f8bSJeremy L Thompson         const CeedInt elems_per_block = CeedIntMax(block_size / (thread_1d * thread_1d), 1);
1942b730f8bSJeremy L Thompson         CeedInt       grid            = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0);
1952b730f8bSJeremy L Thompson         CeedInt       shared_mem      = elems_per_block * thread_1d * thread_1d * sizeof(CeedScalar);
196b2165e7aSSebastian Grimberg 
1979e201c85SYohann         if (t_mode == CEED_TRANSPOSE) {
198eb7e6cafSJeremy L Thompson           CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->GradTranspose, grid, thread_1d, thread_1d, elems_per_block, shared_mem, grad_args));
1999e201c85SYohann         } else {
200eb7e6cafSJeremy L Thompson           CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->Grad, grid, thread_1d, thread_1d, elems_per_block, shared_mem, grad_args));
2019e201c85SYohann         }
2027d8d0e25Snbeams       }
2037d8d0e25Snbeams     } break;
2047d8d0e25Snbeams     case CEED_EVAL_WEIGHT: {
205437930d1SJeremy L Thompson       CeedInt Q_1d;
206437930d1SJeremy L Thompson       CeedInt block_size = data->block_sizes[2];
207b7453713SJeremy L Thompson 
208097cc795SJames Wright       CeedCheck(data->d_q_weight_1d, ceed, CEED_ERROR_BACKEND, "%s not supported; q_weights_1d not set", CeedEvalModes[eval_mode]);
2092b730f8bSJeremy L Thompson       CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
210437930d1SJeremy L Thompson       void *weight_args[] = {(void *)&num_elem, (void *)&data->d_q_weight_1d, &d_v};
211b7453713SJeremy L Thompson 
2127d8d0e25Snbeams       if (dim == 1) {
213437930d1SJeremy L Thompson         const CeedInt opt_elems       = block_size / Q_1d;
214437930d1SJeremy L Thompson         const CeedInt elems_per_block = opt_elems > 0 ? opt_elems : 1;
2152b730f8bSJeremy L Thompson         const CeedInt grid_size       = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0);
216b2165e7aSSebastian Grimberg 
217eb7e6cafSJeremy L Thompson         CeedCallBackend(CeedRunKernelDim_Hip(ceed, data->Weight, grid_size, Q_1d, elems_per_block, 1, weight_args));
2187d8d0e25Snbeams       } else if (dim == 2) {
219437930d1SJeremy L Thompson         const CeedInt opt_elems       = block_size / (Q_1d * Q_1d);
220437930d1SJeremy L Thompson         const CeedInt elems_per_block = opt_elems > 0 ? opt_elems : 1;
2212b730f8bSJeremy L Thompson         const CeedInt grid_size       = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0);
222b2165e7aSSebastian Grimberg 
223eb7e6cafSJeremy L Thompson         CeedCallBackend(CeedRunKernelDim_Hip(ceed, data->Weight, grid_size, Q_1d, Q_1d, elems_per_block, weight_args));
2247d8d0e25Snbeams       } else if (dim == 3) {
2259e201c85SYohann         const CeedInt opt_elems       = block_size / (Q_1d * Q_1d);
2269e201c85SYohann         const CeedInt elems_per_block = opt_elems > 0 ? opt_elems : 1;
2272b730f8bSJeremy L Thompson         const CeedInt grid_size       = num_elem / elems_per_block + ((num_elem / elems_per_block * elems_per_block < num_elem) ? 1 : 0);
228b2165e7aSSebastian Grimberg 
229eb7e6cafSJeremy L Thompson         CeedCallBackend(CeedRunKernelDim_Hip(ceed, data->Weight, grid_size, Q_1d, Q_1d, elems_per_block, weight_args));
2307d8d0e25Snbeams       }
2317d8d0e25Snbeams     } break;
2329ea2cfd9SJeremy L Thompson     case CEED_EVAL_NONE: /* handled separately below */
2339ea2cfd9SJeremy L Thompson       break;
2347d8d0e25Snbeams     // LCOV_EXCL_START
2357d8d0e25Snbeams     case CEED_EVAL_DIV:
2367d8d0e25Snbeams     case CEED_EVAL_CURL:
237bcbe1c99SJeremy L Thompson       return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]);
2387d8d0e25Snbeams       // LCOV_EXCL_STOP
2397d8d0e25Snbeams   }
2407d8d0e25Snbeams 
2419ea2cfd9SJeremy L Thompson   // Restore vectors, cover CEED_EVAL_NONE
2422b730f8bSJeremy L Thompson   CeedCallBackend(CeedVectorRestoreArray(v, &d_v));
2439ea2cfd9SJeremy L Thompson   if (eval_mode == CEED_EVAL_NONE) CeedCallBackend(CeedVectorSetArray(v, CEED_MEM_DEVICE, CEED_COPY_VALUES, (CeedScalar *)d_u));
2449ea2cfd9SJeremy L Thompson   if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u));
245e15f9bd0SJeremy L Thompson   return CEED_ERROR_SUCCESS;
2467d8d0e25Snbeams }
2477d8d0e25Snbeams 
2487d8d0e25Snbeams //------------------------------------------------------------------------------
249*1dda9c1aSJeremy L Thompson // Basis apply - tensor AtPoints
250*1dda9c1aSJeremy L Thompson //------------------------------------------------------------------------------
251*1dda9c1aSJeremy L Thompson int CeedBasisApplyAtPoints_Hip_shared(CeedBasis basis, const CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode,
252*1dda9c1aSJeremy L Thompson                                       CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) {
253*1dda9c1aSJeremy L Thompson   Ceed                  ceed;
254*1dda9c1aSJeremy L Thompson   CeedInt               Q_1d, dim, max_num_points = num_points[0];
255*1dda9c1aSJeremy L Thompson   const CeedInt         is_transpose   = t_mode == CEED_TRANSPOSE;
256*1dda9c1aSJeremy L Thompson   const int             max_block_size = 32;
257*1dda9c1aSJeremy L Thompson   const CeedScalar     *d_x, *d_u;
258*1dda9c1aSJeremy L Thompson   CeedScalar           *d_v;
259*1dda9c1aSJeremy L Thompson   CeedBasis_Hip_shared *data;
260*1dda9c1aSJeremy L Thompson 
261*1dda9c1aSJeremy L Thompson   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
262*1dda9c1aSJeremy L Thompson   CeedCallBackend(CeedBasisGetData(basis, &data));
263*1dda9c1aSJeremy L Thompson   CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
264*1dda9c1aSJeremy L Thompson   CeedCallBackend(CeedBasisGetDimension(basis, &dim));
265*1dda9c1aSJeremy L Thompson 
266*1dda9c1aSJeremy L Thompson   // Check uniform number of points per elem
267*1dda9c1aSJeremy L Thompson   for (CeedInt i = 1; i < num_elem; i++) {
268*1dda9c1aSJeremy L Thompson     CeedCheck(max_num_points == num_points[i], ceed, CEED_ERROR_BACKEND,
269*1dda9c1aSJeremy L Thompson               "BasisApplyAtPoints only supported for the same number of points in each element");
270*1dda9c1aSJeremy L Thompson   }
271*1dda9c1aSJeremy L Thompson 
272*1dda9c1aSJeremy L Thompson   // Weight handled separately
273*1dda9c1aSJeremy L Thompson   if (eval_mode == CEED_EVAL_WEIGHT) {
274*1dda9c1aSJeremy L Thompson     CeedCall(CeedVectorSetValue(v, 1.0));
275*1dda9c1aSJeremy L Thompson     return CEED_ERROR_SUCCESS;
276*1dda9c1aSJeremy L Thompson   }
277*1dda9c1aSJeremy L Thompson 
278*1dda9c1aSJeremy L Thompson   // Build kernels if needed
279*1dda9c1aSJeremy L Thompson   if (data->num_points != max_num_points) {
280*1dda9c1aSJeremy L Thompson     CeedInt P_1d;
281*1dda9c1aSJeremy L Thompson 
282*1dda9c1aSJeremy L Thompson     CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
283*1dda9c1aSJeremy L Thompson     data->num_points = max_num_points;
284*1dda9c1aSJeremy L Thompson 
285*1dda9c1aSJeremy L Thompson     // -- Create interp matrix to Chebyshev coefficients
286*1dda9c1aSJeremy L Thompson     if (!data->d_chebyshev_interp_1d) {
287*1dda9c1aSJeremy L Thompson       CeedSize    interp_bytes;
288*1dda9c1aSJeremy L Thompson       CeedScalar *chebyshev_interp_1d;
289*1dda9c1aSJeremy L Thompson 
290*1dda9c1aSJeremy L Thompson       interp_bytes = P_1d * Q_1d * sizeof(CeedScalar);
291*1dda9c1aSJeremy L Thompson       CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d));
292*1dda9c1aSJeremy L Thompson       CeedCall(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d));
293*1dda9c1aSJeremy L Thompson       CeedCallHip(ceed, hipMalloc((void **)&data->d_chebyshev_interp_1d, interp_bytes));
294*1dda9c1aSJeremy L Thompson       CeedCallHip(ceed, hipMemcpy(data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, hipMemcpyHostToDevice));
295*1dda9c1aSJeremy L Thompson       CeedCallBackend(CeedFree(&chebyshev_interp_1d));
296*1dda9c1aSJeremy L Thompson     }
297*1dda9c1aSJeremy L Thompson 
298*1dda9c1aSJeremy L Thompson     // -- Compile kernels
299*1dda9c1aSJeremy L Thompson     char       *basis_kernel_source;
300*1dda9c1aSJeremy L Thompson     const char *basis_kernel_path;
301*1dda9c1aSJeremy L Thompson     CeedInt     num_comp;
302*1dda9c1aSJeremy L Thompson 
303*1dda9c1aSJeremy L Thompson     if (data->moduleAtPoints) CeedCallHip(ceed, hipModuleUnload(data->moduleAtPoints));
304*1dda9c1aSJeremy L Thompson     CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
305*1dda9c1aSJeremy L Thompson     CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/hip/hip-ref-basis-tensor-at-points.h", &basis_kernel_path));
306*1dda9c1aSJeremy L Thompson     CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n");
307*1dda9c1aSJeremy L Thompson     CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source));
308*1dda9c1aSJeremy L Thompson     CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n");
309*1dda9c1aSJeremy L Thompson     CeedCallBackend(CeedCompile_Hip(ceed, basis_kernel_source, &data->moduleAtPoints, 9, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "BASIS_BUF_LEN",
310*1dda9c1aSJeremy L Thompson                                     num_comp * CeedIntPow(Q_1d > P_1d ? Q_1d : P_1d, dim), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp,
311*1dda9c1aSJeremy L Thompson                                     "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "BASIS_NUM_PTS",
312*1dda9c1aSJeremy L Thompson                                     max_num_points, "POINTS_BUFF_LEN",
313*1dda9c1aSJeremy L Thompson                                     max_num_points * CeedIntPow(Q_1d > max_num_points ? Q_1d : max_num_points, dim - 1)));
314*1dda9c1aSJeremy L Thompson     CeedCallBackend(CeedGetKernel_Hip(ceed, data->moduleAtPoints, "InterpAtPoints", &data->InterpAtPoints));
315*1dda9c1aSJeremy L Thompson     CeedCallBackend(CeedGetKernel_Hip(ceed, data->moduleAtPoints, "GradAtPoints", &data->GradAtPoints));
316*1dda9c1aSJeremy L Thompson     CeedCallBackend(CeedFree(&basis_kernel_path));
317*1dda9c1aSJeremy L Thompson     CeedCallBackend(CeedFree(&basis_kernel_source));
318*1dda9c1aSJeremy L Thompson   }
319*1dda9c1aSJeremy L Thompson 
320*1dda9c1aSJeremy L Thompson   // Get read/write access to u, v
321*1dda9c1aSJeremy L Thompson   CeedCallBackend(CeedVectorGetArrayRead(x_ref, CEED_MEM_DEVICE, &d_x));
322*1dda9c1aSJeremy L Thompson   if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u));
323*1dda9c1aSJeremy L Thompson   else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode");
324*1dda9c1aSJeremy L Thompson   CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v));
325*1dda9c1aSJeremy L Thompson 
326*1dda9c1aSJeremy L Thompson   // Clear v for transpose operation
327*1dda9c1aSJeremy L Thompson   if (is_transpose) {
328*1dda9c1aSJeremy L Thompson     CeedSize length;
329*1dda9c1aSJeremy L Thompson 
330*1dda9c1aSJeremy L Thompson     CeedCallBackend(CeedVectorGetLength(v, &length));
331*1dda9c1aSJeremy L Thompson     CeedCallHip(ceed, hipMemset(d_v, 0, length * sizeof(CeedScalar)));
332*1dda9c1aSJeremy L Thompson   }
333*1dda9c1aSJeremy L Thompson 
334*1dda9c1aSJeremy L Thompson   // Basis action
335*1dda9c1aSJeremy L Thompson   switch (eval_mode) {
336*1dda9c1aSJeremy L Thompson     case CEED_EVAL_INTERP: {
337*1dda9c1aSJeremy L Thompson       void         *interp_args[] = {(void *)&num_elem, (void *)&is_transpose, &data->d_chebyshev_interp_1d, &d_x, &d_u, &d_v};
338*1dda9c1aSJeremy L Thompson       const CeedInt block_size    = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size);
339*1dda9c1aSJeremy L Thompson 
340*1dda9c1aSJeremy L Thompson       CeedCallBackend(CeedRunKernel_Hip(ceed, data->InterpAtPoints, num_elem, block_size, interp_args));
341*1dda9c1aSJeremy L Thompson     } break;
342*1dda9c1aSJeremy L Thompson     case CEED_EVAL_GRAD: {
343*1dda9c1aSJeremy L Thompson       void         *grad_args[] = {(void *)&num_elem, (void *)&is_transpose, &data->d_chebyshev_interp_1d, &d_x, &d_u, &d_v};
344*1dda9c1aSJeremy L Thompson       const CeedInt block_size  = max_block_size;
345*1dda9c1aSJeremy L Thompson 
346*1dda9c1aSJeremy L Thompson       CeedCallBackend(CeedRunKernel_Hip(ceed, data->GradAtPoints, num_elem, block_size, grad_args));
347*1dda9c1aSJeremy L Thompson     } break;
348*1dda9c1aSJeremy L Thompson     case CEED_EVAL_WEIGHT:
349*1dda9c1aSJeremy L Thompson     case CEED_EVAL_NONE: /* handled separately below */
350*1dda9c1aSJeremy L Thompson       break;
351*1dda9c1aSJeremy L Thompson     // LCOV_EXCL_START
352*1dda9c1aSJeremy L Thompson     case CEED_EVAL_DIV:
353*1dda9c1aSJeremy L Thompson     case CEED_EVAL_CURL:
354*1dda9c1aSJeremy L Thompson       return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]);
355*1dda9c1aSJeremy L Thompson       // LCOV_EXCL_STOP
356*1dda9c1aSJeremy L Thompson   }
357*1dda9c1aSJeremy L Thompson 
358*1dda9c1aSJeremy L Thompson   // Restore vectors, cover CEED_EVAL_NONE
359*1dda9c1aSJeremy L Thompson   CeedCallBackend(CeedVectorRestoreArrayRead(x_ref, &d_x));
360*1dda9c1aSJeremy L Thompson   CeedCallBackend(CeedVectorRestoreArray(v, &d_v));
361*1dda9c1aSJeremy L Thompson   if (eval_mode == CEED_EVAL_NONE) CeedCallBackend(CeedVectorSetArray(v, CEED_MEM_DEVICE, CEED_COPY_VALUES, (CeedScalar *)d_u));
362*1dda9c1aSJeremy L Thompson   if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u));
363*1dda9c1aSJeremy L Thompson   return CEED_ERROR_SUCCESS;
364*1dda9c1aSJeremy L Thompson }
365*1dda9c1aSJeremy L Thompson 
366*1dda9c1aSJeremy L Thompson //------------------------------------------------------------------------------
3677d8d0e25Snbeams // Destroy basis
3687d8d0e25Snbeams //------------------------------------------------------------------------------
3697d8d0e25Snbeams static int CeedBasisDestroy_Hip_shared(CeedBasis basis) {
3707d8d0e25Snbeams   Ceed                  ceed;
3717d8d0e25Snbeams   CeedBasis_Hip_shared *data;
372b7453713SJeremy L Thompson 
373b7453713SJeremy L Thompson   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
3742b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisGetData(basis, &data));
3752b730f8bSJeremy L Thompson   CeedCallHip(ceed, hipModuleUnload(data->module));
376*1dda9c1aSJeremy L Thompson   if (data->moduleAtPoints) CeedCallHip(ceed, hipModuleUnload(data->moduleAtPoints));
377097cc795SJames Wright   if (data->d_q_weight_1d) CeedCallHip(ceed, hipFree(data->d_q_weight_1d));
3782b730f8bSJeremy L Thompson   CeedCallHip(ceed, hipFree(data->d_interp_1d));
3792b730f8bSJeremy L Thompson   CeedCallHip(ceed, hipFree(data->d_grad_1d));
3802b730f8bSJeremy L Thompson   CeedCallHip(ceed, hipFree(data->d_collo_grad_1d));
381*1dda9c1aSJeremy L Thompson   CeedCallHip(ceed, hipFree(data->d_chebyshev_interp_1d));
3822b730f8bSJeremy L Thompson   CeedCallBackend(CeedFree(&data));
383e15f9bd0SJeremy L Thompson   return CEED_ERROR_SUCCESS;
3847d8d0e25Snbeams }
3857d8d0e25Snbeams 
3867d8d0e25Snbeams //------------------------------------------------------------------------------
3877d8d0e25Snbeams // Create tensor basis
3887d8d0e25Snbeams //------------------------------------------------------------------------------
3892b730f8bSJeremy L Thompson int CeedBasisCreateTensorH1_Hip_shared(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d,
3906574a04fSJeremy L Thompson                                        const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) {
3917d8d0e25Snbeams   Ceed                  ceed;
39222070f95SJeremy L Thompson   char                 *basis_kernel_source;
39322070f95SJeremy L Thompson   const char           *basis_kernel_path;
394b7453713SJeremy L Thompson   CeedInt               num_comp;
395b7453713SJeremy L Thompson   const CeedInt         q_bytes      = Q_1d * sizeof(CeedScalar);
396397164e9SSebastian Grimberg   const CeedInt         interp_bytes = q_bytes * P_1d;
3977d8d0e25Snbeams   CeedBasis_Hip_shared *data;
398b7453713SJeremy L Thompson 
399b7453713SJeremy L Thompson   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
4002b730f8bSJeremy L Thompson   CeedCallBackend(CeedCalloc(1, &data));
4017d8d0e25Snbeams 
4027d8d0e25Snbeams   // Copy basis data to GPU
403097cc795SJames Wright   if (q_weight_1d) {
404b7453713SJeremy L Thompson     CeedCallHip(ceed, hipMalloc((void **)&data->d_q_weight_1d, q_bytes));
405b7453713SJeremy L Thompson     CeedCallHip(ceed, hipMemcpy(data->d_q_weight_1d, q_weight_1d, q_bytes, hipMemcpyHostToDevice));
406097cc795SJames Wright   }
407397164e9SSebastian Grimberg   CeedCallHip(ceed, hipMalloc((void **)&data->d_interp_1d, interp_bytes));
408397164e9SSebastian Grimberg   CeedCallHip(ceed, hipMemcpy(data->d_interp_1d, interp_1d, interp_bytes, hipMemcpyHostToDevice));
409397164e9SSebastian Grimberg   CeedCallHip(ceed, hipMalloc((void **)&data->d_grad_1d, interp_bytes));
410397164e9SSebastian Grimberg   CeedCallHip(ceed, hipMemcpy(data->d_grad_1d, grad_1d, interp_bytes, hipMemcpyHostToDevice));
4117d8d0e25Snbeams 
4127d8d0e25Snbeams   // Compute collocated gradient and copy to GPU
413437930d1SJeremy L Thompson   data->d_collo_grad_1d    = NULL;
4149e201c85SYohann   bool has_collocated_grad = dim == 3 && Q_1d >= P_1d;
415b7453713SJeremy L Thompson 
4169e201c85SYohann   if (has_collocated_grad) {
417437930d1SJeremy L Thompson     CeedScalar *collo_grad_1d;
418b7453713SJeremy L Thompson 
4192b730f8bSJeremy L Thompson     CeedCallBackend(CeedMalloc(Q_1d * Q_1d, &collo_grad_1d));
4202b730f8bSJeremy L Thompson     CeedCallBackend(CeedBasisGetCollocatedGrad(basis, collo_grad_1d));
421b7453713SJeremy L Thompson     CeedCallHip(ceed, hipMalloc((void **)&data->d_collo_grad_1d, q_bytes * Q_1d));
422b7453713SJeremy L Thompson     CeedCallHip(ceed, hipMemcpy(data->d_collo_grad_1d, collo_grad_1d, q_bytes * Q_1d, hipMemcpyHostToDevice));
4232b730f8bSJeremy L Thompson     CeedCallBackend(CeedFree(&collo_grad_1d));
4247d8d0e25Snbeams   }
4257d8d0e25Snbeams 
4269e31c45bSnbeams   // Set number of threads per block for basis kernels
4272b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
4282b730f8bSJeremy L Thompson   CeedCallBackend(ComputeBasisThreadBlockSizes(dim, P_1d, Q_1d, num_comp, data->block_sizes));
4299e31c45bSnbeams 
4309e31c45bSnbeams   // Compile basis kernels
4312b730f8bSJeremy L Thompson   CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/hip/hip-shared-basis-tensor.h", &basis_kernel_path));
43223d4529eSJeremy L Thompson   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source -----\n");
4332b730f8bSJeremy L Thompson   CeedCallBackend(CeedLoadSourceToBuffer(ceed, basis_kernel_path, &basis_kernel_source));
43423d4529eSJeremy L Thompson   CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Basis Kernel Source Complete! -----\n");
435eb7e6cafSJeremy L Thompson   CeedCallBackend(CeedCompile_Hip(ceed, basis_kernel_source, &data->module, 11, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d, "T_1D",
436eb7e6cafSJeremy L Thompson                                   CeedIntMax(Q_1d, P_1d), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp, "BASIS_NUM_NODES", CeedIntPow(P_1d, dim),
437eb7e6cafSJeremy L Thompson                                   "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "BASIS_INTERP_BLOCK_SIZE", data->block_sizes[0], "BASIS_GRAD_BLOCK_SIZE",
4382b730f8bSJeremy L Thompson                                   data->block_sizes[1], "BASIS_WEIGHT_BLOCK_SIZE", data->block_sizes[2], "BASIS_HAS_COLLOCATED_GRAD",
4392b730f8bSJeremy L Thompson                                   has_collocated_grad));
440eb7e6cafSJeremy L Thompson   CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "Interp", &data->Interp));
441eb7e6cafSJeremy L Thompson   CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "InterpTranspose", &data->InterpTranspose));
442eb7e6cafSJeremy L Thompson   CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "Grad", &data->Grad));
443eb7e6cafSJeremy L Thompson   CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "GradTranspose", &data->GradTranspose));
444eb7e6cafSJeremy L Thompson   CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "Weight", &data->Weight));
4452b730f8bSJeremy L Thompson   CeedCallBackend(CeedFree(&basis_kernel_path));
4462b730f8bSJeremy L Thompson   CeedCallBackend(CeedFree(&basis_kernel_source));
4477d8d0e25Snbeams 
4482b730f8bSJeremy L Thompson   CeedCallBackend(CeedBasisSetData(basis, data));
4497d8d0e25Snbeams 
4507d8d0e25Snbeams   // Register backend functions
4512b730f8bSJeremy L Thompson   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyTensor_Hip_shared));
452*1dda9c1aSJeremy L Thompson   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAtPoints", CeedBasisApplyAtPoints_Hip_shared));
4532b730f8bSJeremy L Thompson   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroy_Hip_shared));
454e15f9bd0SJeremy L Thompson   return CEED_ERROR_SUCCESS;
4557d8d0e25Snbeams }
4562a86cc9dSSebastian Grimberg 
4577d8d0e25Snbeams //------------------------------------------------------------------------------
458