From afbd3b0347c24f8679efb15f4c41b4d164639721 Mon Sep 17 00:00:00 2001 From: ddg93 Date: Fri, 10 Jul 2026 11:18:21 +0200 Subject: [PATCH 1/9] Feature: Jacobi preconditioner for GPU execution. CUDA Unified memory allocation of the preconditioning matrix. --- .../linear_algebra/CPreconditioner.hpp | 28 ++++++- Common/include/linear_algebra/CSysMatrix.hpp | 17 +++++ Common/include/linear_algebra/GPUComms.cuh | 13 ++++ .../include/toolboxes/allocation_toolbox.hpp | 45 +++++++++++ Common/src/linear_algebra/CSysMatrix.cpp | 23 +++++- Common/src/linear_algebra/CSysMatrixGPU.cu | 74 ++++++++++++++++++- 6 files changed, 194 insertions(+), 6 deletions(-) diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index e4fc7cf159f..625caaeda9f 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -114,13 +114,37 @@ class CJacobiPreconditioner final : public CPreconditioner { * \param[out] v - CSysVector that is the result of the preconditioning */ inline void operator()(const CSysVector& u, CSysVector& v) const override { - sparse_matrix.ComputeJacobiPreconditioner(u, v, geometry, config); + if (config->GetCUDA()) { +#ifdef HAVE_CUDA + sparse_matrix.GPUComputeJacobiPreconditioner(u, v, geometry, config); +#else + SU2_MPI::Error( + "\nError in launching sparse matrix Preconditioner Function\nENABLE_CUDA is set to YES\nPlease compile with CUDA " + "options enabled in Meson to access GPU Functions", + CURRENT_FUNCTION); +#endif + } else { + sparse_matrix.ComputeJacobiPreconditioner(u, v, geometry, config); + } } /*! * \note Request the associated matrix to build the preconditioner. */ - inline void Build() override { sparse_matrix.BuildJacobiPreconditioner(); } + inline void Build() override { + if (config->GetCUDA()) { +#ifdef HAVE_CUDA + sparse_matrix.GPUBuildJacobiPreconditioner(); +#else + SU2_MPI::Error( + "\nError in building sparse matrix Preconditioner Function\nENABLE_CUDA is set to YES\nPlease compile with CUDA " + "options enabled in Meson to access GPU Functions", + CURRENT_FUNCTION); +#endif + } else { + sparse_matrix.BuildJacobiPreconditioner(); + } + } }; /*! diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 11baf6edf8e..9aa010ee814 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -177,6 +177,8 @@ class CSysMatrix { CCompressedSparsePatternUL levels_ilu; ScalarType* invM; /*!< \brief Inverse of (Jacobi) preconditioner. */ + ScalarType* d_invM; /*!< \brief Inverse of (Jacobi) preconditioner on device. */ + bool invM_is_managed = false; /*!< \brief Boolean that indicates whether GPU supports Unified Memory or not */ /*--- Temporary (hence mutable) working memory used in the Linelet preconditioner, outer vector is for threads ---*/ mutable vector > @@ -924,6 +926,11 @@ class CSysMatrix { */ void BuildJacobiPreconditioner(); + /*! + * \brief Build the Jacobi preconditioner on GPU. + */ + void GPUBuildJacobiPreconditioner(); + /*! * \brief Multiply CSysVector by the preconditioner * \param[in] vec - CSysVector to be multiplied by the preconditioner. @@ -934,6 +941,16 @@ class CSysMatrix { void ComputeJacobiPreconditioner(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const; + /*! + * \brief Multiply CSysVector by the preconditioner + * \param[in] vec - CSysVector to be multiplied by the preconditioner. + * \param[out] prod - Result of the product A*vec. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void GPUComputeJacobiPreconditioner(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, + const CConfig* config) const; + /*! * \brief Build the ILU preconditioner. */ diff --git a/Common/include/linear_algebra/GPUComms.cuh b/Common/include/linear_algebra/GPUComms.cuh index 13854381877..afea91848e0 100644 --- a/Common/include/linear_algebra/GPUComms.cuh +++ b/Common/include/linear_algebra/GPUComms.cuh @@ -51,3 +51,16 @@ inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=t } #define gpuErrChk(ans) { gpuAssert((ans), __FILE__, __LINE__); } + +/*! + * \brief Prefetch a CUDA Unified Memory array to a device asynchronously + * \param[in] ptr, pointer to the memory we want to prefetch. + * \param[in] size in bytes. + * \tparam ZeroInit, initialize memory to 0. + */ +template +inline void gpu_um_prefetch(T* ptr, size_t size, int device) noexcept { +#ifdef HAVE_CUDA + gpuErrChk(cudaMemPrefetchAsync((void*)ptr, size, device)); +#endif +} \ No newline at end of file diff --git a/Common/include/toolboxes/allocation_toolbox.hpp b/Common/include/toolboxes/allocation_toolbox.hpp index 9c357405b75..bdbe5693365 100644 --- a/Common/include/toolboxes/allocation_toolbox.hpp +++ b/Common/include/toolboxes/allocation_toolbox.hpp @@ -98,6 +98,51 @@ inline void aligned_free(T* ptr) noexcept { } // namespace MemoryAllocation namespace GPUMemoryAllocation { +/*! + * \brief Unified Memory support verification + * \return 1 if UM is supported, false otherwise + */ +inline bool UMSupported() noexcept { +#if defined(HAVE_CUDA) + static int managed = -1; + if (managed == -1) { + int device = 0; //one device per process supposed + gpuErrChk(cudaGetDevice(&device)); + gpuErrChk(cudaDeviceGetAttribute(&managed, cudaDevAttrManagedMemory, device)); + } + return managed == 1; +#else + return false; +#endif +} +inline int GetCurrentDevice() noexcept { +#ifdef HAVE_CUDA + static int device = -1; + if (device == -1) gpuErrChk(cudaGetDevice(&device)); + return device; +#else + return -1; +#endif +} + +/*! + * \brief Memory allocation for variables through CUDA Unified Memory: one pointer valid on host and device + * \param[in] size in bytes. + * \tparam ZeroInit, initialize memory to 0. + * \return Pointer to memory, always use gpu_free to deallocate. + */ +template +inline T* gpu_um_alloc(size_t size) noexcept { + void* ptr = nullptr; + +#if defined(HAVE_CUDA) + gpuErrChk(cudaMallocManaged((void**)(&ptr), size)); + if (ZeroInit) std::memset(ptr, 0, size); +#else + return 0; +#endif + return static_cast(ptr); +} /*! * \brief Memory allocation for variables on the GPU. * \param[in] size in bytes. diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index f02d908a059..610b048a229 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -104,7 +104,6 @@ CSysMatrix::~CSysMatrix() { MemoryAllocation::aligned_free(mat.d); MemoryAllocation::aligned_free(mat.l); MemoryAllocation::aligned_free(mat.u); - MemoryAllocation::aligned_free(invM); if (useCuda) { GPUMemoryAllocation::gpu_free(gpu.d); @@ -116,6 +115,15 @@ CSysMatrix::~CSysMatrix() { GPUMemoryAllocation::gpu_free(gpu.col_ind_u); } + if (invM_is_managed) { + GPUMemoryAllocation::gpu_free(invM); + } else { + MemoryAllocation::aligned_free(invM); + if (useCuda) { + GPUMemoryAllocation::gpu_free(d_invM); + } + } + #ifdef USE_MKL mkl_jit_destroy(MatrixMatrixProductJitter); mkl_jit_destroy(MatrixVectorProductJitterBetaZero); @@ -241,7 +249,18 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi allocAndInit(ilu.u, ilu.nnz_u * nVar * nEqn); } - if (diag_needed) allocAndInit(invM, nPointDomain * nVar * nEqn); + if (diag_needed) { + if (useCuda && GPUMemoryAllocation::UMSupported()) { + invM = GPUMemoryAllocation::gpu_um_alloc(nPointDomain * nVar * nEqn * sizeof(ScalarType)); + d_invM = invM; // temporary alias for testing + invM_is_managed = true; + } else { + allocAndInit(invM, nPointDomain * nVar * nEqn); + if (useCuda) { + d_invM = GPUMemoryAllocation::gpu_alloc(nPointDomain * nVar * nEqn * sizeof(ScalarType)); + } + } + } /*--- Thread parallel initialization. ---*/ diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 1ebba30097c..31b099dbef7 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -1,7 +1,7 @@ /*! * \file CSysMatrixGPU.cu * \brief Implementations of Kernels and Functions for Matrix Operations on the GPU - * \author A. Raj + * \author A. Raj, Jesse Li, D. Di Giusto * \version 8.5.0 "Harrier" * * SU2 Project Website: https://su2code.github.io @@ -26,8 +26,32 @@ */ #include "../../include/linear_algebra/CSysMatrix.hpp" +#include "../../include/linear_algebra/CSysMatrix.inl" #include "../../include/linear_algebra/GPUComms.cuh" +/*! + * \brief Matrix-vector product kernel. + */ +template +__global__ void GPUMatrixVectorProductKernel(matrixType *invM, vectorType* vec, vectorType* prod, unsigned long nPointDomain, unsigned long nVar) +{ + + const unsigned long iPoint = blockIdx.x * blockDim.x + threadIdx.x; + if (iPoint >= nPointDomain) return; + + const auto block = &invM[iPoint * nVar * nVar]; + const auto rhs = &vec[iPoint * nVar]; + auto out = &prod[iPoint * nVar]; + + for (auto iVar = 0; iVar < nVar; ++iVar) { + vectorType sum = vectorType(0); + for (auto jVar = 0; jVar < nVar; ++jVar) { + sum += block[iVar * nVar + jVar] * rhs[jVar]; + } + out[iVar] = sum; + } +} + /*! * \brief Block-LDU SpMV kernel: y[iRow] = (L + D + U) * x per block-row. * One CUDA block per block-row; threadIdx.x indexes output variable (0..nVar-1). @@ -96,4 +120,50 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector prod.DtHTransfer(); } -template class CSysMatrix; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits. +template +void CSysMatrix::GPUComputeJacobiPreconditioner(const CSysVector& vec, + CSysVector& prod, CGeometry* geometry, + const CConfig* config) const { + SU2_ZONE_SCOPED + /*--- Apply Jacobi preconditioner, y = D^{-1} * x, the inverse of the diagonal is already known and synced to device ---*/ + + ScalarType* d_vec = vec.GetDevicePointer(); + ScalarType* d_prod = prod.GetDevicePointer(); + + vec.HtDTransfer(); // this is now the entry point of the cuda section so we always want to copy + prod.GPUSetVal(0.0); + + dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE,1,1); + int gridx = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, nPointDomain); + dim3 gridDim(gridx, 1, 1); + + GPUMatrixVectorProductKernel<<>>(d_invM, d_vec, d_prod, nPointDomain, nVar); + gpuErrChk( cudaPeekAtLastError() ); + + prod.DtHTransfer();//forcely copy back prod to host for MPI synch + + /*--- MPI Parallelization ---*/ + CSysMatrixComms::Initiate(prod, geometry, config); + CSysMatrixComms::Complete(prod, geometry, config); + prod.HtDTransfer();//forcely copy back prod to device to continue calculations + +} + +template +void CSysMatrix::GPUBuildJacobiPreconditioner() { + SU2_ZONE_SCOPED + /*--- Build Jacobi preconditioner (M = D), compute and store the inverses of the diagonal blocks. ---*/ + SU2_OMP_FOR_DYN(omp_heavy_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) + InverseDiagonalBlock(iPoint, &(invM[iPoint * nVar * nVar])); + END_SU2_OMP_FOR + + //copy to device or prefetch + if (invM_is_managed) { + gpu_um_prefetch(d_invM, nPointDomain * nVar * nVar * sizeof(ScalarType), GPUMemoryAllocation::GetCurrentDevice()); + } else { + gpuErrChk(cudaMemcpy(d_invM, invM, nPointDomain * nVar * nVar * sizeof(ScalarType), cudaMemcpyHostToDevice)); + } +} + +template class CSysMatrix; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits. \ No newline at end of file From d142c0948764d1494ac882bec34108773617b3b0 Mon Sep 17 00:00:00 2001 From: ddg93 Date: Tue, 14 Jul 2026 14:59:07 +0200 Subject: [PATCH 2/9] run pre-commit to fix the code style/formatting --- Common/include/linear_algebra/CPreconditioner.hpp | 6 ++++-- Common/include/linear_algebra/CSysMatrix.hpp | 12 ++++++------ Common/include/linear_algebra/GPUComms.cuh | 2 +- Common/include/toolboxes/allocation_toolbox.hpp | 2 +- Common/src/linear_algebra/CSysMatrix.cpp | 2 +- Common/src/linear_algebra/CSysMatrixGPU.cu | 4 ++-- 6 files changed, 15 insertions(+), 13 deletions(-) diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index 625caaeda9f..b73a1f8967c 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -119,7 +119,8 @@ class CJacobiPreconditioner final : public CPreconditioner { sparse_matrix.GPUComputeJacobiPreconditioner(u, v, geometry, config); #else SU2_MPI::Error( - "\nError in launching sparse matrix Preconditioner Function\nENABLE_CUDA is set to YES\nPlease compile with CUDA " + "\nError in launching sparse matrix Preconditioner Function\nENABLE_CUDA is set to YES\nPlease compile with " + "CUDA " "options enabled in Meson to access GPU Functions", CURRENT_FUNCTION); #endif @@ -137,7 +138,8 @@ class CJacobiPreconditioner final : public CPreconditioner { sparse_matrix.GPUBuildJacobiPreconditioner(); #else SU2_MPI::Error( - "\nError in building sparse matrix Preconditioner Function\nENABLE_CUDA is set to YES\nPlease compile with CUDA " + "\nError in building sparse matrix Preconditioner Function\nENABLE_CUDA is set to YES\nPlease compile with " + "CUDA " "options enabled in Meson to access GPU Functions", CURRENT_FUNCTION); #endif diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 9aa010ee814..d11ccbdf370 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -176,9 +176,9 @@ class CSysMatrix { /*!< \brief Level structure for alternative shared memory parallelization of ILU. */ CCompressedSparsePatternUL levels_ilu; - ScalarType* invM; /*!< \brief Inverse of (Jacobi) preconditioner. */ - ScalarType* d_invM; /*!< \brief Inverse of (Jacobi) preconditioner on device. */ - bool invM_is_managed = false; /*!< \brief Boolean that indicates whether GPU supports Unified Memory or not */ + ScalarType* invM; /*!< \brief Inverse of (Jacobi) preconditioner. */ + ScalarType* d_invM; /*!< \brief Inverse of (Jacobi) preconditioner on device. */ + bool invM_is_managed = false; /*!< \brief Boolean that indicates whether GPU supports Unified Memory or not */ /*--- Temporary (hence mutable) working memory used in the Linelet preconditioner, outer vector is for threads ---*/ mutable vector > @@ -948,9 +948,9 @@ class CSysMatrix { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void GPUComputeJacobiPreconditioner(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, - const CConfig* config) const; - + void GPUComputeJacobiPreconditioner(const CSysVector& vec, CSysVector& prod, + CGeometry* geometry, const CConfig* config) const; + /*! * \brief Build the ILU preconditioner. */ diff --git a/Common/include/linear_algebra/GPUComms.cuh b/Common/include/linear_algebra/GPUComms.cuh index afea91848e0..1f8f1202135 100644 --- a/Common/include/linear_algebra/GPUComms.cuh +++ b/Common/include/linear_algebra/GPUComms.cuh @@ -63,4 +63,4 @@ inline void gpu_um_prefetch(T* ptr, size_t size, int device) noexcept { #ifdef HAVE_CUDA gpuErrChk(cudaMemPrefetchAsync((void*)ptr, size, device)); #endif -} \ No newline at end of file +} diff --git a/Common/include/toolboxes/allocation_toolbox.hpp b/Common/include/toolboxes/allocation_toolbox.hpp index bdbe5693365..35eeda15353 100644 --- a/Common/include/toolboxes/allocation_toolbox.hpp +++ b/Common/include/toolboxes/allocation_toolbox.hpp @@ -106,7 +106,7 @@ inline bool UMSupported() noexcept { #if defined(HAVE_CUDA) static int managed = -1; if (managed == -1) { - int device = 0; //one device per process supposed + int device = 0; // one device per process supposed gpuErrChk(cudaGetDevice(&device)); gpuErrChk(cudaDeviceGetAttribute(&managed, cudaDevAttrManagedMemory, device)); } diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 610b048a229..a865859633d 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -252,7 +252,7 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi if (diag_needed) { if (useCuda && GPUMemoryAllocation::UMSupported()) { invM = GPUMemoryAllocation::gpu_um_alloc(nPointDomain * nVar * nEqn * sizeof(ScalarType)); - d_invM = invM; // temporary alias for testing + d_invM = invM; // temporary alias for testing invM_is_managed = true; } else { allocAndInit(invM, nPointDomain * nVar * nEqn); diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 31b099dbef7..a750659e3d5 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -35,7 +35,7 @@ template __global__ void GPUMatrixVectorProductKernel(matrixType *invM, vectorType* vec, vectorType* prod, unsigned long nPointDomain, unsigned long nVar) { - + const unsigned long iPoint = blockIdx.x * blockDim.x + threadIdx.x; if (iPoint >= nPointDomain) return; @@ -166,4 +166,4 @@ void CSysMatrix::GPUBuildJacobiPreconditioner() { } } -template class CSysMatrix; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits. \ No newline at end of file +template class CSysMatrix; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits. From 83efe791b85d5f70d8fd295247c2d12f116895e9 Mon Sep 17 00:00:00 2001 From: ddg93 Date: Sat, 18 Jul 2026 15:24:33 +0200 Subject: [PATCH 3/9] implemented multidot on GPU with custom cuda kernel. Validated on 1 test-case --- Common/include/linear_algebra/CSysVector.hpp | 13 ++ Common/src/linear_algebra/CSysVector.cpp | 7 ++ Common/src/linear_algebra/CSysVectorGPU.cu | 123 +++++++++++++++++++ 3 files changed, 143 insertions(+) diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 1498b549bbb..ad2b609687f 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -384,6 +384,19 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> static const su2matrix& multiDot(const std::vector& V, size_t i0, size_t n, const std::vector& W, size_t m); + /*! + * \brief Computes the product of V^T W on the GPU, where V and W are tall matrices stored as vectors of CSysVector. + * \param[in] V - Tall matrix. + * \param[in] i0 - First column of V to consider. + * \param[in] n - Number of columns to consider from V starting at i0. + * \param[in] W - Tall matrix. + * \param[in] m - Number of columns to consider from W. + * \return n by m matrix with the result of the product. + */ + static const su2matrix& multiDotGPU(const std::vector>& V, const size_t i0, + const size_t n, const std::vector>& W, + const size_t m); + /*! * \brief Squared L2 norm of the vector (via dot with self). * \return Squared L2 norm. diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp index f7df34633a0..76a58712290 100644 --- a/Common/src/linear_algebra/CSysVector.cpp +++ b/Common/src/linear_algebra/CSysVector.cpp @@ -73,6 +73,13 @@ const su2matrix& CSysVector::multiDot(const std::vector< const std::vector>& W, const size_t m) { SU2_ZONE_SCOPED + +#ifdef HAVE_CUDA + if (V[i0].GetDevicePointer() != nullptr) { + return multiDotGPU(V, i0, n, W, m); + } +#endif + static constexpr size_t BLOCK_SIZE = 1024; static su2matrix shared; diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 94ec17bb88f..2081118979c 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -46,4 +46,127 @@ void CSysVector::GPUSetVal(ScalarType val, bool trigger) const if(trigger) gpuErrChk(cudaMemset((void*)(d_vec_val), val, (sizeof(ScalarType)*nElm))); } +/*! + * \brief multi vector product CUDA kernel one line of blocks per pair V[i0+i],W[j]; + * Configurable multiple blocks reducing over the size of the vectors + */ +template +__global__ void GPUmultiDot(const ScalarType* const* __restrict__ d_V, const size_t n, + const ScalarType* const* __restrict__ d_W, const size_t m, const size_t size, + ScalarType* __restrict__ d_local) +{ + // Map each x,y block to the specific (i,j) dot product + const size_t pair_idx = blockIdx.y; + if (pair_idx >= n * m) return; + + const size_t i = pair_idx / m; + const size_t j = pair_idx % m; + + //get the corresponding vectors + const ScalarType* __restrict__ vi = d_V[i]; + const ScalarType* __restrict__ wj = d_W[j]; + + // grid strided loop over the vector elements + ScalarType local_sum = 0.0; + const size_t tid = blockIdx.x * blockDim.x + threadIdx.x; + const size_t stride = gridDim.x * blockDim.x; + + for (size_t k = tid; k < size; k += stride) + { + local_sum += vi[k] * wj[k]; + } + + // shared memory reduction within the block + extern __shared__ char shared_mem[]; + ScalarType* sdata = reinterpret_cast(shared_mem); + + sdata[threadIdx.x] = local_sum; + __syncthreads(); + + // parallel reduction on the block + for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) + { + if (threadIdx.x < s) + { + sdata[threadIdx.x] += sdata[threadIdx.x + s]; + } + __syncthreads(); + } + + // atomic add of each block partial sum to the output matrix, operated by thread 0 of each block + if (threadIdx.x == 0) + { + atomicAdd(&d_local[i * m + j], sdata[0]); + } +} + +template +const su2matrix& CSysVector::multiDotGPU(const std::vector>& V, + const size_t i0, const size_t n, + const std::vector>& W, + const size_t m) { + SU2_ZONE_SCOPED + + static su2matrix shared; + if (n == 0 || m == 0) return shared; + + const size_t size = V[0].nElmDomain; + + // ensure all vectors are synced on device and get the pointers + std::vector h_V_ptrs(n), h_W_ptrs(m); + for (size_t i = 0; i < n; ++i){ + V[i0 + i].HtDTransfer(); + h_V_ptrs[i] = V[i0 + i].GetDevicePointer(); + } + for (size_t j = 0; j < m; ++j){ + W[j].HtDTransfer(); + h_W_ptrs[j] = W[j].GetDevicePointer(); + } + + //copy the pointers to the device arrays of pointers + const ScalarType** d_V_ptrs; + const ScalarType** d_W_ptrs; + gpuErrChk(cudaMalloc(&d_V_ptrs, n * sizeof(ScalarType*))); + gpuErrChk(cudaMalloc(&d_W_ptrs, m * sizeof(ScalarType*))); + gpuErrChk(cudaMemcpy(d_V_ptrs, h_V_ptrs.data(), n * sizeof(ScalarType*), cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpy(d_W_ptrs, h_W_ptrs.data(), m * sizeof(ScalarType*), cudaMemcpyHostToDevice)); + + // allocate result buffer, zero it + ScalarType* d_local; + gpuErrChk(cudaMalloc(&d_local, n * m * sizeof(ScalarType))); + gpuErrChk(cudaMemset(d_local, 0, n * m * sizeof(ScalarType))); + + // launch + int threads = 256; + int numBlocksPerPair = std::min((size + threads - 1) / threads, size_t(1024)); + dim3 grid(numBlocksPerPair, n * m); + GPUmultiDot<<>>(d_V_ptrs, n, d_W_ptrs, m, size, d_local); + gpuErrChk(cudaGetLastError()); + + // copy result to host, MPI reduce + su2matrix local(n,m); + gpuErrChk(cudaMemcpy(local.data(), d_local, n * m * sizeof(ScalarType), cudaMemcpyDeviceToHost)); + + /*--- Single AllReduce of the result, only the master thread communicates. ---*/ + SU2_OMP_MASTER { + shared.resize(n, m); + + const auto mpi_type = (sizeof(ScalarType) < sizeof(double)) ? MPI_FLOAT : MPI_DOUBLE; + SelectMPIWrapper::W::Allreduce(local.data(), shared.data(), n * m, mpi_type, MPI_SUM, + SU2_MPI::GetComm()); + } + END_SU2_OMP_MASTER + + /*--- All threads have the same view of the result. ---*/ + SU2_OMP_BARRIER + + + // cleanup (or cache these allocations) + gpuErrChk(cudaFree(d_local)); + gpuErrChk(cudaFree(d_V_ptrs)); + gpuErrChk(cudaFree(d_W_ptrs)); + + return shared; +} + template class CSysVector; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits. From 60419c624b2e8570495de7dc2a8899f0a3a0102f Mon Sep 17 00:00:00 2001 From: ddg93 Date: Sun, 19 Jul 2026 22:17:40 +0200 Subject: [PATCH 4/9] implemented linear combination on GPU with custom cuda kernel. Validated on 1 test-case. --- Common/include/linear_algebra/CSysVector.hpp | 11 +++++ Common/src/linear_algebra/CSysSolve.cpp | 12 +++++ Common/src/linear_algebra/CSysVectorGPU.cu | 52 ++++++++++++++++++++ 3 files changed, 75 insertions(+) diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index ad2b609687f..5abcf874de7 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -397,6 +397,17 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> const size_t n, const std::vector>& W, const size_t m); + /*! + * \brief Computes v = vs * ws or v += vs * ws with unrolling of up to 4 iterations on the GPU + * \param[in] n - number of vectors to consider + * \param[in] ws - array of scalar weights corresponding to the device pointers to vectors + * \param[in] vs_ptrs - array of device pointers + * \param[in] v - target vector + * \param[in] inc - If true, adds results to target vector. If false, overwrites + */ + static void LinearCombinationGPU(const unsigned long n, const std::vector>& vs, const ScalarType* ws, + CSysVector& v, bool inc=false); + /*! * \brief Squared L2 norm of the vector (via dot with self). * \return Squared L2 norm. diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index aae5d9ce707..f5e0319d946 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -33,6 +33,7 @@ #include "../../include/linear_algebra/CSysMatrix.hpp" #include "../../include/linear_algebra/CMatrixVectorProduct.hpp" #include "../../include/linear_algebra/CPreconditioner.hpp" +#include "../../include/linear_algebra/CSysVector.hpp" SU2_IGNORE_WARNING("-Wmaybe-uninitialized") #include "Eigen/Eigenvalues" @@ -107,6 +108,17 @@ void LinearCombinationImpl(const unsigned long n, const Vectors& vs, const Weigh template void LinearCombinationImpl(const unsigned long n, const std::vector>& vs, const Weights& ws, CSysVector& v, bool inc = false) { +#ifdef HAVE_CUDA + if (v.GetDevicePointer() != nullptr) { + std::vector ws_host(n); // collect weights into simple host array + for (unsigned long i = 0; i < n; ++i) { + ws_host[i] = static_cast(ws(i)); + } + CSysVector::LinearCombinationGPU(n, vs, ws_host.data(), v, inc); + return; + } +#endif + LinearCombinationImpl( n, [&vs](auto i) -> auto& { return vs[i]; }, ws, v, inc); } diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 2081118979c..303f45897df 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -169,4 +169,56 @@ const su2matrix& CSysVector::multiDotGPU(const std::vect return shared; } +template +struct WeightedVecs { + const ScalarType* ptrs[N]; + ScalarType weights[N]; +}; + +template +__global__ void LinearCombinationKernel(ScalarType* __restrict__ v, WeightedVecs wv, + int n, unsigned long nElm, bool inc) +{ + const unsigned long k = blockIdx.x * blockDim.x + threadIdx.x; + if (k >= nElm) return; + + //handle overwriting or combination with existing + ScalarType result = inc ? v[k] : ScalarType(0); + + #pragma unroll + for (int i = 0; i < N; ++i) // N is known at compile time (4), this unrolls to if (i < n) result += weight[i] * vector[i][k]; i<4 + if (i < n) result += wv.weights[i] * wv.ptrs[i][k]; + v[k] = result; +} + +template +void CSysVector::LinearCombinationGPU(const unsigned long n, const std::vector>& vs, const ScalarType* ws, + CSysVector& v, bool inc) +{ + + const unsigned long nElm = v.nElmDomain; + constexpr unsigned threads = 256; + const unsigned blocks = (nElm + threads - 1) / threads; + + // ensure v is on device before first kernel + v.HtDTransfer(); + ScalarType* d_v = v.GetDevicePointer(); + + for (unsigned long i = 0; i < n; i += 4) { + const int rem = static_cast(std::min(n - i, 4ul)); + //prepare vectors pointers and corresponding weights, passing them by value + WeightedVecs vs_ws = {}; + for (int j = 0; j < rem; ++j) { + vs_ws.ptrs[j] = vs[i + j].GetDevicePointer(); // already on device from multiDot + vs_ws.weights[j] = ws[i + j]; // plain array indexing, not ws(k) + } + //calculate the linear combination on GPU, handle more than 4 vectors through inc || i > 0 + LinearCombinationKernel<<>>(d_v, vs_ws, rem, nElm, inc || i > 0); + gpuErrChk(cudaPeekAtLastError()); + } + + // bring result back to host + v.DtHTransfer(); +} + template class CSysVector; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits. From e070ec7b6a1b559801a875a60765a0d35a54b8a5 Mon Sep 17 00:00:00 2001 From: ddg93 Date: Mon, 20 Jul 2026 00:34:51 +0200 Subject: [PATCH 5/9] implemented CUDA Unified Memory for the CsysVector and cleaned CUDA logic from device pointers. Validated on test-case turb_SA_RAE2822 comparing residuals over 100 iterations. --- Common/include/linear_algebra/CSysVector.hpp | 13 ++++-- Common/src/linear_algebra/CSysMatrixGPU.cu | 17 +++---- Common/src/linear_algebra/CSysSolve.cpp | 12 ++--- Common/src/linear_algebra/CSysVector.cpp | 34 +++++++++++--- Common/src/linear_algebra/CSysVectorGPU.cu | 48 +++++++++----------- 5 files changed, 71 insertions(+), 53 deletions(-) diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 5abcf874de7..cbddd51339d 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -77,6 +77,8 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> unsigned long nVar = 1; /*!< \brief Number of elements in a block. */ ScalarType* d_vec_val = nullptr; /*!< \brief Device Pointer to store the vector values on the GPU. */ + bool vec_is_managed = false; /*!< \brief Boolean that indicates whether GPU supports Unified Memory or not */ + bool useCuda = false; /*!< \brief Whether CUDA is enabled. */ #ifdef HAVE_OMP mutable std::unique_ptr @@ -245,6 +247,11 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ inline ScalarType* GetDevicePointer() const { return d_vec_val; } + /*! + * \brief return pointer that points to the CSysVector values in CPU memory + */ + ScalarType* data() const { return vec_val; } + /*! * \brief return the number of local elements in the CSysVector */ @@ -403,10 +410,10 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] ws - array of scalar weights corresponding to the device pointers to vectors * \param[in] vs_ptrs - array of device pointers * \param[in] v - target vector - * \param[in] inc - If true, adds results to target vector. If false, overwrites + * \param[in] inc - If true, adds results to target vector. If false, overwrites */ - static void LinearCombinationGPU(const unsigned long n, const std::vector>& vs, const ScalarType* ws, - CSysVector& v, bool inc=false); + static void LinearCombinationGPU(const unsigned long n, const std::vector>& vs, + const ScalarType* ws, CSysVector& v, bool inc = false); /*! * \brief Squared L2 norm of the vector (via dot with self). diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index a750659e3d5..29c97002357 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -28,6 +28,7 @@ #include "../../include/linear_algebra/CSysMatrix.hpp" #include "../../include/linear_algebra/CSysMatrix.inl" #include "../../include/linear_algebra/GPUComms.cuh" +#include "../../include/linear_algebra/CSysVector.hpp" /*! * \brief Matrix-vector product kernel. @@ -108,7 +109,7 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector ScalarType* d_vec = vec.GetDevicePointer(); ScalarType* d_prod = prod.GetDevicePointer(); - vec.HtDTransfer(); + //vec.HtDTransfer(); dim3 blockDim(static_cast(nVar), 1, 1); dim3 gridDim(static_cast(nPointDomain), 1, 1); @@ -117,7 +118,7 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector gpu.row_ptr_u, gpu.col_ind_u, gpu.u, d_vec, d_prod); gpuErrChk(cudaGetLastError()); - prod.DtHTransfer(); + //prod.DtHTransfer(); } template @@ -127,25 +128,21 @@ void CSysMatrix::GPUComputeJacobiPreconditioner(const CSysVector>>(d_invM, d_vec, d_prod, nPointDomain, nVar); + GPUMatrixVectorProductKernel<<>>(d_invM, vec.data(), prod.data(), nPointDomain, nVar); gpuErrChk( cudaPeekAtLastError() ); - - prod.DtHTransfer();//forcely copy back prod to host for MPI synch + gpuErrChk(cudaDeviceSynchronize()); /*--- MPI Parallelization ---*/ CSysMatrixComms::Initiate(prod, geometry, config); CSysMatrixComms::Complete(prod, geometry, config); - prod.HtDTransfer();//forcely copy back prod to device to continue calculations + + gpuErrChk(cudaDeviceSynchronize()); } diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index f5e0319d946..14fa8d7b6de 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -109,14 +109,12 @@ template void LinearCombinationImpl(const unsigned long n, const std::vector>& vs, const Weights& ws, CSysVector& v, bool inc = false) { #ifdef HAVE_CUDA - if (v.GetDevicePointer() != nullptr) { - std::vector ws_host(n); // collect weights into simple host array - for (unsigned long i = 0; i < n; ++i) { - ws_host[i] = static_cast(ws(i)); - } - CSysVector::LinearCombinationGPU(n, vs, ws_host.data(), v, inc); - return; + std::vector ws_host(n); // collect weights into simple host array + for (unsigned long i = 0; i < n; ++i) { + ws_host[i] = static_cast(ws(i)); } + CSysVector::LinearCombinationGPU(n, vs, ws_host.data(), v, inc); + return; #endif LinearCombinationImpl( diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp index 76a58712290..09d3740fb0c 100644 --- a/Common/src/linear_algebra/CSysVector.cpp +++ b/Common/src/linear_algebra/CSysVector.cpp @@ -50,9 +50,21 @@ void CSysVector::Initialize(unsigned long numBlk, unsigned long numB omp_chunk_size = computeStaticChunkSize(nElm, omp_get_max_threads(), OMP_MAX_SIZE); - if (vec_val == nullptr) vec_val = MemoryAllocation::aligned_alloc(64, nElm * sizeof(ScalarType)); + // if (vec_val == nullptr) vec_val = MemoryAllocation::aligned_alloc(64, nElm * sizeof(ScalarType)); +#ifdef HAVE_CUDA + useCuda = true; +#endif - d_vec_val = GPUMemoryAllocation::gpu_alloc(nElm * sizeof(ScalarType)); + if (useCuda && GPUMemoryAllocation::UMSupported()) { + vec_val = GPUMemoryAllocation::gpu_um_alloc(nElm * sizeof(ScalarType)); + d_vec_val = vec_val; // temporary alias for testing + vec_is_managed = true; + } else { + vec_val = MemoryAllocation::aligned_alloc(64, nElm * sizeof(ScalarType)); + if (useCuda) { + d_vec_val = GPUMemoryAllocation::gpu_alloc(nElm * sizeof(ScalarType)); + } + } #ifdef HAVE_OMP dot_scratch.reset(new ScalarType[omp_get_max_threads()]); @@ -75,9 +87,7 @@ const su2matrix& CSysVector::multiDot(const std::vector< SU2_ZONE_SCOPED #ifdef HAVE_CUDA - if (V[i0].GetDevicePointer() != nullptr) { - return multiDotGPU(V, i0, n, W, m); - } + return multiDotGPU(V, i0, n, W, m); #endif static constexpr size_t BLOCK_SIZE = 1024; @@ -144,9 +154,19 @@ CSysVector::~CSysVector() { if constexpr (!std::is_trivial_v) { for (auto i = 0ul; i < nElm; i++) vec_val[i].~ScalarType(); } - MemoryAllocation::aligned_free(vec_val); + // MemoryAllocation::aligned_free(vec_val); - GPUMemoryAllocation::gpu_free(d_vec_val); + // GPUMemoryAllocation::gpu_free(d_vec_val); + + if (useCuda && GPUMemoryAllocation::UMSupported()) { + GPUMemoryAllocation::gpu_free(vec_val); + d_vec_val = nullptr; + } else { + MemoryAllocation::aligned_free(vec_val); + if (useCuda) { + GPUMemoryAllocation::gpu_free(d_vec_val); + } + } } /*--- Explicit instantiations ---*/ diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 303f45897df..ba590daa277 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -112,18 +112,17 @@ const su2matrix& CSysVector::multiDotGPU(const std::vect const size_t size = V[0].nElmDomain; - // ensure all vectors are synced on device and get the pointers + // get all the device pointers std::vector h_V_ptrs(n), h_W_ptrs(m); for (size_t i = 0; i < n; ++i){ - V[i0 + i].HtDTransfer(); - h_V_ptrs[i] = V[i0 + i].GetDevicePointer(); + h_V_ptrs[i] = V[i0 + i].data(); } for (size_t j = 0; j < m; ++j){ - W[j].HtDTransfer(); - h_W_ptrs[j] = W[j].GetDevicePointer(); + //gpuErrChk(cudaMemAdvise(vec_val, nElm * sizeof(ScalarType), cudaMemAdviseSetReadMostly, device_id)); //if read only, could be good + h_W_ptrs[j] = W[j].data(); } - //copy the pointers to the device arrays of pointers + // copy the pointers to the device arrays of pointers const ScalarType** d_V_ptrs; const ScalarType** d_W_ptrs; gpuErrChk(cudaMalloc(&d_V_ptrs, n * sizeof(ScalarType*))); @@ -136,14 +135,14 @@ const su2matrix& CSysVector::multiDotGPU(const std::vect gpuErrChk(cudaMalloc(&d_local, n * m * sizeof(ScalarType))); gpuErrChk(cudaMemset(d_local, 0, n * m * sizeof(ScalarType))); - // launch - int threads = 256; - int numBlocksPerPair = std::min((size + threads - 1) / threads, size_t(1024)); - dim3 grid(numBlocksPerPair, n * m); - GPUmultiDot<<>>(d_V_ptrs, n, d_W_ptrs, m, size, d_local); + dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE,1,1); + int numBlocksPerPair = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, size); + dim3 gridDim(numBlocksPerPair, n * m, 1); + + GPUmultiDot<<>>(d_V_ptrs, n, d_W_ptrs, m, size, d_local); gpuErrChk(cudaGetLastError()); - // copy result to host, MPI reduce + // copy result to host for MPI reduce su2matrix local(n,m); gpuErrChk(cudaMemcpy(local.data(), d_local, n * m * sizeof(ScalarType), cudaMemcpyDeviceToHost)); @@ -161,7 +160,7 @@ const su2matrix& CSysVector::multiDotGPU(const std::vect SU2_OMP_BARRIER - // cleanup (or cache these allocations) + // clean allocations gpuErrChk(cudaFree(d_local)); gpuErrChk(cudaFree(d_V_ptrs)); gpuErrChk(cudaFree(d_W_ptrs)); @@ -176,7 +175,7 @@ struct WeightedVecs { }; template -__global__ void LinearCombinationKernel(ScalarType* __restrict__ v, WeightedVecs wv, +__global__ void LinearCombinationKernel(ScalarType* __restrict__ v, WeightedVecs wv, int n, unsigned long nElm, bool inc) { const unsigned long k = blockIdx.x * blockDim.x + threadIdx.x; @@ -186,7 +185,7 @@ __global__ void LinearCombinationKernel(ScalarType* __restrict__ v, WeightedVecs ScalarType result = inc ? v[k] : ScalarType(0); #pragma unroll - for (int i = 0; i < N; ++i) // N is known at compile time (4), this unrolls to if (i < n) result += weight[i] * vector[i][k]; i<4 + for (int i = 0; i < N; ++i) // N is known at compile time (4), this unrolls to: if (i < n) result += weight[i] * vector[i][k]; i<4 if (i < n) result += wv.weights[i] * wv.ptrs[i][k]; v[k] = result; } @@ -197,28 +196,25 @@ void CSysVector::LinearCombinationGPU(const unsigned long n, const s { const unsigned long nElm = v.nElmDomain; - constexpr unsigned threads = 256; - const unsigned blocks = (nElm + threads - 1) / threads; - - // ensure v is on device before first kernel - v.HtDTransfer(); - ScalarType* d_v = v.GetDevicePointer(); + dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE,1,1); + int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, nElm); + dim3 gridDim(numBlocks, 1, 1); for (unsigned long i = 0; i < n; i += 4) { const int rem = static_cast(std::min(n - i, 4ul)); //prepare vectors pointers and corresponding weights, passing them by value WeightedVecs vs_ws = {}; for (int j = 0; j < rem; ++j) { - vs_ws.ptrs[j] = vs[i + j].GetDevicePointer(); // already on device from multiDot - vs_ws.weights[j] = ws[i + j]; // plain array indexing, not ws(k) + vs_ws.ptrs[j] = vs[i + j].data(); // already on device from multiDot + vs_ws.weights[j] = ws[i + j]; // plain array indexing, not ws(k) } //calculate the linear combination on GPU, handle more than 4 vectors through inc || i > 0 - LinearCombinationKernel<<>>(d_v, vs_ws, rem, nElm, inc || i > 0); + LinearCombinationKernel<<>>(v.data(), vs_ws, rem, nElm, inc || i > 0); gpuErrChk(cudaPeekAtLastError()); } - // bring result back to host - v.DtHTransfer(); + gpuErrChk(cudaDeviceSynchronize()); // this is now the exit point of the cuda section so we want to synchronize + } template class CSysVector; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits. From 7ca84162303e17a87dce22db93c6337189cef4b9 Mon Sep 17 00:00:00 2001 From: ddg93 Date: Sun, 26 Jul 2026 19:59:22 +0200 Subject: [PATCH 6/9] Implemented vector-scalar operators with templates, implemented dot product through custom CUDA kernel for norm of GPU, validated on rae2822 test-case comparing residuals between CPU and GPU --- Common/include/linear_algebra/CSysVector.hpp | 56 ++++++++-- Common/src/linear_algebra/CSysVectorGPU.cu | 111 ++++++++++++++++++- 2 files changed, 157 insertions(+), 10 deletions(-) diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index cbddd51339d..28d3121f2db 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -224,6 +224,18 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> END_CSYSVEC_PARFOR } + /*! + * \brief enum listing the possible scalar operators on GPU + */ + enum class GPUScalarOp { SET, ADD, SUB, MUL, DIV }; + + /*! + * \brief method to launch the generic Unary Operation Kernel and apply given functors on GPU + * \param[in] op - Unitary operation + * \param[in] val - scalar value + */ + void GPUUnaryOperation(GPUScalarOp op, ScalarType val); + /*! * \brief Performs the memory copy from host to device. * \param[in] trigger - boolean value that decides whether to conduct the transfer or not. True by default. @@ -318,7 +330,21 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \brief Compound assignement operations with scalars and expressions. * \param[in] val/expr - Scalar value or expression. */ -#define MAKE_COMPOUND(OP) \ +#ifdef HAVE_CUDA +#define MAKE_COMPOUND(OP, TAG) \ + CSysVector& operator OP(ScalarType val) { \ + GPUUnaryOperation(GPUScalarOp::TAG, val); \ + return *this; \ + } \ + template \ + CSysVector& operator OP(const VecExpr::CVecExpr& expr) { \ + CSYSVEC_PARFOR \ + for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP expr.derived()[i]; \ + END_CSYSVEC_PARFOR \ + return *this; \ + } +#else +#define MAKE_COMPOUND(OP, TAG) \ CSysVector& operator OP(ScalarType val) { \ CSYSVEC_PARFOR \ for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val; \ @@ -332,11 +358,13 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> END_CSYSVEC_PARFOR \ return *this; \ } - MAKE_COMPOUND(=) - MAKE_COMPOUND(+=) - MAKE_COMPOUND(-=) - MAKE_COMPOUND(*=) - MAKE_COMPOUND(/=) +#endif + + MAKE_COMPOUND(=, SET) + MAKE_COMPOUND(+=, ADD) + MAKE_COMPOUND(-=, SUB) + MAKE_COMPOUND(*=, MUL) + MAKE_COMPOUND(/=, DIV) #undef MAKE_COMPOUND /*! @@ -353,9 +381,12 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> ScalarType dot(const VecExpr::CVecExpr& expr) const { /*--- All threads get the same "view" of the vectors. ---*/ SU2_OMP_BARRIER + ScalarType sum = 0.0; +#ifdef HAVE_CUDA + sum = GPUDot(static_cast(expr.derived())); //assuming expr is a vector +#else /*--- Local dot product for each thread. ---*/ - ScalarType sum = 0.0; CSYSVEC_PARFOR for (auto i = 0ul; i < nElmDomain; ++i) { @@ -364,11 +395,13 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> END_CSYSVEC_PARFOR dot_scratch[omp_get_thread_num()] = sum; +#endif BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { +#ifndef HAVE_CUDA /*--- Reduce over all threads in an ordered way to ensure a deterministic result. ---*/ for (int i = 1; i < omp_get_num_threads(); ++i) sum += dot_scratch[i]; - +#endif /*--- Reduce across all mpi ranks, only the master thread communicates. ---*/ const auto mpi_type = (sizeof(ScalarType) < sizeof(double)) ? MPI_FLOAT : MPI_DOUBLE; SelectMPIWrapper::W::Allreduce(&sum, &dot_scratch[0], 1, mpi_type, MPI_SUM, SU2_MPI::GetComm()); @@ -379,6 +412,13 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> return dot_scratch[0]; } + /*! + * \brief Dot product between "this" and another vector. + * \param[in] other - the other CSysVector. + * \return Result of dot product + */ + ScalarType GPUDot(const CSysVector& other) const; + /*! * \brief Computes the product of V^T W efficiencly, where V and W are tall matrices stored as vectors of CSysVector. * \param[in] V - Tall matrix. diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index ba590daa277..103b41a8851 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -1,7 +1,7 @@ /*! * \file CSysVectorGPU.cu * \brief Implementations of Kernels and Functions for Vector Operations on the GPU - * \author A. Raj + * \author A. Raj, D. Di giusto * \version 8.5.0 "Harrier" * * SU2 Project Website: https://su2code.github.io @@ -28,6 +28,113 @@ #include "../../include/linear_algebra/CSysVector.hpp" #include "../../include/linear_algebra/GPUComms.cuh" +/*! + * \brief block-level reduction of elementwise products, accumulated into a single device scalar. + */ +template +__global__ void GPUDotKernel(const ScalarType* __restrict__ a, const ScalarType* __restrict__ b, + unsigned long n, ScalarType* __restrict__ result) { + //shared memory + extern __shared__ unsigned char smem_raw[]; + ScalarType* sdata = reinterpret_cast(smem_raw); + + //local and global thread indexes for access and block reduction + const unsigned long tid = threadIdx.x; + unsigned long idx = blockIdx.x * blockDim.x + tid; + const unsigned long stride = blockDim.x * gridDim.x; + + //thread reduction + ScalarType local = ScalarType(0); + for (; idx < n; idx += stride) local += a[idx] * b[idx]; + + sdata[tid] = local; + __syncthreads(); + + //block reduction + for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) { + if (tid < s) sdata[tid] += sdata[tid + s]; + __syncthreads(); + } + + //final atomic add per block + if (tid == 0) atomicAdd(result, sdata[0]); +} + +template +ScalarType CSysVector::GPUDot(const CSysVector& other) const { + + dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); + int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain); + dim3 gridDim(numBlocks, 1, 1); + + // allocate and zero the result scalar + ScalarType* d_dot_result; + gpuErrChk(cudaMalloc(&d_dot_result, sizeof(ScalarType))); + gpuErrChk(cudaMemset(d_dot_result, 0, sizeof(ScalarType))); + + const size_t sharedBytes = KernelParameters::MVP_BLOCK_SIZE * sizeof(ScalarType); + GPUDotKernel<<>>(this->d_vec_val, other.d_vec_val, this->nElmDomain, d_dot_result); + gpuErrChk(cudaPeekAtLastError()); + + ScalarType result; + gpuErrChk(cudaMemcpy(&result, d_dot_result, sizeof(ScalarType), cudaMemcpyDeviceToHost)); + gpuErrChk(cudaFree(d_dot_result)); + + + + return result; +} + +/*! + * \brief namespace defining the scalar operators for GPU + */ +namespace { +template struct OpSetScalar { T val; __device__ __forceinline__ T operator()(T x) const { return val; } }; +template struct OpAddScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x + val; } }; +template struct OpSubScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x - val; } }; +template struct OpMulScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x * val; } }; +template struct OpDivScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x / val; } }; +} // namespace + +/*! + * \brief generic Unary Operation Kernel to apply given functors on GPU + */ +template +__global__ void GPUUnaryOperationKernel(ScalarType* __restrict__ vec, unsigned long n, Operator op) { + const unsigned long idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < n) vec[idx] = op(vec[idx]); +} + +template +void CSysVector::GPUUnaryOperation(GPUScalarOp op, ScalarType val) { + + dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); + int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain); + dim3 gridDim(numBlocks, 1, 1); + + switch (op) { + case GPUScalarOp::SET: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSetScalar{val}); + break; + case GPUScalarOp::ADD: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpAddScalar{val}); + break; + case GPUScalarOp::SUB: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSubScalar{val}); + break; + case GPUScalarOp::MUL: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpMulScalar{val}); + break; + case GPUScalarOp::DIV: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpDivScalar{val}); + break; + } + gpuErrChk(cudaPeekAtLastError()); + gpuErrChk(cudaDeviceSynchronize()); + + +} + template void CSysVector::HtDTransfer(bool trigger) const { @@ -213,7 +320,7 @@ void CSysVector::LinearCombinationGPU(const unsigned long n, const s gpuErrChk(cudaPeekAtLastError()); } - gpuErrChk(cudaDeviceSynchronize()); // this is now the exit point of the cuda section so we want to synchronize + gpuErrChk(cudaDeviceSynchronize()); } From 781d19d6f62895dba01c5a891437e2e063527c17 Mon Sep 17 00:00:00 2001 From: ddg93 Date: Sun, 26 Jul 2026 22:12:57 +0200 Subject: [PATCH 7/9] implemented vector-vector operators for GPU explicitly, not the negative assignment, validated on rae2822 --- Common/include/linear_algebra/CSysVector.hpp | 38 ++- Common/src/linear_algebra/CSysVectorGPU.cu | 246 +++++++++++-------- 2 files changed, 178 insertions(+), 106 deletions(-) diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 28d3121f2db..3306b952507 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -225,10 +225,15 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> } /*! - * \brief enum listing the possible scalar operators on GPU + * \brief enum listing the possible vector-scalar operators on GPU */ enum class GPUScalarOp { SET, ADD, SUB, MUL, DIV }; + /*! + * \brief enum listing the possible vector-vector operators on GPU + */ + enum class GPUVectorOp { SET, ADD, SUB, NEG }; + /*! * \brief method to launch the generic Unary Operation Kernel and apply given functors on GPU * \param[in] op - Unitary operation @@ -236,6 +241,13 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ void GPUUnaryOperation(GPUScalarOp op, ScalarType val); + /*! + * \brief method to launch a generic Binary Operation Kernel and apply given functors on GPU + * \param[in] op - Binary operation + * \param[in] other - other CsysVector + */ + void GPUBinaryOperation(GPUVectorOp op, const CSysVector& other); + /*! * \brief Performs the memory copy from host to device. * \param[in] trigger - boolean value that decides whether to conduct the transfer or not. True by default. @@ -320,12 +332,30 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] other - Another vector. */ CSysVector& operator=(const CSysVector& other) { +#ifdef HAVE_CUDA + GPUBinaryOperation(GPUVectorOp::SET, other); +#else CSYSVEC_PARFOR for (auto i = 0ul; i < nElm; ++i) vec_val[i] = other.vec_val[i]; END_CSYSVEC_PARFOR +#endif return *this; } +#ifdef HAVE_CUDA + /*! + * \brief CSysVector-CSysVector overloaded operators on GPUs + */ + CSysVector& operator+=(const CSysVector& other) { + GPUBinaryOperation(GPUVectorOp::ADD, other); + return *this; + } + CSysVector& operator-=(const CSysVector& other) { + GPUBinaryOperation(GPUVectorOp::SUB, other); + return *this; + } +#endif + /*! * \brief Compound assignement operations with scalars and expressions. * \param[in] val/expr - Scalar value or expression. @@ -344,7 +374,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> return *this; \ } #else -#define MAKE_COMPOUND(OP, TAG) \ +#define MAKE_COMPOUND(OP, TAG) \ CSysVector& operator OP(ScalarType val) { \ CSYSVEC_PARFOR \ for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val; \ @@ -360,7 +390,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> } #endif - MAKE_COMPOUND(=, SET) + MAKE_COMPOUND(=, SET) MAKE_COMPOUND(+=, ADD) MAKE_COMPOUND(-=, SUB) MAKE_COMPOUND(*=, MUL) @@ -384,7 +414,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> ScalarType sum = 0.0; #ifdef HAVE_CUDA - sum = GPUDot(static_cast(expr.derived())); //assuming expr is a vector + sum = GPUDot(static_cast(expr.derived())); // assuming expr is a vector #else /*--- Local dot product for each thread. ---*/ diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 103b41a8851..ad1c7258c31 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -29,60 +29,102 @@ #include "../../include/linear_algebra/GPUComms.cuh" /*! + * \brief namespace defining the vector-vector operators for GPU + */ +namespace { +template struct OpSetVec { __device__ __forceinline__ T operator()(T, T b) const { return b; } }; +template struct OpAddVec { __device__ __forceinline__ T operator()(T a, T b) const { return a + b; } }; +template struct OpSubVec { __device__ __forceinline__ T operator()(T a, T b) const { return a - b; } }; +} + +template +__global__ void GPUBinaryOperationKernel(ScalarType* __restrict__ out, const ScalarType* __restrict__ other, + unsigned long n, Operator op) { + const unsigned long idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < n) out[idx] = op(out[idx], other[idx]); +} + +/*! + * \brief generic Binary Operation Kernel to apply given functors on GPU + */ +template +void CSysVector::GPUBinaryOperation(GPUVectorOp op, const CSysVector& other) { + + dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); + int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain); + dim3 gridDim(numBlocks, 1, 1); + + switch (op) { + case GPUVectorOp::SET: + GPUBinaryOperationKernel<<>>(this->d_vec_val, other.d_vec_val, this->nElmDomain, OpSetVec{}); + break; + case GPUVectorOp::ADD: + GPUBinaryOperationKernel<<>>(this->d_vec_val, other.d_vec_val, this->nElmDomain, OpAddVec{}); + break; + case GPUVectorOp::SUB: + GPUBinaryOperationKernel<<>>(this->d_vec_val, other.d_vec_val, this->nElmDomain, OpSubVec{}); + break; + } + gpuErrChk(cudaPeekAtLastError()); +} + +/*! + * \brief GPU dot prodcut kernel * \brief block-level reduction of elementwise products, accumulated into a single device scalar. */ template __global__ void GPUDotKernel(const ScalarType* __restrict__ a, const ScalarType* __restrict__ b, unsigned long n, ScalarType* __restrict__ result) { - //shared memory - extern __shared__ unsigned char smem_raw[]; - ScalarType* sdata = reinterpret_cast(smem_raw); - - //local and global thread indexes for access and block reduction - const unsigned long tid = threadIdx.x; - unsigned long idx = blockIdx.x * blockDim.x + tid; - const unsigned long stride = blockDim.x * gridDim.x; - - //thread reduction - ScalarType local = ScalarType(0); - for (; idx < n; idx += stride) local += a[idx] * b[idx]; - - sdata[tid] = local; - __syncthreads(); - - //block reduction - for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) { - if (tid < s) sdata[tid] += sdata[tid + s]; - __syncthreads(); - } - - //final atomic add per block - if (tid == 0) atomicAdd(result, sdata[0]); + //shared memory + extern __shared__ unsigned char smem_raw[]; + ScalarType* sdata = reinterpret_cast(smem_raw); + + //local and global thread indexes for access and block reduction + const unsigned long tid = threadIdx.x; + unsigned long idx = blockIdx.x * blockDim.x + tid; + const unsigned long stride = blockDim.x * gridDim.x; + + //thread reduction + ScalarType local = ScalarType(0); + for (; idx < n; idx += stride) local += a[idx] * b[idx]; + + sdata[tid] = local; + __syncthreads(); + + //block reduction + for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) { + if (tid < s) sdata[tid] += sdata[tid + s]; + __syncthreads(); + } + + //final atomic add per block + if (tid == 0) atomicAdd(result, sdata[0]); } +/*! + * \brief GPU dot product method between this and other CSysVector + */ template ScalarType CSysVector::GPUDot(const CSysVector& other) const { - dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); - int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain); - dim3 gridDim(numBlocks, 1, 1); - - // allocate and zero the result scalar - ScalarType* d_dot_result; - gpuErrChk(cudaMalloc(&d_dot_result, sizeof(ScalarType))); - gpuErrChk(cudaMemset(d_dot_result, 0, sizeof(ScalarType))); - - const size_t sharedBytes = KernelParameters::MVP_BLOCK_SIZE * sizeof(ScalarType); - GPUDotKernel<<>>(this->d_vec_val, other.d_vec_val, this->nElmDomain, d_dot_result); - gpuErrChk(cudaPeekAtLastError()); + dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); + int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain); + dim3 gridDim(numBlocks, 1, 1); - ScalarType result; - gpuErrChk(cudaMemcpy(&result, d_dot_result, sizeof(ScalarType), cudaMemcpyDeviceToHost)); - gpuErrChk(cudaFree(d_dot_result)); + // allocate and zero the result scalar + ScalarType* d_dot_result; + gpuErrChk(cudaMalloc(&d_dot_result, sizeof(ScalarType))); + gpuErrChk(cudaMemset(d_dot_result, 0, sizeof(ScalarType))); + const size_t sharedBytes = KernelParameters::MVP_BLOCK_SIZE * sizeof(ScalarType); + GPUDotKernel<<>>(this->d_vec_val, other.d_vec_val, this->nElmDomain, d_dot_result); + gpuErrChk(cudaPeekAtLastError()); + ScalarType result; + gpuErrChk(cudaMemcpy(&result, d_dot_result, sizeof(ScalarType), cudaMemcpyDeviceToHost)); + gpuErrChk(cudaFree(d_dot_result)); - return result; + return result; } /*! @@ -101,36 +143,36 @@ template struct OpDivScalar { T val; __device__ __forceinline__ T oper */ template __global__ void GPUUnaryOperationKernel(ScalarType* __restrict__ vec, unsigned long n, Operator op) { - const unsigned long idx = blockIdx.x * blockDim.x + threadIdx.x; - if (idx < n) vec[idx] = op(vec[idx]); + const unsigned long idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < n) vec[idx] = op(vec[idx]); } template void CSysVector::GPUUnaryOperation(GPUScalarOp op, ScalarType val) { - dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); - int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain); - dim3 gridDim(numBlocks, 1, 1); - - switch (op) { - case GPUScalarOp::SET: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSetScalar{val}); - break; - case GPUScalarOp::ADD: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpAddScalar{val}); - break; - case GPUScalarOp::SUB: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSubScalar{val}); - break; - case GPUScalarOp::MUL: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpMulScalar{val}); - break; - case GPUScalarOp::DIV: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpDivScalar{val}); - break; - } - gpuErrChk(cudaPeekAtLastError()); - gpuErrChk(cudaDeviceSynchronize()); + dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); + int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain); + dim3 gridDim(numBlocks, 1, 1); + + switch (op) { + case GPUScalarOp::SET: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSetScalar{val}); + break; + case GPUScalarOp::ADD: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpAddScalar{val}); + break; + case GPUScalarOp::SUB: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSubScalar{val}); + break; + case GPUScalarOp::MUL: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpMulScalar{val}); + break; + case GPUScalarOp::DIV: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpDivScalar{val}); + break; + } + gpuErrChk(cudaPeekAtLastError()); + gpuErrChk(cudaDeviceSynchronize()); } @@ -162,49 +204,49 @@ __global__ void GPUmultiDot(const ScalarType* const* __restrict__ d_V, const siz const ScalarType* const* __restrict__ d_W, const size_t m, const size_t size, ScalarType* __restrict__ d_local) { - // Map each x,y block to the specific (i,j) dot product - const size_t pair_idx = blockIdx.y; - if (pair_idx >= n * m) return; + // Map each x,y block to the specific (i,j) dot product + const size_t pair_idx = blockIdx.y; + if (pair_idx >= n * m) return; - const size_t i = pair_idx / m; - const size_t j = pair_idx % m; + const size_t i = pair_idx / m; + const size_t j = pair_idx % m; - //get the corresponding vectors - const ScalarType* __restrict__ vi = d_V[i]; - const ScalarType* __restrict__ wj = d_W[j]; + //get the corresponding vectors + const ScalarType* __restrict__ vi = d_V[i]; + const ScalarType* __restrict__ wj = d_W[j]; - // grid strided loop over the vector elements - ScalarType local_sum = 0.0; - const size_t tid = blockIdx.x * blockDim.x + threadIdx.x; - const size_t stride = gridDim.x * blockDim.x; + // grid strided loop over the vector elements + ScalarType local_sum = 0.0; + const size_t tid = blockIdx.x * blockDim.x + threadIdx.x; + const size_t stride = gridDim.x * blockDim.x; - for (size_t k = tid; k < size; k += stride) - { - local_sum += vi[k] * wj[k]; - } + for (size_t k = tid; k < size; k += stride) + { + local_sum += vi[k] * wj[k]; + } - // shared memory reduction within the block - extern __shared__ char shared_mem[]; - ScalarType* sdata = reinterpret_cast(shared_mem); - - sdata[threadIdx.x] = local_sum; - __syncthreads(); - - // parallel reduction on the block - for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) - { - if (threadIdx.x < s) - { - sdata[threadIdx.x] += sdata[threadIdx.x + s]; - } - __syncthreads(); - } + // shared memory reduction within the block + extern __shared__ char shared_mem[]; + ScalarType* sdata = reinterpret_cast(shared_mem); + + sdata[threadIdx.x] = local_sum; + __syncthreads(); + + // parallel reduction on the block + for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) + { + if (threadIdx.x < s) + { + sdata[threadIdx.x] += sdata[threadIdx.x + s]; + } + __syncthreads(); + } - // atomic add of each block partial sum to the output matrix, operated by thread 0 of each block - if (threadIdx.x == 0) - { - atomicAdd(&d_local[i * m + j], sdata[0]); - } + // atomic add of each block partial sum to the output matrix, operated by thread 0 of each block + if (threadIdx.x == 0) + { + atomicAdd(&d_local[i * m + j], sdata[0]); + } } template From 5fe25b814ff1c9b256a00a568f4c5e41be5c96ee Mon Sep 17 00:00:00 2001 From: ddg93 Date: Wed, 29 Jul 2026 11:20:36 +0200 Subject: [PATCH 8/9] introduced GPU compatibility in vector_expressions through abstract syntax tree. Binary vector expressions are coded into the tree and passed to the kernel for execution. Negative assignment is possible. validated on rae2822 --- Common/include/linear_algebra/CSysVector.hpp | 81 +++++---- Common/include/linear_algebra/gpu_ast.hpp | 71 ++++++++ .../linear_algebra/vector_expressions.hpp | 103 ++++++++--- Common/src/linear_algebra/CSysMatrixGPU.cu | 5 +- Common/src/linear_algebra/CSysVector.cpp | 1 - Common/src/linear_algebra/CSysVectorGPU.cu | 163 ++++++++++-------- 6 files changed, 295 insertions(+), 129 deletions(-) create mode 100644 Common/include/linear_algebra/gpu_ast.hpp diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 3306b952507..843ca29c1a6 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -34,6 +34,7 @@ #include "../parallelization/mpi_structure.hpp" #include "../parallelization/omp_structure.hpp" #include "../parallelization/vectorization.hpp" +#include "gpu_ast.hpp" #include "vector_expressions.hpp" #include "../../include/CConfig.hpp" @@ -225,15 +226,42 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> } /*! - * \brief enum listing the possible vector-scalar operators on GPU + * \brief enumeratethe possible vector-scalar operators on GPU */ enum class GPUScalarOp { SET, ADD, SUB, MUL, DIV }; /*! - * \brief enum listing the possible vector-vector operators on GPU + * \brief enumerate the possible vector-vector operators on GPU */ enum class GPUVectorOp { SET, ADD, SUB, NEG }; + /*! + * \brief map GPU scalar Ops to the appropriate GPU Operation code in the tree + * \note prevents SET=assignement operations raising an error + */ + inline GPUOpType ToASTCombineOp(GPUScalarOp op) { + switch (op) { + case GPUScalarOp::ADD: return GPUOpType::ADD; + case GPUScalarOp::SUB: return GPUOpType::SUB; + case GPUScalarOp::MUL: return GPUOpType::MUL; + case GPUScalarOp::DIV: return GPUOpType::DIV; + default: + SU2_MPI::Error("ToASTCombineOp called with GPUScalarOp::SET.", CURRENT_FUNCTION); + return GPUOpType::ADD; + } + } + + /*! + * \brief appends a new vector node to the Abstract Syntax Tree for the given payload + */ + int BuildAST(GPUASTPayload& payload) const { + //SyncToDevice(); + int idx = payload.NewNode(); + payload.nodes[idx].op = GPUOpType::VEC; + payload.nodes[idx].d_ptr = this->d_vec_val; + return idx; + } + /*! * \brief method to launch the generic Unary Operation Kernel and apply given functors on GPU * \param[in] op - Unitary operation @@ -241,13 +269,6 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ void GPUUnaryOperation(GPUScalarOp op, ScalarType val); - /*! - * \brief method to launch a generic Binary Operation Kernel and apply given functors on GPU - * \param[in] op - Binary operation - * \param[in] other - other CsysVector - */ - void GPUBinaryOperation(GPUVectorOp op, const CSysVector& other); - /*! * \brief Performs the memory copy from host to device. * \param[in] trigger - boolean value that decides whether to conduct the transfer or not. True by default. @@ -332,32 +353,14 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] other - Another vector. */ CSysVector& operator=(const CSysVector& other) { -#ifdef HAVE_CUDA - GPUBinaryOperation(GPUVectorOp::SET, other); -#else CSYSVEC_PARFOR for (auto i = 0ul; i < nElm; ++i) vec_val[i] = other.vec_val[i]; END_CSYSVEC_PARFOR -#endif return *this; } -#ifdef HAVE_CUDA /*! - * \brief CSysVector-CSysVector overloaded operators on GPUs - */ - CSysVector& operator+=(const CSysVector& other) { - GPUBinaryOperation(GPUVectorOp::ADD, other); - return *this; - } - CSysVector& operator-=(const CSysVector& other) { - GPUBinaryOperation(GPUVectorOp::SUB, other); - return *this; - } -#endif - - /*! - * \brief Compound assignement operations with scalars and expressions. + * \brief Compound assignement operations with scalars and expressions, GPU or CPU * \param[in] val/expr - Scalar value or expression. */ #ifdef HAVE_CUDA @@ -368,9 +371,18 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> } \ template \ CSysVector& operator OP(const VecExpr::CVecExpr& expr) { \ - CSYSVEC_PARFOR \ - for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP expr.derived()[i]; \ - END_CSYSVEC_PARFOR \ + GPUASTPayload payload; \ + int rhs_idx = expr.derived().BuildAST(payload); \ + int root_idx = rhs_idx; \ + if (GPUScalarOp::TAG != GPUScalarOp::SET) { \ + int self_idx = this->BuildAST(payload); \ + root_idx = payload.NewNode(); \ + payload.nodes[root_idx].op = ToASTCombineOp(GPUScalarOp::TAG); \ + payload.nodes[root_idx].left = self_idx; \ + payload.nodes[root_idx].right = rhs_idx; \ + } \ + payload.root_idx = root_idx; \ + LaunchGPUAST(payload, this->nElmDomain); \ return *this; \ } #else @@ -449,6 +461,13 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ ScalarType GPUDot(const CSysVector& other) const; + /*! + * \brief launch method for the generic Abstract Syntax Tree GPU kernel + * \param[in] payload - payload is the abstract syntax tree + * \param[in] nElm - number of elements in the vector + */ + void LaunchGPUAST(const GPUASTPayload& payload, unsigned long nElm); + /*! * \brief Computes the product of V^T W efficiencly, where V and W are tall matrices stored as vectors of CSysVector. * \param[in] V - Tall matrix. diff --git a/Common/include/linear_algebra/gpu_ast.hpp b/Common/include/linear_algebra/gpu_ast.hpp new file mode 100644 index 00000000000..b16f8970839 --- /dev/null +++ b/Common/include/linear_algebra/gpu_ast.hpp @@ -0,0 +1,71 @@ +/*! + * \file gpu_ast.hpp + * \brief definition of the Abstract Syntax Tree for evaluating CVecExpr expression trees on GPU with a single generic interpreter kernel + * \author D. Di giusto + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include +#include "../parallelization/mpi_structure.hpp" + +/*! + * \brief enumeration of the possible node-operations on GPU to build abstract syntax trees + */ +enum class GPUOpType { VEC, SCALAR, ADD, SUB, MUL, DIV, NEG }; + +/*! + * \brief GPU Abstract Syntax Tree node structure + */ +template +struct GPUASTNode { + GPUOpType op = GPUOpType::SCALAR; // operation for this node + int left = -1; // index of left child node-operation, -1 for leaf + int right = -1; // index of right child node-operation, -1 for leaf + const ScalarType* d_ptr = nullptr; // device pointer for vector + ScalarType val = ScalarType(0); // constant value for scalar +}; + +/*! + * \brief arbitrary max limit on the number of nodes in one Abstract Syntax tree + */ +constexpr int MAX_AST_NODES = 16; + +/*! + * \brief struct for the abstract syntax tree payload sent from CPU to GPU during one kernel launch + */ +template +struct GPUASTPayload { + GPUASTNode nodes[MAX_AST_NODES]; // array hosting all nodes for a single expression-operator + int root_idx = 0; + int node_count = 0; + + // allocation method + int NewNode() { + if (node_count >= MAX_AST_NODES) { + SU2_MPI::Error("GPU AST node budget exceeded, you should increase MAX_AST_NODES!", __FUNCTION__); + } + return node_count++; + } +}; \ No newline at end of file diff --git a/Common/include/linear_algebra/vector_expressions.hpp b/Common/include/linear_algebra/vector_expressions.hpp index a0d0ce28901..04662162bfe 100644 --- a/Common/include/linear_algebra/vector_expressions.hpp +++ b/Common/include/linear_algebra/vector_expressions.hpp @@ -35,6 +35,10 @@ #include #include +#ifdef HAVE_CUDA +#include "gpu_ast.hpp" +#endif + namespace VecExpr { /// \addtogroup VecExpr /// @{ @@ -78,6 +82,19 @@ class Bcast : public CVecExpr, Scalar> { static constexpr bool StoreAsRef = false; FORCEINLINE Bcast(const Scalar& x_) : x(x_) {} FORCEINLINE const Scalar& operator[](size_t) const { return x; } + +/*! + * \brief appends a new scalar node to the abstract expression tree of a given payload + */ +#ifdef HAVE_CUDA + int BuildAST(GPUASTPayload& payload) const { + int idx = payload.NewNode(); + payload.nodes[idx].op = GPUOpType::SCALAR; + payload.nodes[idx].val = x; + return idx; + } +#endif + }; /*! @@ -117,10 +134,29 @@ namespace math = ::std; #define RETURNS(...) \ ->decltype(__VA_ARGS__) { return __VA_ARGS__; } +/*--- Macro to inject GPU abstract syntax trees into unary functions + * The corresponding NONE macro is for operators not yet existing ---*/ + + +#ifdef HAVE_CUDA +#define GPU_UNARY_AST(OPTYPE) \ + int BuildAST(GPUASTPayload& payload) const { \ + int child = u.BuildAST(payload); \ + int idx = payload.NewNode(); \ + payload.nodes[idx].op = GPUOpType::OPTYPE; \ + payload.nodes[idx].left = child; \ + return idx; \ + } +#define GPU_UNARY_AST_NONE() +#else +#define GPU_UNARY_AST(OPTYPE) +#define GPU_UNARY_AST_NONE() +#endif + /*--- Macro to create expression classes (EXPR) and overloads (FUN) for unary * functions, based on their coefficient-wise implementation (IMPL). ---*/ -#define MAKE_UNARY_FUN(FUN, EXPR, IMPL) \ +#define MAKE_UNARY_FUN(FUN, EXPR, IMPL, GPU_AST) \ /*!--- Expression class. ---*/ \ template \ class EXPR : public CVecExpr, Scalar> { \ @@ -130,25 +166,49 @@ namespace math = ::std; static constexpr bool StoreAsRef = false; \ FORCEINLINE EXPR(const U& u_) : u(u_) {} \ FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i])) \ + GPU_AST \ }; \ /*!--- Function overload, returns an expression object. ---*/ \ template \ FORCEINLINE auto FUN(const CVecExpr& u) RETURNS(EXPR(u.derived())) #define sign_impl(x) Scalar(1 - 2 * (x < 0)) -MAKE_UNARY_FUN(operator-, minus_, -) -MAKE_UNARY_FUN(abs, abs_, math::abs) -MAKE_UNARY_FUN(exp, exp_, math::exp) -MAKE_UNARY_FUN(sqrt, sqrt_, math::sqrt) -MAKE_UNARY_FUN(sign, sign_, sign_impl) +MAKE_UNARY_FUN(operator-, minus_, -, GPU_UNARY_AST(NEG)) +MAKE_UNARY_FUN(abs, abs_, math::abs, GPU_UNARY_AST_NONE()) +MAKE_UNARY_FUN(exp, exp_, math::exp, GPU_UNARY_AST_NONE()) +MAKE_UNARY_FUN(sqrt, sqrt_, math::sqrt, GPU_UNARY_AST_NONE()) +MAKE_UNARY_FUN(sign, sign_, sign_impl, GPU_UNARY_AST_NONE()) #undef sign_impl #undef MAKE_UNARY_FUN +#undef GPU_UNARY_AST +#undef GPU_UNARY_AST_NONE + + +/*--- Macro to inject GPU abstract syntax trees into binary functions + * The corresponding NONE macro is for operators not yet existing ---*/ + +#ifdef HAVE_CUDA +#define GPU_BINARY_AST(OPTYPE) \ + int BuildAST(GPUASTPayload& payload) const { \ + int l = u.BuildAST(payload); \ + int r = v.BuildAST(payload); \ + int idx = payload.NewNode(); \ + payload.nodes[idx].op = GPUOpType::OPTYPE; \ + payload.nodes[idx].left = l; \ + payload.nodes[idx].right = r; \ + return idx; \ + } +#define GPU_BINARY_AST_NONE() +#else +#define GPU_BINARY_AST(OPTYPE) +#define GPU_BINARY_AST_NONE() +#endif /*--- Macro to create expressions and overloads for binary functions. ---*/ // clang-format off -#define MAKE_BINARY_FUN(FUN, EXPR, IMPL) \ +#define MAKE_BINARY_FUN(FUN, EXPR, IMPL, GPU_AST) \ /*!--- Expression class. ---*/ \ template \ class EXPR : public CVecExpr, Scalar> { \ @@ -159,6 +219,7 @@ MAKE_UNARY_FUN(sign, sign_, sign_impl) static constexpr bool StoreAsRef = false; \ FORCEINLINE EXPR(const U& u_, const V& v_) : u(u_), v(v_) {} \ FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i], v[i])) \ + GPU_AST \ }; \ /*!--- Vector with vector function overload. ---*/ \ template \ @@ -191,9 +252,9 @@ using std::fmax; using std::fmin; #undef MAKE_FMINMAX_OVERLOADS -MAKE_BINARY_FUN(fmax, max_, fmax) -MAKE_BINARY_FUN(fmin, min_, fmin) -MAKE_BINARY_FUN(pow, pow_, math::pow) +MAKE_BINARY_FUN(fmax, max_, fmax, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(fmin, min_, fmin, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(pow, pow_, math::pow, GPU_BINARY_AST_NONE()) /*--- sts::plus and co. were tried, the code was horrendous (due to the forced * conversion between different types) and creating functions for these ops @@ -203,10 +264,10 @@ MAKE_BINARY_FUN(pow, pow_, math::pow) #define sub_impl(a, b) a - b #define mul_impl(a, b) a* b #define div_impl(a, b) a / b -MAKE_BINARY_FUN(operator+, add_, add_impl) -MAKE_BINARY_FUN(operator-, sub_, sub_impl) -MAKE_BINARY_FUN(operator*, mul_, mul_impl) -MAKE_BINARY_FUN(operator/, div_, div_impl) +MAKE_BINARY_FUN(operator+, add_, add_impl, GPU_BINARY_AST(ADD)) +MAKE_BINARY_FUN(operator-, sub_, sub_impl, GPU_BINARY_AST(SUB)) +MAKE_BINARY_FUN(operator*, mul_, mul_impl, GPU_BINARY_AST(MUL)) +MAKE_BINARY_FUN(operator/, div_, div_impl, GPU_BINARY_AST(DIV)) #undef add_impl #undef sub_impl #undef mul_impl @@ -224,12 +285,12 @@ MAKE_BINARY_FUN(operator/, div_, div_impl) #define ne_impl(a, b) TO_PASSIVE(a != b) #define lt_impl(a, b) TO_PASSIVE(a < b) #define gt_impl(a, b) TO_PASSIVE(a > b) -MAKE_BINARY_FUN(operator<=, le_, le_impl) -MAKE_BINARY_FUN(operator>=, ge_, ge_impl) -MAKE_BINARY_FUN(operator==, eq_, eq_impl) -MAKE_BINARY_FUN(operator!=, ne_, ne_impl) -MAKE_BINARY_FUN(operator<, lt_, lt_impl) -MAKE_BINARY_FUN(operator>, gt_, gt_impl) +MAKE_BINARY_FUN(operator<=, le_, le_impl, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(operator>=, ge_, ge_impl, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(operator==, eq_, eq_impl, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(operator!=, ne_, ne_impl, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(operator<, lt_, lt_impl, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(operator>, gt_, gt_impl, GPU_BINARY_AST_NONE()) #undef TO_PASSIVE #undef le_impl #undef ge_impl @@ -239,6 +300,8 @@ MAKE_BINARY_FUN(operator>, gt_, gt_impl) #undef gt_impl #undef MAKE_BINARY_FUN +#undef GPU_BINARY_AST +#undef GPU_BINARY_AST_NONE /// @} } // namespace VecExpr diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 29c97002357..712ae29357f 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -109,16 +109,15 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector ScalarType* d_vec = vec.GetDevicePointer(); ScalarType* d_prod = prod.GetDevicePointer(); - //vec.HtDTransfer(); - + dim3 blockDim(static_cast(nVar), 1, 1); dim3 gridDim(static_cast(nPointDomain), 1, 1); BlockLDU_SpMV_kernel<<>>( nPointDomain, nVar, gpu.row_ptr_l, gpu.col_ind_l, gpu.l, gpu.d, gpu.row_ptr_u, gpu.col_ind_u, gpu.u, d_vec, d_prod); gpuErrChk(cudaGetLastError()); + gpuErrChk(cudaDeviceSynchronize()); - //prod.DtHTransfer(); } template diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp index 09d3740fb0c..8dcb5dad909 100644 --- a/Common/src/linear_algebra/CSysVector.cpp +++ b/Common/src/linear_algebra/CSysVector.cpp @@ -50,7 +50,6 @@ void CSysVector::Initialize(unsigned long numBlk, unsigned long numB omp_chunk_size = computeStaticChunkSize(nElm, omp_get_max_threads(), OMP_MAX_SIZE); - // if (vec_val == nullptr) vec_val = MemoryAllocation::aligned_alloc(64, nElm * sizeof(ScalarType)); #ifdef HAVE_CUDA useCuda = true; #endif diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index ad1c7258c31..40d74bce52f 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -27,47 +27,110 @@ #include "../../include/linear_algebra/CSysVector.hpp" #include "../../include/linear_algebra/GPUComms.cuh" +#include "../../include/linear_algebra/gpu_ast.hpp" + /*! - * \brief namespace defining the vector-vector operators for GPU + * \brief evaluates the abstract syntax tree iteratively. + * The tree is unrolled to prevent recursion and undefinite memory mapping. */ -namespace { -template struct OpSetVec { __device__ __forceinline__ T operator()(T, T b) const { return b; } }; -template struct OpAddVec { __device__ __forceinline__ T operator()(T a, T b) const { return a + b; } }; -template struct OpSubVec { __device__ __forceinline__ T operator()(T a, T b) const { return a - b; } }; +template +__device__ __forceinline__ ScalarType EvaluateASTIterative(const GPUASTPayload& payload, + unsigned long elem_idx) { + ScalarType values[MAX_AST_NODES]; + +#pragma unroll + for (int idx = 0; idx < payload.node_count; ++idx) { + const GPUASTNode& node = payload.nodes[idx]; + switch (node.op) { + case GPUOpType::VEC: values[idx] = node.d_ptr[elem_idx]; break; + case GPUOpType::SCALAR: values[idx] = node.val; break; + case GPUOpType::ADD: values[idx] = values[node.left] + values[node.right]; break; + case GPUOpType::SUB: values[idx] = values[node.left] - values[node.right]; break; + case GPUOpType::MUL: values[idx] = values[node.left] * values[node.right]; break; + case GPUOpType::DIV: values[idx] = values[node.left] / values[node.right]; break; + case GPUOpType::NEG: values[idx] = -values[node.left]; break; + } + } + return values[payload.root_idx]; } -template -__global__ void GPUBinaryOperationKernel(ScalarType* __restrict__ out, const ScalarType* __restrict__ other, - unsigned long n, Operator op) { - const unsigned long idx = blockIdx.x * blockDim.x + threadIdx.x; - if (idx < n) out[idx] = op(out[idx], other[idx]); +/*! + * \brief kernel evaluating the abstract syntax tree for each thread + */ +template +__global__ void GPUASTInterpreterKernel(ScalarType* __restrict__ d_out, GPUASTPayload payload, unsigned long n) { + unsigned long idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < n) { + d_out[idx] = EvaluateASTIterative(payload, idx); + } +} + + +template +void CSysVector::LaunchGPUAST(const GPUASTPayload& payload, unsigned long nElm) { + + dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); + int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, nElm); + dim3 gridDim(numBlocks, 1, 1); + + GPUASTInterpreterKernel<<>>(this->d_vec_val, payload, nElm); + gpuErrChk(cudaPeekAtLastError()); + + gpuErrChk(cudaDeviceSynchronize()); + } /*! - * \brief generic Binary Operation Kernel to apply given functors on GPU + * \brief namespace defining the scalar operators for GPU */ +namespace { +template struct OpSetScalar { T val; __device__ __forceinline__ T operator()(T x) const { return val; } }; +template struct OpAddScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x + val; } }; +template struct OpSubScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x - val; } }; +template struct OpMulScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x * val; } }; +template struct OpDivScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x / val; } }; +} // namespace + +/*! + * \brief generic Unary Operation Kernel to apply given functors on GPU + */ +template +__global__ void GPUUnaryOperationKernel(ScalarType* __restrict__ vec, unsigned long n, Operator op) { + const unsigned long idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < n) vec[idx] = op(vec[idx]); +} + template -void CSysVector::GPUBinaryOperation(GPUVectorOp op, const CSysVector& other) { +void CSysVector::GPUUnaryOperation(GPUScalarOp op, ScalarType val) { dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain); dim3 gridDim(numBlocks, 1, 1); + switch (op) { - case GPUVectorOp::SET: - GPUBinaryOperationKernel<<>>(this->d_vec_val, other.d_vec_val, this->nElmDomain, OpSetVec{}); + case GPUScalarOp::SET: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSetScalar{val}); break; - case GPUVectorOp::ADD: - GPUBinaryOperationKernel<<>>(this->d_vec_val, other.d_vec_val, this->nElmDomain, OpAddVec{}); + case GPUScalarOp::ADD: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpAddScalar{val}); break; - case GPUVectorOp::SUB: - GPUBinaryOperationKernel<<>>(this->d_vec_val, other.d_vec_val, this->nElmDomain, OpSubVec{}); + case GPUScalarOp::SUB: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSubScalar{val}); + break; + case GPUScalarOp::MUL: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpMulScalar{val}); + break; + case GPUScalarOp::DIV: + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpDivScalar{val}); break; } gpuErrChk(cudaPeekAtLastError()); + gpuErrChk(cudaDeviceSynchronize()); } + /*! * \brief GPU dot prodcut kernel * \brief block-level reduction of elementwise products, accumulated into a single device scalar. @@ -103,6 +166,7 @@ __global__ void GPUDotKernel(const ScalarType* __restrict__ a, const ScalarType* /*! * \brief GPU dot product method between this and other CSysVector + * \brief this is a vectors read-only kernel returning a scalar on host */ template ScalarType CSysVector::GPUDot(const CSysVector& other) const { @@ -127,56 +191,6 @@ ScalarType CSysVector::GPUDot(const CSysVector& other) const { return result; } -/*! - * \brief namespace defining the scalar operators for GPU - */ -namespace { -template struct OpSetScalar { T val; __device__ __forceinline__ T operator()(T x) const { return val; } }; -template struct OpAddScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x + val; } }; -template struct OpSubScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x - val; } }; -template struct OpMulScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x * val; } }; -template struct OpDivScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x / val; } }; -} // namespace - -/*! - * \brief generic Unary Operation Kernel to apply given functors on GPU - */ -template -__global__ void GPUUnaryOperationKernel(ScalarType* __restrict__ vec, unsigned long n, Operator op) { - const unsigned long idx = blockIdx.x * blockDim.x + threadIdx.x; - if (idx < n) vec[idx] = op(vec[idx]); -} - -template -void CSysVector::GPUUnaryOperation(GPUScalarOp op, ScalarType val) { - - dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); - int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain); - dim3 gridDim(numBlocks, 1, 1); - - switch (op) { - case GPUScalarOp::SET: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSetScalar{val}); - break; - case GPUScalarOp::ADD: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpAddScalar{val}); - break; - case GPUScalarOp::SUB: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSubScalar{val}); - break; - case GPUScalarOp::MUL: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpMulScalar{val}); - break; - case GPUScalarOp::DIV: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpDivScalar{val}); - break; - } - gpuErrChk(cudaPeekAtLastError()); - gpuErrChk(cudaDeviceSynchronize()); - - -} - template void CSysVector::HtDTransfer(bool trigger) const { @@ -249,6 +263,9 @@ __global__ void GPUmultiDot(const ScalarType* const* __restrict__ d_V, const siz } } +/*! + * \brief multi vector dot produt method for GPU, this is a vectors-read only method that returns an array of scalars + */ template const su2matrix& CSysVector::multiDotGPU(const std::vector>& V, const size_t i0, const size_t n, @@ -267,7 +284,6 @@ const su2matrix& CSysVector::multiDotGPU(const std::vect h_V_ptrs[i] = V[i0 + i].data(); } for (size_t j = 0; j < m; ++j){ - //gpuErrChk(cudaMemAdvise(vec_val, nElm * sizeof(ScalarType), cudaMemAdviseSetReadMostly, device_id)); //if read only, could be good h_W_ptrs[j] = W[j].data(); } @@ -343,7 +359,6 @@ template void CSysVector::LinearCombinationGPU(const unsigned long n, const std::vector>& vs, const ScalarType* ws, CSysVector& v, bool inc) { - const unsigned long nElm = v.nElmDomain; dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE,1,1); int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, nElm); @@ -354,16 +369,16 @@ void CSysVector::LinearCombinationGPU(const unsigned long n, const s //prepare vectors pointers and corresponding weights, passing them by value WeightedVecs vs_ws = {}; for (int j = 0; j < rem; ++j) { - vs_ws.ptrs[j] = vs[i + j].data(); // already on device from multiDot - vs_ws.weights[j] = ws[i + j]; // plain array indexing, not ws(k) + //vs[i + j].HtDTransfer(); //ensure is on GPU, was copied surely in multiDot + vs_ws.ptrs[j] = vs[i + j].data(); //.GetDevicePointer(); // get the pointer + vs_ws.weights[j] = ws[i + j]; // plain array indexing, not ws(k) } //calculate the linear combination on GPU, handle more than 4 vectors through inc || i > 0 LinearCombinationKernel<<>>(v.data(), vs_ws, rem, nElm, inc || i > 0); gpuErrChk(cudaPeekAtLastError()); } - gpuErrChk(cudaDeviceSynchronize()); - + } template class CSysVector; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits. From 4b43fa0d365475d5e500ff50a3fa1c73c0d13976 Mon Sep 17 00:00:00 2001 From: ddg93 Date: Wed, 29 Jul 2026 20:54:05 +0200 Subject: [PATCH 9/9] fell back to CUDA MAnaged memory and introduced synchronization to Host in all the Host routines that access vectors to ensure correctness. Validated against rae2822 --- Common/include/linear_algebra/CSysVector.hpp | 161 ++++++++++++++++-- Common/include/linear_algebra/gpu_ast.hpp | 21 ++- .../linear_algebra/vector_expressions.hpp | 94 +++++----- Common/src/linear_algebra/CSysMatrixGPU.cu | 16 +- Common/src/linear_algebra/CSysVector.cpp | 27 +-- Common/src/linear_algebra/CSysVectorGPU.cu | 69 +++++--- 6 files changed, 264 insertions(+), 124 deletions(-) diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 843ca29c1a6..94aac89e113 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -80,6 +80,12 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> ScalarType* d_vec_val = nullptr; /*!< \brief Device Pointer to store the vector values on the GPU. */ bool vec_is_managed = false; /*!< \brief Boolean that indicates whether GPU supports Unified Memory or not */ bool useCuda = false; /*!< \brief Whether CUDA is enabled. */ + enum class GPUMemState { + HOST, + DEVICE, + SYNCED + }; /*!< \brief enumerate the possible memory synchronization states for this */ + mutable GPUMemState memState = GPUMemState::HOST; /*!< \brief mutable GPUMemState initialized to HOST */ #ifdef HAVE_OMP mutable std::unique_ptr @@ -167,7 +173,12 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \note Not defined for expressions because we do not know their sizes. * \param[in] u - Vector being copied. */ - CSysVector(const CSysVector& u) { Initialize(u.GetNBlk(), u.GetNBlkDomain(), u.nVar, u.vec_val, true); } + CSysVector(const CSysVector& u) { +#ifdef HAVE_CUDA + u.SyncToHost(); +#endif + Initialize(u.GetNBlk(), u.GetNBlkDomain(), u.nVar, u.vec_val, true); + } /*! * \brief Swap contents with another vector. @@ -180,6 +191,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> std::swap(nElmDomain, other.nElmDomain); std::swap(nVar, other.nVar); std::swap(dot_scratch, other.dot_scratch); + std::swap(memState, other.memState); } /*! @@ -213,7 +225,9 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> void PassiveCopy(const CSysVector& other) { /*--- This is a method and not the overload of an operator to make sure who * calls it knows the consequence to the derivative information (lost) ---*/ - +#ifdef HAVE_CUDA + other.SyncToHost(); +#endif /*--- check if self-assignment, otherwise perform deep copy ---*/ if ((const void*)this == (const void*)&other) return; @@ -223,6 +237,9 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> CSYSVEC_PARFOR for (auto i = 0ul; i < nElm; i++) vec_val[i] = SU2_TYPE::GetValue(other[i]); END_CSYSVEC_PARFOR +#ifdef HAVE_CUDA + MarkHostDirty(); +#endif } /*! @@ -241,10 +258,14 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ inline GPUOpType ToASTCombineOp(GPUScalarOp op) { switch (op) { - case GPUScalarOp::ADD: return GPUOpType::ADD; - case GPUScalarOp::SUB: return GPUOpType::SUB; - case GPUScalarOp::MUL: return GPUOpType::MUL; - case GPUScalarOp::DIV: return GPUOpType::DIV; + case GPUScalarOp::ADD: + return GPUOpType::ADD; + case GPUScalarOp::SUB: + return GPUOpType::SUB; + case GPUScalarOp::MUL: + return GPUOpType::MUL; + case GPUScalarOp::DIV: + return GPUOpType::DIV; default: SU2_MPI::Error("ToASTCombineOp called with GPUScalarOp::SET.", CURRENT_FUNCTION); return GPUOpType::ADD; @@ -255,7 +276,6 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \brief appends a new vector node to the Abstract Syntax Tree for the given payload */ int BuildAST(GPUASTPayload& payload) const { - //SyncToDevice(); int idx = payload.NewNode(); payload.nodes[idx].op = GPUOpType::VEC; payload.nodes[idx].d_ptr = this->d_vec_val; @@ -295,7 +315,12 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> /*! * \brief return pointer that points to the CSysVector values in CPU memory */ - ScalarType* data() const { return vec_val; } + ScalarType* data() const { +#ifdef HAVE_CUDA + SyncToHost(); +#endif + return vec_val; + } /*! * \brief return the number of local elements in the CSysVector @@ -327,14 +352,35 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] i - Local index to access. * \return Value at position i. */ - inline ScalarType& operator[](unsigned long i) { return vec_val[i]; } - inline const ScalarType& operator[](unsigned long i) const { return vec_val[i]; } + inline ScalarType& operator[](unsigned long i) { +#ifdef HAVE_CUDA + SyncToHost(); + MarkHostDirty(); +#endif + return vec_val[i]; + } + inline const ScalarType& operator[](unsigned long i) const { +#ifdef HAVE_CUDA + SyncToHost(); +#endif + return vec_val[i]; + } /*! * \brief Iterators for range for loops. */ - inline const ScalarType* begin() const { return vec_val; } - inline const ScalarType* end() const { return vec_val + nElm; } + inline const ScalarType* begin() const { +#ifdef HAVE_CUDA + SyncToHost(); +#endif + return vec_val; + } + inline const ScalarType* end() const { +#ifdef HAVE_CUDA + SyncToHost(); +#endif + return vec_val + nElm; + } /*! * \brief Access operator with assignment permitted block version. @@ -342,8 +388,17 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] iVar - Index of variable. * \return Value at position (i,j). */ - inline ScalarType& operator()(unsigned long iPoint, unsigned long iVar) { return vec_val[iPoint * nVar + iVar]; } + inline ScalarType& operator()(unsigned long iPoint, unsigned long iVar) { +#ifdef HAVE_CUDA + SyncToHost(); + MarkHostDirty(); +#endif + return vec_val[iPoint * nVar + iVar]; + } inline const ScalarType& operator()(unsigned long iPoint, unsigned long iVar) const { +#ifdef HAVE_CUDA + SyncToHost(); +#endif return vec_val[iPoint * nVar + iVar]; } @@ -353,9 +408,15 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] other - Another vector. */ CSysVector& operator=(const CSysVector& other) { +#ifdef HAVE_CUDA + other.SyncToHost(); +#endif CSYSVEC_PARFOR for (auto i = 0ul; i < nElm; ++i) vec_val[i] = other.vec_val[i]; END_CSYSVEC_PARFOR +#ifdef HAVE_CUDA + this->MarkHostDirty(); +#endif return *this; } @@ -382,7 +443,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> payload.nodes[root_idx].right = rhs_idx; \ } \ payload.root_idx = root_idx; \ - LaunchGPUAST(payload, this->nElmDomain); \ + LaunchGPUAST(payload, this->nElm); \ return *this; \ } #else @@ -504,6 +565,37 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> static void LinearCombinationGPU(const unsigned long n, const std::vector>& vs, const ScalarType* ws, CSysVector& v, bool inc = false); + /*! + * \brief method to mark the vector as lastly modified by device + */ + void MarkDeviceDirty() const { memState = GPUMemState::DEVICE; } + + /*! + * \brief method to mark the vector as lastly modified by host + */ + void MarkHostDirty() const { memState = GPUMemState::HOST; } + + /*! + * \brief wrapper around HtDTransfer to manage memory state + */ + void SyncToDevice() const { + if (memState == GPUMemState::HOST) { + HtDTransfer(); + memState = GPUMemState::SYNCED; + } + } + /*! + * \brief wrapper around DtHTransfer to manage memory state + */ + void SyncToHost() const { + if (memState == GPUMemState::DEVICE) { + DtHTransfer(); + memState = GPUMemState::SYNCED; + } + } + + void print_memory_state() const { cout << "memory state: " << static_cast(memState) << endl; } + /*! * \brief Squared L2 norm of the vector (via dot with self). * \return Squared L2 norm. @@ -521,15 +613,32 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] iPoint - Index of block. * \return Pointer to start of block. */ - inline ScalarType* GetBlock(unsigned long iPoint) { return &vec_val[iPoint * nVar]; } - inline const ScalarType* GetBlock(unsigned long iPoint) const { return &vec_val[iPoint * nVar]; } + inline ScalarType* GetBlock(unsigned long iPoint) { +#ifdef HAVE_CUDA + SyncToHost(); + MarkHostDirty(); +#endif + return &vec_val[iPoint * nVar]; + } + inline const ScalarType* GetBlock(unsigned long iPoint) const { +#ifdef HAVE_CUDA + SyncToHost(); +#endif + return &vec_val[iPoint * nVar]; + } /*! * \brief Set the values to zero for one block. * \param[in] iPoint - Index of the block being set to zero. */ inline void SetBlock_Zero(unsigned long iPoint) { +#ifdef HAVE_CUDA + SyncToHost(); +#endif for (auto iVar = 0ul; iVar < nVar; iVar++) vec_val[iPoint * nVar + iVar] = 0.0; +#ifdef HAVE_CUDA + MarkHostDirty(); +#endif } /*! @@ -541,11 +650,17 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ template FORCEINLINE void SetBlock(unsigned long iPoint, const VectorType& block, ScalarType alpha = 1) { +#ifdef HAVE_CUDA + SyncToHost(); +#endif if (Overwrite) { for (auto i = 0ul; i < nVar; ++i) vec_val[iPoint * nVar + i] = alpha * block[i]; } else { for (auto i = 0ul; i < nVar; ++i) vec_val[iPoint * nVar + i] += alpha * block[i]; } +#ifdef HAVE_CUDA + MarkHostDirty(); +#endif } /*! @@ -588,13 +703,18 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> assert(nVar == this->nVar); ScalarType vec[N][nVar]; UnpackBlock(vector, mask, vec); - +#ifdef HAVE_CUDA + SyncToHost(); +#endif /*--- Update one by one skipping if mask is 0. ---*/ for (size_t k = 0; k < N; ++k) { if (mask[k] == 0) continue; SU2_OMP_SIMD for (size_t i = 0; i < nVar; ++i) vec_val[iPoint[k] * nVar + i] = vec[k][i]; } +#ifdef HAVE_CUDA + MarkHostDirty(); +#endif } /*! @@ -609,7 +729,9 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> assert(nVar == this->nVar); ScalarType vec[N][nVar]; UnpackBlock(vector, mask, vec); - +#ifdef HAVE_CUDA + SyncToHost(); +#endif /*--- Update one by one skipping if mask is 0. ---*/ for (size_t k = 0; k < N; ++k) { if (mask[k] == 0) continue; @@ -619,6 +741,9 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> vec_val[jPoint[k] * nVar + i] -= vec[k][i]; } } +#ifdef HAVE_CUDA + MarkHostDirty(); +#endif } }; diff --git a/Common/include/linear_algebra/gpu_ast.hpp b/Common/include/linear_algebra/gpu_ast.hpp index b16f8970839..9047d5dfae5 100644 --- a/Common/include/linear_algebra/gpu_ast.hpp +++ b/Common/include/linear_algebra/gpu_ast.hpp @@ -1,8 +1,7 @@ /*! * \file gpu_ast.hpp - * \brief definition of the Abstract Syntax Tree for evaluating CVecExpr expression trees on GPU with a single generic interpreter kernel - * \author D. Di giusto - * \version 8.5.0 "Harrier" + * \brief definition of the Abstract Syntax Tree for evaluating CVecExpr expression trees on GPU with a single generic + * interpreter kernel \author D. Di giusto \version 8.5.0 "Harrier" * * SU2 Project Website: https://su2code.github.io * @@ -40,11 +39,11 @@ enum class GPUOpType { VEC, SCALAR, ADD, SUB, MUL, DIV, NEG }; */ template struct GPUASTNode { - GPUOpType op = GPUOpType::SCALAR; // operation for this node - int left = -1; // index of left child node-operation, -1 for leaf - int right = -1; // index of right child node-operation, -1 for leaf - const ScalarType* d_ptr = nullptr; // device pointer for vector - ScalarType val = ScalarType(0); // constant value for scalar + GPUOpType op = GPUOpType::SCALAR; // operation for this node + int left = -1; // index of left child node-operation, -1 for leaf + int right = -1; // index of right child node-operation, -1 for leaf + const ScalarType* d_ptr = nullptr; // device pointer for vector + ScalarType val = ScalarType(0); // constant value for scalar }; /*! @@ -57,15 +56,15 @@ constexpr int MAX_AST_NODES = 16; */ template struct GPUASTPayload { - GPUASTNode nodes[MAX_AST_NODES]; // array hosting all nodes for a single expression-operator + GPUASTNode nodes[MAX_AST_NODES]; // array hosting all nodes for a single expression-operator int root_idx = 0; int node_count = 0; - // allocation method + // allocation method int NewNode() { if (node_count >= MAX_AST_NODES) { SU2_MPI::Error("GPU AST node budget exceeded, you should increase MAX_AST_NODES!", __FUNCTION__); } return node_count++; } -}; \ No newline at end of file +}; diff --git a/Common/include/linear_algebra/vector_expressions.hpp b/Common/include/linear_algebra/vector_expressions.hpp index 04662162bfe..e35379a3609 100644 --- a/Common/include/linear_algebra/vector_expressions.hpp +++ b/Common/include/linear_algebra/vector_expressions.hpp @@ -94,7 +94,6 @@ class Bcast : public CVecExpr, Scalar> { return idx; } #endif - }; /*! @@ -134,18 +133,17 @@ namespace math = ::std; #define RETURNS(...) \ ->decltype(__VA_ARGS__) { return __VA_ARGS__; } -/*--- Macro to inject GPU abstract syntax trees into unary functions +/*--- Macro to inject GPU abstract syntax trees into unary functions * The corresponding NONE macro is for operators not yet existing ---*/ - #ifdef HAVE_CUDA -#define GPU_UNARY_AST(OPTYPE) \ - int BuildAST(GPUASTPayload& payload) const { \ - int child = u.BuildAST(payload); \ - int idx = payload.NewNode(); \ - payload.nodes[idx].op = GPUOpType::OPTYPE; \ - payload.nodes[idx].left = child; \ - return idx; \ +#define GPU_UNARY_AST(OPTYPE) \ + int BuildAST(GPUASTPayload& payload) const { \ + int child = u.BuildAST(payload); \ + int idx = payload.NewNode(); \ + payload.nodes[idx].op = GPUOpType::OPTYPE; \ + payload.nodes[idx].left = child; \ + return idx; \ } #define GPU_UNARY_AST_NONE() #else @@ -156,48 +154,46 @@ namespace math = ::std; /*--- Macro to create expression classes (EXPR) and overloads (FUN) for unary * functions, based on their coefficient-wise implementation (IMPL). ---*/ -#define MAKE_UNARY_FUN(FUN, EXPR, IMPL, GPU_AST) \ - /*!--- Expression class. ---*/ \ - template \ - class EXPR : public CVecExpr, Scalar> { \ - store_t u; \ - \ - public: \ - static constexpr bool StoreAsRef = false; \ - FORCEINLINE EXPR(const U& u_) : u(u_) {} \ - FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i])) \ - GPU_AST \ - }; \ - /*!--- Function overload, returns an expression object. ---*/ \ - template \ +#define MAKE_UNARY_FUN(FUN, EXPR, IMPL, GPU_AST) \ + /*!--- Expression class. ---*/ \ + template \ + class EXPR : public CVecExpr, Scalar> { \ + store_t u; \ + \ + public: \ + static constexpr bool StoreAsRef = false; \ + FORCEINLINE EXPR(const U& u_) : u(u_) {} \ + FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i])) GPU_AST \ + }; \ + /*!--- Function overload, returns an expression object. ---*/ \ + template \ FORCEINLINE auto FUN(const CVecExpr& u) RETURNS(EXPR(u.derived())) #define sign_impl(x) Scalar(1 - 2 * (x < 0)) -MAKE_UNARY_FUN(operator-, minus_, -, GPU_UNARY_AST(NEG)) -MAKE_UNARY_FUN(abs, abs_, math::abs, GPU_UNARY_AST_NONE()) -MAKE_UNARY_FUN(exp, exp_, math::exp, GPU_UNARY_AST_NONE()) -MAKE_UNARY_FUN(sqrt, sqrt_, math::sqrt, GPU_UNARY_AST_NONE()) -MAKE_UNARY_FUN(sign, sign_, sign_impl, GPU_UNARY_AST_NONE()) +MAKE_UNARY_FUN(operator-, minus_, -, GPU_UNARY_AST(NEG)) +MAKE_UNARY_FUN(abs, abs_, math::abs, GPU_UNARY_AST_NONE()) +MAKE_UNARY_FUN(exp, exp_, math::exp, GPU_UNARY_AST_NONE()) +MAKE_UNARY_FUN(sqrt, sqrt_, math::sqrt, GPU_UNARY_AST_NONE()) +MAKE_UNARY_FUN(sign, sign_, sign_impl, GPU_UNARY_AST_NONE()) #undef sign_impl #undef MAKE_UNARY_FUN #undef GPU_UNARY_AST #undef GPU_UNARY_AST_NONE - -/*--- Macro to inject GPU abstract syntax trees into binary functions +/*--- Macro to inject GPU abstract syntax trees into binary functions * The corresponding NONE macro is for operators not yet existing ---*/ #ifdef HAVE_CUDA -#define GPU_BINARY_AST(OPTYPE) \ - int BuildAST(GPUASTPayload& payload) const { \ - int l = u.BuildAST(payload); \ - int r = v.BuildAST(payload); \ - int idx = payload.NewNode(); \ - payload.nodes[idx].op = GPUOpType::OPTYPE; \ - payload.nodes[idx].left = l; \ - payload.nodes[idx].right = r; \ - return idx; \ +#define GPU_BINARY_AST(OPTYPE) \ + int BuildAST(GPUASTPayload& payload) const { \ + int l = u.BuildAST(payload); \ + int r = v.BuildAST(payload); \ + int idx = payload.NewNode(); \ + payload.nodes[idx].op = GPUOpType::OPTYPE; \ + payload.nodes[idx].left = l; \ + payload.nodes[idx].right = r; \ + return idx; \ } #define GPU_BINARY_AST_NONE() #else @@ -252,9 +248,9 @@ using std::fmax; using std::fmin; #undef MAKE_FMINMAX_OVERLOADS -MAKE_BINARY_FUN(fmax, max_, fmax, GPU_BINARY_AST_NONE()) -MAKE_BINARY_FUN(fmin, min_, fmin, GPU_BINARY_AST_NONE()) -MAKE_BINARY_FUN(pow, pow_, math::pow, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(fmax, max_, fmax, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(fmin, min_, fmin, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(pow, pow_, math::pow, GPU_BINARY_AST_NONE()) /*--- sts::plus and co. were tried, the code was horrendous (due to the forced * conversion between different types) and creating functions for these ops @@ -285,12 +281,12 @@ MAKE_BINARY_FUN(operator/, div_, div_impl, GPU_BINARY_AST(DIV)) #define ne_impl(a, b) TO_PASSIVE(a != b) #define lt_impl(a, b) TO_PASSIVE(a < b) #define gt_impl(a, b) TO_PASSIVE(a > b) -MAKE_BINARY_FUN(operator<=, le_, le_impl, GPU_BINARY_AST_NONE()) -MAKE_BINARY_FUN(operator>=, ge_, ge_impl, GPU_BINARY_AST_NONE()) -MAKE_BINARY_FUN(operator==, eq_, eq_impl, GPU_BINARY_AST_NONE()) -MAKE_BINARY_FUN(operator!=, ne_, ne_impl, GPU_BINARY_AST_NONE()) -MAKE_BINARY_FUN(operator<, lt_, lt_impl, GPU_BINARY_AST_NONE()) -MAKE_BINARY_FUN(operator>, gt_, gt_impl, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(operator<=, le_, le_impl, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(operator>=, ge_, ge_impl, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(operator==, eq_, eq_impl, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(operator!=, ne_, ne_impl, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(operator<, lt_, lt_impl, GPU_BINARY_AST_NONE()) +MAKE_BINARY_FUN(operator>, gt_, gt_impl, GPU_BINARY_AST_NONE()) #undef TO_PASSIVE #undef le_impl #undef ge_impl diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 712ae29357f..0cd31d966c6 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -107,17 +107,18 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector SU2_MPI::Error("CUDA CSysMatrix block-LDU SpMV requires square blocks.", CURRENT_FUNCTION); } + vec.SyncToDevice(); ScalarType* d_vec = vec.GetDevicePointer(); ScalarType* d_prod = prod.GetDevicePointer(); - + dim3 blockDim(static_cast(nVar), 1, 1); dim3 gridDim(static_cast(nPointDomain), 1, 1); BlockLDU_SpMV_kernel<<>>( nPointDomain, nVar, gpu.row_ptr_l, gpu.col_ind_l, gpu.l, gpu.d, gpu.row_ptr_u, gpu.col_ind_u, gpu.u, d_vec, d_prod); gpuErrChk(cudaGetLastError()); - gpuErrChk(cudaDeviceSynchronize()); + prod.MarkDeviceDirty(); } template @@ -127,21 +128,26 @@ void CSysMatrix::GPUComputeJacobiPreconditioner(const CSysVector>>(d_invM, vec.data(), prod.data(), nPointDomain, nVar); + GPUMatrixVectorProductKernel<<>>(d_invM, d_vec, d_prod, nPointDomain, nVar); gpuErrChk( cudaPeekAtLastError() ); - gpuErrChk(cudaDeviceSynchronize()); + + prod.MarkDeviceDirty(); // this triggers () operators to sync data from device when preparing mpi buffers /*--- MPI Parallelization ---*/ CSysMatrixComms::Initiate(prod, geometry, config); CSysMatrixComms::Complete(prod, geometry, config); - gpuErrChk(cudaDeviceSynchronize()); + prod.SyncToDevice(); // this will trigger migration if mpi comms updated the host data } diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp index 8dcb5dad909..cbbdf1e3e4a 100644 --- a/Common/src/linear_algebra/CSysVector.cpp +++ b/Common/src/linear_algebra/CSysVector.cpp @@ -50,19 +50,16 @@ void CSysVector::Initialize(unsigned long numBlk, unsigned long numB omp_chunk_size = computeStaticChunkSize(nElm, omp_get_max_threads(), OMP_MAX_SIZE); + if (vec_val == nullptr) vec_val = MemoryAllocation::aligned_alloc(64, nElm * sizeof(ScalarType)); + #ifdef HAVE_CUDA useCuda = true; #endif - - if (useCuda && GPUMemoryAllocation::UMSupported()) { - vec_val = GPUMemoryAllocation::gpu_um_alloc(nElm * sizeof(ScalarType)); - d_vec_val = vec_val; // temporary alias for testing - vec_is_managed = true; - } else { - vec_val = MemoryAllocation::aligned_alloc(64, nElm * sizeof(ScalarType)); - if (useCuda) { + if (useCuda) { + if (d_vec_val == nullptr) { d_vec_val = GPUMemoryAllocation::gpu_alloc(nElm * sizeof(ScalarType)); } + MarkHostDirty(); } #ifdef HAVE_OMP @@ -153,18 +150,10 @@ CSysVector::~CSysVector() { if constexpr (!std::is_trivial_v) { for (auto i = 0ul; i < nElm; i++) vec_val[i].~ScalarType(); } - // MemoryAllocation::aligned_free(vec_val); - - // GPUMemoryAllocation::gpu_free(d_vec_val); + MemoryAllocation::aligned_free(vec_val); - if (useCuda && GPUMemoryAllocation::UMSupported()) { - GPUMemoryAllocation::gpu_free(vec_val); - d_vec_val = nullptr; - } else { - MemoryAllocation::aligned_free(vec_val); - if (useCuda) { - GPUMemoryAllocation::gpu_free(d_vec_val); - } + if (useCuda) { + GPUMemoryAllocation::gpu_free(d_vec_val); } } diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 40d74bce52f..be451800599 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -28,10 +28,11 @@ #include "../../include/linear_algebra/CSysVector.hpp" #include "../../include/linear_algebra/GPUComms.cuh" #include "../../include/linear_algebra/gpu_ast.hpp" - +#include // Needed for backtrace and backtrace_symbols +#include /*! - * \brief evaluates the abstract syntax tree iteratively. + * \brief evaluates the abstract syntax tree iteratively. * The tree is unrolled to prevent recursion and undefinite memory mapping. */ template @@ -66,19 +67,22 @@ __global__ void GPUASTInterpreterKernel(ScalarType* __restrict__ d_out, GPUASTPa } } - +/*! + * \brief method to configure and launch the kernel executing the abstract syntax tree for a given vector operators + */ template void CSysVector::LaunchGPUAST(const GPUASTPayload& payload, unsigned long nElm) { - + dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, nElm); dim3 gridDim(numBlocks, 1, 1); + this->SyncToDevice(); + GPUASTInterpreterKernel<<>>(this->d_vec_val, payload, nElm); gpuErrChk(cudaPeekAtLastError()); - gpuErrChk(cudaDeviceSynchronize()); - + this->MarkDeviceDirty(); } /*! @@ -101,33 +105,38 @@ __global__ void GPUUnaryOperationKernel(ScalarType* __restrict__ vec, unsigned l if (idx < n) vec[idx] = op(vec[idx]); } +/*! + * \brief dispatching method for the GPU unary operations + */ template void CSysVector::GPUUnaryOperation(GPUScalarOp op, ScalarType val) { dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1); - int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain); + int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElm); dim3 gridDim(numBlocks, 1, 1); + this->SyncToDevice(); switch (op) { case GPUScalarOp::SET: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSetScalar{val}); + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElm, OpSetScalar{val}); break; case GPUScalarOp::ADD: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpAddScalar{val}); + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElm, OpAddScalar{val}); break; case GPUScalarOp::SUB: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpSubScalar{val}); + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElm, OpSubScalar{val}); break; case GPUScalarOp::MUL: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpMulScalar{val}); + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElm, OpMulScalar{val}); break; case GPUScalarOp::DIV: - GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElmDomain, OpDivScalar{val}); + GPUUnaryOperationKernel<<>>(this->d_vec_val, this->nElm, OpDivScalar{val}); break; } gpuErrChk(cudaPeekAtLastError()); - gpuErrChk(cudaDeviceSynchronize()); + + this->MarkDeviceDirty(); } @@ -166,7 +175,7 @@ __global__ void GPUDotKernel(const ScalarType* __restrict__ a, const ScalarType* /*! * \brief GPU dot product method between this and other CSysVector - * \brief this is a vectors read-only kernel returning a scalar on host + * \note this is a vectors read-only kernel returning a scalar on host */ template ScalarType CSysVector::GPUDot(const CSysVector& other) const { @@ -175,6 +184,9 @@ ScalarType CSysVector::GPUDot(const CSysVector& other) const { int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain); dim3 gridDim(numBlocks, 1, 1); + SyncToDevice(); + other.SyncToDevice(); + // allocate and zero the result scalar ScalarType* d_dot_result; gpuErrChk(cudaMalloc(&d_dot_result, sizeof(ScalarType))); @@ -264,7 +276,8 @@ __global__ void GPUmultiDot(const ScalarType* const* __restrict__ d_V, const siz } /*! - * \brief multi vector dot produt method for GPU, this is a vectors-read only method that returns an array of scalars + * \brief multi vector dot produt method for GPU + * \note this is a vectors-read only method that returns an array of scalars */ template const su2matrix& CSysVector::multiDotGPU(const std::vector>& V, @@ -281,10 +294,12 @@ const su2matrix& CSysVector::multiDotGPU(const std::vect // get all the device pointers std::vector h_V_ptrs(n), h_W_ptrs(m); for (size_t i = 0; i < n; ++i){ - h_V_ptrs[i] = V[i0 + i].data(); + V[i0 + i].SyncToDevice(); + h_V_ptrs[i] = V[i0 + i].GetDevicePointer(); } for (size_t j = 0; j < m; ++j){ - h_W_ptrs[j] = W[j].data(); + W[j].SyncToDevice(); + h_W_ptrs[j] = W[j].GetDevicePointer(); } // copy the pointers to the device arrays of pointers @@ -339,6 +354,9 @@ struct WeightedVecs { ScalarType weights[N]; }; +/*! + * \brief linear combination kernel to calculate the next vector v from weights and vectors + */ template __global__ void LinearCombinationKernel(ScalarType* __restrict__ v, WeightedVecs wv, int n, unsigned long nElm, bool inc) @@ -355,6 +373,9 @@ __global__ void LinearCombinationKernel(ScalarType* __restrict__ v, WeightedVecs v[k] = result; } +/*! + * \brief dispatcher for the linear combination kernel on GPU + */ template void CSysVector::LinearCombinationGPU(const unsigned long n, const std::vector>& vs, const ScalarType* ws, CSysVector& v, bool inc) @@ -364,21 +385,25 @@ void CSysVector::LinearCombinationGPU(const unsigned long n, const s int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, nElm); dim3 gridDim(numBlocks, 1, 1); + v.SyncToDevice(); + ScalarType* d_v = v.GetDevicePointer(); + for (unsigned long i = 0; i < n; i += 4) { const int rem = static_cast(std::min(n - i, 4ul)); //prepare vectors pointers and corresponding weights, passing them by value WeightedVecs vs_ws = {}; for (int j = 0; j < rem; ++j) { - //vs[i + j].HtDTransfer(); //ensure is on GPU, was copied surely in multiDot - vs_ws.ptrs[j] = vs[i + j].data(); //.GetDevicePointer(); // get the pointer + vs[i + j].SyncToDevice(); //ensure is on GPU, was copied in multiDot + vs_ws.ptrs[j] = vs[i + j].GetDevicePointer(); // get the pointer vs_ws.weights[j] = ws[i + j]; // plain array indexing, not ws(k) } //calculate the linear combination on GPU, handle more than 4 vectors through inc || i > 0 - LinearCombinationKernel<<>>(v.data(), vs_ws, rem, nElm, inc || i > 0); + LinearCombinationKernel<<>>(d_v, vs_ws, rem, nElm, inc || i > 0); gpuErrChk(cudaPeekAtLastError()); } - gpuErrChk(cudaDeviceSynchronize()); - + + v.MarkDeviceDirty(); + } template class CSysVector; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits.