1 // Copyright (c) 2017-2024, 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 /// @file 9 /// Internal header for CUDA operator full assembly 10 11 #include <ceed.h> 12 13 #if USE_CEEDSIZE 14 typedef CeedSize IndexType; 15 #else 16 typedef CeedInt IndexType; 17 #endif 18 19 //------------------------------------------------------------------------------ 20 // Matrix assembly kernel 21 //------------------------------------------------------------------------------ 22 extern "C" __launch_bounds__(BLOCK_SIZE) __global__ 23 void LinearAssemble(const CeedInt num_elem, const CeedScalar *B_in, const CeedScalar *B_out, const bool *orients_in, 24 const CeedInt8 *curl_orients_in, const bool *orients_out, const CeedInt8 *curl_orients_out, 25 const CeedScalar *__restrict__ qf_array, CeedScalar *__restrict__ values_array) { 26 extern __shared__ CeedScalar s_CT[]; 27 CeedScalar *s_C = s_CT + NUM_NODES_OUT * NUM_NODES_IN; 28 29 const int l = threadIdx.x; // The output column index of each B^T D B operation 30 // such that we have (Bout^T)_ij D_jk Bin_kl = C_il 31 32 // Strides for final output ordering, determined by the reference (interface) implementation of the symbolic assembly, slowest --> fastest: e, 33 // comp_in, comp_out, node_row, node_col 34 const IndexType comp_out_stride = NUM_NODES_OUT * NUM_NODES_IN; 35 const IndexType comp_in_stride = comp_out_stride * NUM_COMP_OUT; 36 const IndexType e_stride = comp_in_stride * NUM_COMP_IN; 37 38 // Strides for QF array, slowest --> fastest: e_in, comp_in, e_out, comp_out, e, q 39 const IndexType q_e_stride = NUM_QPTS; 40 const IndexType q_comp_out_stride = num_elem * q_e_stride; 41 const IndexType q_eval_mode_out_stride = q_comp_out_stride * NUM_COMP_OUT; 42 const IndexType q_comp_in_stride = q_eval_mode_out_stride * NUM_EVAL_MODES_OUT; 43 const IndexType q_eval_mode_in_stride = q_comp_in_stride * NUM_COMP_IN; 44 45 // Loop over each element (if necessary) 46 for (IndexType e = blockIdx.x * blockDim.z + threadIdx.z; e < num_elem; e += gridDim.x * blockDim.z) { 47 for (IndexType comp_in = 0; comp_in < NUM_COMP_IN; comp_in++) { 48 for (IndexType comp_out = 0; comp_out < NUM_COMP_OUT; comp_out++) { 49 for (IndexType i = threadIdx.y; i < NUM_NODES_OUT; i += BLOCK_SIZE_Y) { 50 CeedScalar result = 0.0; 51 IndexType qf_index_comp = q_comp_in_stride * comp_in + q_comp_out_stride * comp_out + q_e_stride * e; 52 53 for (IndexType e_in = 0; e_in < NUM_EVAL_MODES_IN; e_in++) { 54 IndexType b_in_index = e_in * NUM_QPTS * NUM_NODES_IN; 55 56 for (IndexType e_out = 0; e_out < NUM_EVAL_MODES_OUT; e_out++) { 57 IndexType b_out_index = e_out * NUM_QPTS * NUM_NODES_OUT; 58 IndexType qf_index = qf_index_comp + q_eval_mode_out_stride * e_out + q_eval_mode_in_stride * e_in; 59 60 // Perform the B^T D B operation for this 'chunk' of D (the qf_array) 61 for (IndexType j = 0; j < NUM_QPTS; j++) { 62 result += B_out[b_out_index + j * NUM_NODES_OUT + i] * qf_array[qf_index + j] * B_in[b_in_index + j * NUM_NODES_IN + l]; 63 } 64 } // end of out eval mode 65 } // end of in eval mode 66 if (orients_in) { 67 result *= orients_in[NUM_NODES_IN * e + l] ? -1.0 : 1.0; 68 } 69 if (orients_out) { 70 result *= orients_out[NUM_NODES_OUT * e + i] ? -1.0 : 1.0; 71 } 72 if (!curl_orients_in && !curl_orients_out) { 73 IndexType val_index = e_stride * e + comp_in_stride * comp_in + comp_out_stride * comp_out + NUM_NODES_IN * i + l; 74 75 values_array[val_index] = result; 76 } else if (curl_orients_in) { 77 s_C[NUM_NODES_IN * threadIdx.y + l] = result; 78 __syncthreads(); 79 s_CT[NUM_NODES_IN * i + l] = 80 (l > 0 ? s_C[NUM_NODES_IN * threadIdx.y + l - 1] * curl_orients_in[3 * NUM_NODES_IN * e + 3 * l - 1] : 0.0) + 81 s_C[NUM_NODES_IN * threadIdx.y + l] * curl_orients_in[3 * NUM_NODES_IN * e + 3 * l + 1] + 82 (l < (NUM_NODES_IN - 1) ? s_C[NUM_NODES_IN * threadIdx.y + l + 1] * curl_orients_in[3 * NUM_NODES_IN * e + 3 * l + 3] : 0.0); 83 } else { 84 s_CT[NUM_NODES_IN * i + l] = result; 85 } 86 } // end of loop over element node index, i 87 if (curl_orients_in || curl_orients_out) { 88 // Compute and store the final T^T (B^T D B T) using the fully computed C T product in shared memory 89 if (curl_orients_out) __syncthreads(); 90 for (IndexType i = threadIdx.y; i < NUM_NODES_OUT; i += BLOCK_SIZE_Y) { 91 IndexType val_index = e_stride * e + comp_in_stride * comp_in + comp_out_stride * comp_out + NUM_NODES_IN * i + l; 92 93 if (curl_orients_out) { 94 values_array[val_index] = 95 (i > 0 ? s_CT[NUM_NODES_IN * (i - 1) + l] * curl_orients_out[3 * NUM_NODES_OUT * e + 3 * i - 1] : 0.0) + 96 s_CT[NUM_NODES_IN * i + l] * curl_orients_out[3 * NUM_NODES_OUT * e + 3 * i + 1] + 97 (i < (NUM_NODES_OUT - 1) ? s_CT[NUM_NODES_IN * (i + 1) + l] * curl_orients_out[3 * NUM_NODES_OUT * e + 3 * i + 3] : 0.0); 98 } else { 99 values_array[val_index] = s_CT[NUM_NODES_IN * i + l]; 100 } 101 } 102 } 103 } // end of out component 104 } // end of in component 105 } // end of element loop 106 } 107