diff --git a/AUTHORS.md b/AUTHORS.md index fea157beba2..67ecbf541e5 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -69,6 +69,7 @@ Christian Bauer Clark Pederson Daumantas Kavolis Dave Taflin +Davide Di Giusto Eduardo Molina Edwin van der Weide Eitan Aberman diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index 2b6c67701ff..8bd7d515710 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -260,13 +260,13 @@ class CLU_SGSPreconditioner final : public CPreconditioner { * \param[out] v - CSysVector that is the result of the preconditioning. */ inline void operator()(const CSysVector& u, CSysVector& v) const override { - ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputeLU_SGSPreconditioner(u, v, geometry, config); }); + sparse_matrix.ComputeLU_SGSPreconditioner(u, v, geometry, config); } /*! * \note Also serves Q_LU_SGS: quantizes the diagonal blocks, no-op for plain LU_SGS. */ - inline void Build() override { sparse_matrix.QuantizeDiagonalBlocks(); } + inline void Build() override { sparse_matrix.BuildLU_SGSPreconditioner(); } }; /*! diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index abecc767e26..384b0bf6212 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -303,7 +303,7 @@ class CSysMatrix { LDU gpu; /*!< \brief Device matrix (all pointers to GPU memory). */ LDU ilu; /*!< \brief ILU factorization, host (values owned; pattern from geometry). */ LDU gpu_ilu; /*!< \brief ILU factorization, device (values and pattern in GPU memory). */ - ScalarType* d_invM = nullptr; /*!< \brief Device inverse diagonal blocks for the Jacobi preconditioner. */ + ScalarType* d_invM = nullptr; /*!< \brief Device inverse diagonal blocks for the Jacobi or LU-SGS preconditioner. */ /*--- Quantized off-diagonal storage (used when quantized_mode == true). ---*/ using QuantType = int8_t; @@ -356,7 +356,7 @@ class CSysMatrix { * rows in level k only depend on rows in levels < k. The same table drives the forward * (increasing level) and backward (decreasing level) substitution, because the U pattern is * the transpose of the L pattern. Used directly by the host/OMP substitution, and flattened - * into ilu_level_ptr / d_ilu_level_idx below for the GPU triangular solves. */ + * into ilu_level_ptr / d_precond_level_idx below for the GPU triangular solves. */ CCompressedSparsePatternUL levels_ilu; /*!< \brief Coloring of the (domain-only) ILU dependency graph, used only by the GPU iterative @@ -373,19 +373,22 @@ class CSysMatrix { vector ilu_color_ptr; /*!< \brief Start of each color in d_ilu_color_idx, size nColors+1. */ su2uint* d_ilu_color_idx = nullptr; /*!< \brief Row indices, grouped by color. */ - vector ilu_level_ptr; /*!< \brief Start of each level in d_ilu_level_idx, size nLevels+1. */ - su2uint* d_ilu_level_idx = nullptr; /*!< \brief Row indices, grouped by level. */ + vector precond_level_ptr; /*!< \brief Start of each level in d_precond_level_idx, size nLevels+1. */ + su2uint* d_precond_level_idx = nullptr; /*!< \brief Row indices, grouped by level. */ /*--- The per-color (factorization) and per-level (triangular solves) kernel launch sequences * are identical on every call: same grid/block sizes, same device pointers (all fixed members, * allocated once). Each is captured once into a CUDA graph and replayed to remove * host-side launch overhead without changing the parallelization. ---*/ mutable struct CUgraphExec_st* ilu_build_graph_exec = nullptr; - mutable struct CUgraphExec_st* ilu_apply_graph_exec = nullptr; - mutable const ScalarType* ilu_apply_graph_vec = nullptr; /*!< \brief Pointers the apply graph - * was captured with, to detect when - * it must be recaptured. */ - mutable ScalarType* ilu_apply_graph_prod = nullptr; + mutable struct CUgraphExec_st* precond_fwd_graph_exec = nullptr; // ILU or LU-SGS forward only + mutable struct CUgraphExec_st* precond_bwd_graph_exec = nullptr; // LU-SGS backward only + mutable const ScalarType* precond_fwd_graph_vec = nullptr; /*!< \brief Pointers the apply graph + * was captured with, to detect when + * it must be recaptured. */ + mutable ScalarType* precond_fwd_graph_prod = nullptr; + mutable ScalarType* precond_bwd_graph_prod = nullptr; + /*--- Non-default stream, needed for two mutually exclusive uses that never overlap on a given * matrix (quantized_mode and ILU are alternative preconditioner choices, decided once in * Initialize()): (1) the ILU build/apply CUDA graphs below, since the legacy default stream @@ -657,6 +660,31 @@ class CSysMatrix { */ void ComputeILUPreconditionerGPU(const CSysVector& vec, CSysVector& prod) const; + /*! + * \brief Build the LU-SGS preconditioner on the device + */ + void BuildLU_SGSPreconditionerGPU(); + + /*! + * \brief Apply the LU-SGS preconditioner forward pass on the device + */ + void ComputeLU_SGSForwardGPU(const CSysVector& vec, CSysVector& prod) const; + + /*! + * \brief Apply the LU-SGS preconditioner backward pass on the device + */ + void ComputeLU_SGSBackwardGPU(CSysVector& prod) const; + + /*! + * \brief Apply the forward pass of the LU-SGS preconditioner + */ + void ComputeLU_SGSPreconditionerForward(const CSysVector& vec, CSysVector& prod) const; + + /*! + * \brief Apply the backward pass of the LU-SGS preconditioner + */ + void ComputeLU_SGSPreconditionerBackward(CSysVector& prod) const; + public: /*! * \brief Constructor of the class. @@ -1230,6 +1258,11 @@ class CSysMatrix { void ComputeILUPreconditioner(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const; + /*! + * \brief Build the LU-SGS preconditioner. + */ + void BuildLU_SGSPreconditioner(); + /*! * \brief Multiply CSysVector by the preconditioner * \param[in] vec - CSysVector to be multiplied by the preconditioner. diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 84f30a2aa3d..e6292b4c10c 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -159,10 +159,11 @@ CSysMatrix::~CSysMatrix() { freeLDU(d_q_blocks); GPUMemoryAllocation::gpu_free(d_invM); GPUMemoryAllocation::gpu_free(d_ilu_color_idx); - GPUMemoryAllocation::gpu_free(d_ilu_level_idx); + GPUMemoryAllocation::gpu_free(d_precond_level_idx); #ifdef SU2_ENABLE_CUDA_KERNELS if (ilu_build_graph_exec != nullptr) cudaGraphExecDestroy(ilu_build_graph_exec); - if (ilu_apply_graph_exec != nullptr) cudaGraphExecDestroy(ilu_apply_graph_exec); + if (precond_fwd_graph_exec != nullptr) cudaGraphExecDestroy(precond_fwd_graph_exec); + if (precond_bwd_graph_exec != nullptr) cudaGraphExecDestroy(precond_bwd_graph_exec); if (aux_stream != nullptr) cudaStreamDestroy(aux_stream); if (htd_event != nullptr) cudaEventDestroy(htd_event); #endif @@ -217,6 +218,7 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi const bool ilu_needed = (prec == ILU); const bool diag_needed = (prec == JACOBI) || (prec == Q_JACOBI) || (prec == LINELET); + const bool lu_sgs_on_device = useCuda && (prec == LU_SGS || prec == Q_LU_SGS); /*--- Linelet also builds the Jacobi preconditioner but reads the inverse diagonal blocks on * the host, so only plain (or quantized) Jacobi can keep them exclusively on the device. ---*/ @@ -408,62 +410,72 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi if (diag_needed) allocAndInit(invM, nPointDomain * nVar * nEqn); - if (jacobi_on_device) { + const bool any_precond_on_device = useCuda && (jacobi_on_device || lu_sgs_on_device || ilu_needed); + + if (any_precond_on_device) { if (nVar != nEqn) { - SU2_MPI::Error("CUDA Jacobi preconditioner requires square blocks.", CURRENT_FUNCTION); + SU2_MPI::Error("CUDA preconditioners require square blocks.", CURRENT_FUNCTION); } if (nVar * nVar > 1024) { - SU2_MPI::Error("CUDA Jacobi preconditioner uses one thread per block entry, nVar is too large.", - CURRENT_FUNCTION); + SU2_MPI::Error("CUDA preconditioners use one thread per block entry, nVar is too large.", CURRENT_FUNCTION); } - d_invM = GPUMemoryAllocation::gpu_alloc(nPointDomain * nVar * nEqn * sizeof(ScalarType)); - } - if (useCuda && ilu_needed) { - if (nVar != nEqn) { - SU2_MPI::Error("CUDA ILU factorization requires square blocks.", CURRENT_FUNCTION); - } - if (nVar * nVar > 1024) { - SU2_MPI::Error("CUDA ILU factorization uses one thread per block entry, nVar is too large.", CURRENT_FUNCTION); + if (jacobi_on_device || lu_sgs_on_device) { + d_invM = GPUMemoryAllocation::gpu_alloc(nPointDomain * nVar * nEqn * sizeof(ScalarType)); } - /*--- The factors are built and used on the device, only the pattern and the level table - * are uploaded (once, here) because they do not change. ---*/ - gpu_ilu.nnz_l = ilu.nnz_l; - gpu_ilu.nnz_u = ilu.nnz_u; - GPUAllocAndInit(gpu_ilu.d, nPointDomain * nVar * nEqn); - GPUAllocAndInit(gpu_ilu.l, ilu.nnz_l * nVar * nEqn); - GPUAllocAndInit(gpu_ilu.u, ilu.nnz_u * nVar * nEqn); - GPUAllocAndCopy(gpu_ilu.row_ptr_l, ilu.row_ptr_l, nPointDomain + 1); - GPUAllocAndCopy(gpu_ilu.col_ind_l, ilu.col_ind_l, ilu.nnz_l); - GPUAllocAndCopy(gpu_ilu.row_ptr_u, ilu.row_ptr_u, nPointDomain + 1); - GPUAllocAndCopy(gpu_ilu.col_ind_u, ilu.col_ind_u, ilu.nnz_u); - - /*--- Flatten the coloring, the index type differs from the one of the pattern. It drives - * the factorization on the device. ---*/ - std::vector color_idx; - color_idx.reserve(nPointDomain); - ilu_color_ptr.clear(); - ilu_color_ptr.push_back(0); - for (auto color = 0ul; color < color_ilu.getOuterSize(); ++color) { - for (auto k = 0ul; k < color_ilu.getNumNonZeros(color); ++k) { - color_idx.push_back(static_cast(color_ilu.getInnerIdx(color, k))); + + /*--- Flattens a grouped sparse pattern (levels, colors) into a host ptr and device index arrays. + * Used in ILU levels and colors, and LU-SGS levels---*/ + auto FlattenGroupToDevice = [](const auto& grouped, std::vector& group_ptr, unsigned long reserveHint, + unsigned long bound = ~0ul) { + std::vector flat_idx; + flat_idx.reserve(reserveHint); + group_ptr.clear(); + group_ptr.push_back(0); + for (auto group = 0ul; group < grouped.getOuterSize(); ++group) { + for (auto k = 0ul; k < grouped.getNumNonZeros(group); ++k) { + auto idx = grouped.getInnerIdx(group, k); + if (static_cast(idx) >= bound) + continue; // prevent out of bounds in LU-SGS kernels if more than 1 mpi task + flat_idx.push_back(static_cast(idx)); + } + group_ptr.push_back(static_cast(flat_idx.size())); } - ilu_color_ptr.push_back(static_cast(color_idx.size())); + return GPUMemoryAllocation::gpu_alloc_cpy(flat_idx.data(), flat_idx.size() * sizeof(su2uint)); + }; + + if (lu_sgs_on_device) { + // get the zero-filled sparse pattern for the LU-SGS + const auto& pat_lusgs = geometry->GetSparsePattern(type, 0); + + /*--- Compute the levels using the lower pattern for the forward pass and + * reverse the levels for the backward pass. This works if L and U are symmetric, to be verified ---*/ + auto levels_lusgs = computeLevels(pat_lusgs.l); + + /*--- Flatten levels_lusgs. It drives both triangular solves on the device. ---*/ + d_precond_level_idx = FlattenGroupToDevice(levels_lusgs, precond_level_ptr, nPointDomain, nPointDomain); } - d_ilu_color_idx = GPUMemoryAllocation::gpu_alloc_cpy(color_idx.data(), color_idx.size() * sizeof(su2uint)); - /*--- Flatten levels_ilu the same way. It drives both triangular solves on the device. ---*/ - std::vector level_idx; - level_idx.reserve(nPointDomain); - ilu_level_ptr.clear(); - ilu_level_ptr.push_back(0); - for (auto level = 0ul; level < levels_ilu.getOuterSize(); ++level) { - for (auto k = 0ul; k < levels_ilu.getNumNonZeros(level); ++k) { - level_idx.push_back(static_cast(levels_ilu.getInnerIdx(level, k))); - } - ilu_level_ptr.push_back(static_cast(level_idx.size())); + if (ilu_needed) { + /*--- The factors are built and used on the device, only the pattern and the level table + * are uploaded (once, here) because they do not change. ---*/ + gpu_ilu.nnz_l = ilu.nnz_l; + gpu_ilu.nnz_u = ilu.nnz_u; + GPUAllocAndInit(gpu_ilu.d, nPointDomain * nVar * nEqn); + GPUAllocAndInit(gpu_ilu.l, ilu.nnz_l * nVar * nEqn); + GPUAllocAndInit(gpu_ilu.u, ilu.nnz_u * nVar * nEqn); + GPUAllocAndCopy(gpu_ilu.row_ptr_l, ilu.row_ptr_l, nPointDomain + 1); + GPUAllocAndCopy(gpu_ilu.col_ind_l, ilu.col_ind_l, ilu.nnz_l); + GPUAllocAndCopy(gpu_ilu.row_ptr_u, ilu.row_ptr_u, nPointDomain + 1); + GPUAllocAndCopy(gpu_ilu.col_ind_u, ilu.col_ind_u, ilu.nnz_u); + + /*--- Flatten the coloring, the index type differs from the one of the pattern. It drives + * the factorization on the device. ---*/ + d_ilu_color_idx = FlattenGroupToDevice(color_ilu, ilu_color_ptr, nPointDomain); + + /*--- Flatten levels_ilu the same way. It drives both triangular solves on the device. ---*/ + d_precond_level_idx = FlattenGroupToDevice(levels_ilu, precond_level_ptr, nPointDomain); } - d_ilu_level_idx = GPUMemoryAllocation::gpu_alloc_cpy(level_idx.data(), level_idx.size() * sizeof(su2uint)); } /*--- Thread parallel initialization. ---*/ @@ -1284,11 +1296,69 @@ void CSysMatrix::ComputeILUPreconditioner(const CSysVector +void CSysMatrix::BuildLU_SGSPreconditioner() { + SU2_ZONE_SCOPED + + /*--- Quantize diagonal blocks if mode is active ---*/ + QuantizeDiagonalBlocks(); + + /*--- if on GPU, precompute the inverse of the diagonal D. Otherwise, this is a no-op ---*/ + if (useCuda) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + SU2_DEVICE_REGION(BuildLU_SGSPreconditionerGPU();) + return; + } else { + GPUNotAvailable(CURRENT_FUNCTION); + } +#else + GPUNotAvailable(CURRENT_FUNCTION); +#endif + } +} + template void CSysMatrix::ComputeLU_SGSPreconditioner(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const { SU2_ZONE_SCOPED + + /*--- First part of the symmetric iteration: (D+L).x* = b ---*/ + ComputeLU_SGSPreconditionerForward(vec, prod); + + /*--- MPI Parallelization ---*/ + + CSysMatrixComms::Initiate(prod, geometry, config); + CSysMatrixComms::Complete(prod, geometry, config); + + /*--- Second part of the symmetric iteration: (D+U).x_(1) = D.x* ---*/ + ComputeLU_SGSPreconditionerBackward(prod); + + /*--- MPI Parallelization ---*/ + + CSysMatrixComms::Initiate(prod, geometry, config); + CSysMatrixComms::Complete(prod, geometry, config); +} + +template +void CSysMatrix::ComputeLU_SGSPreconditionerForward(const CSysVector& vec, + CSysVector& prod) const { + SU2_ZONE_SCOPED + + if (useCuda) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + SU2_DEVICE_REGION(ComputeLU_SGSForwardGPU(vec, prod);) + return; + } else { + GPUNotAvailable(CURRENT_FUNCTION); + } +#else + GPUNotAvailable(CURRENT_FUNCTION); +#endif + } + /*--- First part of the symmetric iteration: (D+L).x* = b ---*/ /*--- Coherent view of vectors. ---*/ @@ -1325,14 +1395,27 @@ void CSysMatrix::ComputeLU_SGSPreconditioner(const CSysVector +void CSysMatrix::ComputeLU_SGSPreconditionerBackward(CSysVector& prod) const { + SU2_ZONE_SCOPED /*--- Second part of the symmetric iteration: (D+U).x_(1) = D.x* ---*/ + if (useCuda) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + SU2_DEVICE_REGION(ComputeLU_SGSBackwardGPU(prod);) + return; + } else { + GPUNotAvailable(CURRENT_FUNCTION); + } +#else + GPUNotAvailable(CURRENT_FUNCTION); +#endif + } + /*--- OpenMP Parallelization ---*/ SU2_OMP_FOR_STAT(1) for (unsigned long thread = 0; thread < omp_num_parts; ++thread) { @@ -1363,11 +1446,6 @@ void CSysMatrix::ComputeLU_SGSPreconditioner(const CSysVector diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 56b8a100631..723c8fb8268 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -308,6 +308,75 @@ __global__ void IluFactorColorKernel(const su2uint* __restrict__ color_idx, unsi M.d[iRow * blockSize + tid] = Lij[tid]; } +/*! + * \brief Compute blk[iVar,jVar].x[col,jVar] and sum over neighbor rows on device + * \note used for L.x* and U.x* in ILU and LU-SGS preconditioners + */ +template +__device__ FORCEINLINE ScalarType DeviceSparseBlockMatVec(unsigned long iRow, unsigned long iVar, unsigned long jVar, unsigned long nVar, + const su2uint* __restrict__ row_ptr, const su2uint* __restrict__ col_ind, + const ScalarType* __restrict__ blk, const ScalarType* __restrict__ x, + unsigned long nRows= ~0ul) { + + const auto blockSize = nVar * nVar; + // compute blk[iVar,jVar].x[col,jVar] and sum over row + ScalarType acc = 0; + for (auto k = row_ptr[iRow]; k < row_ptr[iRow + 1]; ++k) { + const unsigned long jPoint = col_ind[k]; + if (jPoint >= nRows) break; //default is largest possible value thus skipped by default + acc += blk[k * blockSize + iVar * nVar + jVar] * x[jPoint * nVar + jVar]; + } + return acc; +} + +/*! + * \brief Compute Quantized blk[iVar,jVar].x[col,jVar] and sum over neighbor rows on device + */ +template +__device__ FORCEINLINE ScalarType QuantizedDeviceSparseBlockMatVec(unsigned long iRow, unsigned long iVar, unsigned long jVar, unsigned long nVar, + const su2uint* __restrict__ row_ptr, const su2uint* __restrict__ col_ind, + const QuantType* __restrict__ q_blk, const QuantScaleType* __restrict__ q_scale, + const ScalarType* __restrict__ x, unsigned long nRows= ~0ul) { + + const auto blockSize = nVar * nVar; + ScalarType acc = 0; + + for (auto k = row_ptr[iRow]; k < row_ptr[iRow + 1]; ++k) { + const unsigned long jPoint = col_ind[k]; + if (jPoint >= nRows) break; //default is largest possible value thus skipped by default + const float scale = DecodeQuantScale(q_scale[k * nVar + iVar]); + ScalarType q_val = static_cast(q_blk[k * blockSize + iVar * nVar + jVar]); // directly cast to ScalarType + acc += scale * q_val * x[jPoint * nVar + jVar]; + } + return acc; +} + + +/*! + * \brief Compute the partial sum across a row on device + * \note used after DeviceSparseBlockMatVec, it completes the dot product for a given iVar + */ +template +__device__ FORCEINLINE ScalarType DeviceReduceBlockRow(const ScalarType* __restrict__ x, unsigned long iVar, unsigned long nVar) { + ScalarType sum = 0; + for (auto j = 0ul; j < nVar; ++j) sum += x[iVar * nVar + j]; + return sum; +} + +/*! + * \brief Compute the block by vector multiplication + */ +template + __device__ FORCEINLINE ScalarType DeviceDenseBlockMatVec(const ScalarType* __restrict__ blk, const ScalarType* __restrict__ x, + ScalarType* __restrict__ partial, unsigned long tid, + unsigned long iVar, unsigned long jVar, unsigned long nVar) { + + // Compute blk.x + partial[tid] = blk[iVar * nVar + jVar] * x[jVar]; + __syncthreads(); + return DeviceReduceBlockRow(partial, iVar, nVar); +} + /*! * \brief Exact forward substitution for the rows of one level, (L+I).prod = vec. * \note Every row in a level only depends on rows in earlier levels, which are already @@ -333,20 +402,10 @@ __global__ void IluForwardKernel(const su2uint* __restrict__ level_idx, unsigned extern __shared__ __align__(sizeof(double)) char smem[]; auto* partial = reinterpret_cast(smem); - ScalarType acc = 0; - for (auto kl = M.row_ptr_l[iRow]; kl < M.row_ptr_l[iRow + 1]; ++kl) { - const unsigned long jPoint = M.col_ind_l[kl]; - const auto* blk = M.l + kl * nVar * nVar; - acc += blk[iVar * nVar + jVar] * prod[jPoint * nVar + jVar]; - } - partial[tid] = acc; + partial[tid] = DeviceSparseBlockMatVec(iRow, iVar, jVar, nVar, M.row_ptr_l, M.col_ind_l, M.l, prod); __syncthreads(); - if (jVar == 0) { - ScalarType sum = vec[iRow * nVar + iVar]; - for (auto j = 0ul; j < nVar; ++j) sum -= partial[iVar * nVar + j]; - prod[iRow * nVar + iVar] = sum; - } + if (jVar == 0) prod[iRow * nVar + iVar] = vec[iRow * nVar + iVar] - DeviceReduceBlockRow(partial, iVar, nVar); } /*! @@ -375,30 +434,14 @@ __global__ void IluBackwardKernel(const su2uint* __restrict__ level_idx, unsigne auto* partial = reinterpret_cast(smem); auto* aux = partial + blockSize; - ScalarType acc = 0; - for (auto ku = M.row_ptr_u[iRow]; ku < M.row_ptr_u[iRow + 1]; ++ku) { - const unsigned long jPoint = M.col_ind_u[ku]; - if (jPoint >= nRows) break; - const auto* blk = M.u + ku * blockSize; - acc += blk[iVar * nVar + jVar] * prod[jPoint * nVar + jVar]; - } - partial[tid] = acc; + partial[tid] = DeviceSparseBlockMatVec(iRow, iVar, jVar, nVar, M.row_ptr_u, M.col_ind_u, M.u, prod, nRows); __syncthreads(); - if (jVar == 0) { - ScalarType sum = prod[iRow * nVar + iVar]; - for (auto j = 0ul; j < nVar; ++j) sum -= partial[iVar * nVar + j]; - aux[iVar] = sum; - } + if (jVar == 0) aux[iVar] = prod[iRow * nVar + iVar] - DeviceReduceBlockRow(partial, iVar, nVar); __syncthreads(); - if (jVar == 0) { - /*--- The diagonal blocks are stored inverted by the factorization. ---*/ - const auto* invUii = M.d + iRow * blockSize; - ScalarType out = 0; - for (auto k = 0ul; k < nVar; ++k) out += invUii[iVar * nVar + k] * aux[k]; - prod[iRow * nVar + iVar] = out; - } + ScalarType out = DeviceDenseBlockMatVec(M.d + iRow * blockSize, aux, partial, tid, iVar, jVar, nVar); + if (jVar == 0) prod[iRow * nVar + iVar] = out; } /*! @@ -636,7 +679,7 @@ void CSysMatrix::ComputeILUPreconditionerGPU(const CSysVector::ComputeILUPreconditionerGPU(const CSysVector::ComputeILUPreconditionerGPU(const CSysVector - <<>>(d_ilu_level_idx, begin, size, nVar, M, d_vec, d_prod); + <<>>(d_precond_level_idx, begin, size, nVar, M, d_vec, d_prod); } /*--- Backward substitution: one exact pass over the levels in decreasing order, * U.prod = prod, see IluBackwardKernel. ---*/ for (auto level = nLevels; level > 0;) { --level; - const auto begin = ilu_level_ptr[level]; - const auto size = ilu_level_ptr[level + 1] - begin; + const auto begin = precond_level_ptr[level]; + const auto size = precond_level_ptr[level + 1] - begin; if (size == 0) continue; IluBackwardKernel - <<>>(d_ilu_level_idx, begin, size, nPointDomain, nVar, M, d_prod); + <<>>(d_precond_level_idx, begin, size, nPointDomain, nVar, M, d_prod); } gpuErrChk(cudaStreamEndCapture(aux_stream, &graph)); - gpuErrChk(cudaGraphInstantiate(&ilu_apply_graph_exec, graph, nullptr, nullptr, 0)); + gpuErrChk(cudaGraphInstantiate(&precond_fwd_graph_exec, graph, nullptr, nullptr, 0)); gpuErrChk(cudaGraphDestroy(graph)); - ilu_apply_graph_vec = d_vec; - ilu_apply_graph_prod = d_prod; + precond_fwd_graph_vec = d_vec; + precond_fwd_graph_prod = d_prod; } - gpuErrChk(cudaGraphLaunch(ilu_apply_graph_exec, aux_stream)); + gpuErrChk(cudaGraphLaunch(precond_fwd_graph_exec, aux_stream)); gpuErrChk(cudaStreamSynchronize(aux_stream)); gpuErrChk(cudaGetLastError()); } +/*! + * \brief Exact forward substitution for the rows of one level, x* = D^{-1}.(b-Lx*) + * \note See notes in IluForwardKernel for more details. + */ +template +__global__ void LU_SGS_ForwardKernel(const su2uint* __restrict__ level_idx, unsigned long level_begin, + unsigned long level_size, unsigned long nVar, DeviceLDU M, + const QuantType* __restrict__ q_l, const QuantScaleType* __restrict__ q_scale_l, + const ScalarType* __restrict__ invD, const ScalarType* __restrict__ vec, + ScalarType* __restrict__ prod, bool quantized_mode) { + if (blockIdx.x >= level_size) return; + + const unsigned long iRow = level_idx[level_begin + blockIdx.x]; + const auto blockSize = nVar * nVar; + const unsigned long tid = threadIdx.x; + const auto iVar = tid / nVar, jVar = tid % nVar; + + extern __shared__ __align__(sizeof(double)) char smem[]; + auto* partial = reinterpret_cast(smem); // serves nVar * nVar threads + auto* aux = partial + blockSize; // skip nVar * nVar threads, serves nVar threads + + // Compute L.x* + if (quantized_mode) { + partial[tid] = QuantizedDeviceSparseBlockMatVec(iRow, iVar, jVar, nVar, M.row_ptr_l, M.col_ind_l, q_l, q_scale_l, prod); + } else { + partial[tid] = DeviceSparseBlockMatVec(iRow, iVar, jVar, nVar, M.row_ptr_l, M.col_ind_l, M.l, prod); + } + __syncthreads(); + + // Compute y = b - L.x* + if (jVar == 0) aux[iVar] = vec[iRow * nVar + iVar] - DeviceReduceBlockRow(partial, iVar, nVar); + __syncthreads(); + + // Compute x* - D^{-1}.y + ScalarType out = DeviceDenseBlockMatVec(invD + iRow * blockSize, aux, partial, tid, iVar, jVar, nVar); + if (jVar == 0) prod[iRow * nVar + iVar] = out; +} + + +/*! + * \brief Exact backward substitution for the rows of one level, x* = D^{-1}.(D.x* - U.x) = x* - D^{-1}.U.x + * \note See notes in IluBackwardKernel for more details + */ +template +__global__ void LU_SGS_BackwardKernel(const su2uint* __restrict__ level_idx, unsigned long level_begin, + unsigned long level_size, unsigned long nRows, unsigned long nVar, + DeviceLDU M, const QuantType* __restrict__ q_u, + const QuantScaleType* __restrict__ q_scale_u, const ScalarType* __restrict__ invD, + ScalarType* __restrict__ prod, bool quantized_mode) { + if (blockIdx.x >= level_size) return; + + const unsigned long iRow = level_idx[level_begin + blockIdx.x]; + const auto blockSize = nVar * nVar; + const unsigned long tid = threadIdx.x; + const auto iVar = tid / nVar, jVar = tid % nVar; + + extern __shared__ __align__(sizeof(double)) char smem[]; + auto* partial = reinterpret_cast(smem); // serves nVar * nVar threads + auto* aux = partial + blockSize; // skip nVar * nVar threads, serves nVar threads + + // Compute U.x + if (quantized_mode) { + partial[tid] = QuantizedDeviceSparseBlockMatVec(iRow, iVar, jVar, nVar, M.row_ptr_u, M.col_ind_u, q_u, q_scale_u, prod, nRows); + } else { + partial[tid] = DeviceSparseBlockMatVec(iRow, iVar, jVar, nVar, M.row_ptr_u, M.col_ind_u, M.u, prod, nRows); + } + __syncthreads(); + + + if (jVar == 0) aux[iVar] = DeviceReduceBlockRow(partial, iVar, nVar); + __syncthreads(); + + // Compute x* - D^{-1}.(U.x) + ScalarType correction = DeviceDenseBlockMatVec(invD + iRow * blockSize, aux, partial, tid, iVar, jVar, nVar); + if (jVar == 0) prod[iRow * nVar + iVar] -= correction; + +} + +/*! + * \brief Pre-calculates the inverse of the diagonal matrix D, same as for the Jacobi preconditioner + */ +template +void CSysMatrix::BuildLU_SGSPreconditionerGPU() { + SU2_ZONE_SCOPED + if (d_invM == nullptr) { + SU2_MPI::Error("CUDA LU-SGS preconditioner used without device storage.", CURRENT_FUNCTION); + } + if (nPointDomain == 0) return; + + /*--- The matrix is expected to be on the device already, it is uploaded once per solve by + * CSysMatrixVectorProduct, which is created before the preconditioner is built. ---*/ + const auto blockSize = static_cast(nVar * nVar); + InvertDiagonalBlocksKernel<<(nPointDomain), blockSize, + 2 * blockSize * sizeof(ScalarType)>>>(nPointDomain, nVar, gpu.d, d_invM); + /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ + gpuErrChk(cudaStreamSynchronize(nullptr)); + gpuErrChk(cudaGetLastError()); +} + +/*! + * \brief Compute the LU-SGS preconditioner forward pass + */ +template +void CSysMatrix::ComputeLU_SGSForwardGPU(const CSysVector& vec, + CSysVector& prod) const { + SU2_ZONE_SCOPED + + if (d_invM == nullptr) { + SU2_MPI::Error("CUDA LU-SGS preconditioner used without device storage.", CURRENT_FUNCTION); + } + if (nPointDomain == 0) return; + + const DeviceLDU M{gpu.d, gpu.l, gpu.u, gpu.row_ptr_l, + gpu.col_ind_l, gpu.row_ptr_u, gpu.col_ind_u}; + + auto* d_vec = vec.GetDevicePointer(); + auto* d_prod = prod.GetDevicePointer(); + + /*--- One thread per block entry, as done in ILU preconditioner ---*/ + const auto threads = static_cast(nVar * nVar); + const auto sharedForward = (threads + nVar) * sizeof(ScalarType); + + if (aux_stream == nullptr) gpuErrChk(cudaStreamCreate(&aux_stream)); + + /*--- First part of the symmetric iteration: (D+L).x* = b ---*/ + if (precond_fwd_graph_exec == nullptr || precond_fwd_graph_vec != d_vec || precond_fwd_graph_prod != d_prod) { + if (precond_fwd_graph_exec != nullptr) { + gpuErrChk(cudaGraphExecDestroy(precond_fwd_graph_exec)); + precond_fwd_graph_exec = nullptr; + } + + cudaGraph_t graph; + gpuErrChk(cudaStreamBeginCapture(aux_stream, cudaStreamCaptureModeThreadLocal)); + + const auto nLevels = precond_level_ptr.size() - 1; + /*--- Forward substitution: compute x* = D^{-1}.(vec - L.x*) ---*/ + for (auto level = 0ul; level < nLevels; ++level) { + const auto begin = precond_level_ptr[level]; + const auto size = precond_level_ptr[level + 1] - begin; + if (size == 0) continue; + LU_SGS_ForwardKernel<<>>(d_precond_level_idx, begin, size, nVar, M, d_q_blocks.l, d_q_scale.l, d_invM, d_vec, d_prod, quantized_mode); + } + + gpuErrChk(cudaStreamEndCapture(aux_stream, &graph)); + gpuErrChk(cudaGraphInstantiate(&precond_fwd_graph_exec, graph, nullptr, nullptr, 0)); + gpuErrChk(cudaGraphDestroy(graph)); + precond_fwd_graph_vec = d_vec; + precond_fwd_graph_prod = d_prod; + + } + + gpuErrChk(cudaGraphLaunch(precond_fwd_graph_exec, aux_stream)); + gpuErrChk(cudaStreamSynchronize(aux_stream)); + gpuErrChk(cudaGetLastError()); + +} + +/*! + * \brief Compute the LU-SGS preconditioner forward pass + */ +template +void CSysMatrix::ComputeLU_SGSBackwardGPU(CSysVector& prod) const { + SU2_ZONE_SCOPED + + if (d_invM == nullptr) { + SU2_MPI::Error("CUDA LU-SGS preconditioner used without device storage.", CURRENT_FUNCTION); + } + if (nPointDomain == 0) return; + + const DeviceLDU M{gpu.d, gpu.l, gpu.u, gpu.row_ptr_l, + gpu.col_ind_l, gpu.row_ptr_u, gpu.col_ind_u}; + + auto* d_prod = prod.GetDevicePointer(); + + /*--- One thread per block entry, as done in ILU preconditioner ---*/ + const auto threads = static_cast(nVar * nVar); + const auto sharedBackward = (threads + nVar) * sizeof(ScalarType); + + if (aux_stream == nullptr) gpuErrChk(cudaStreamCreate(&aux_stream)); + + /*--- Second part of the symmetric iteration: (D+U).x_(1) = D.x* ---*/ + if (precond_bwd_graph_exec == nullptr || precond_bwd_graph_prod != d_prod) { + if (precond_bwd_graph_exec != nullptr) { + gpuErrChk(cudaGraphExecDestroy(precond_bwd_graph_exec)); + precond_bwd_graph_exec = nullptr; + } + + cudaGraph_t graph; + gpuErrChk(cudaStreamBeginCapture(aux_stream, cudaStreamCaptureModeThreadLocal)); + + const auto nLevels = precond_level_ptr.size() - 1; + /*--- Backward substitution: compute x* = D^{-1}.(D.x* - U.x) = x* - D^{-1}.U.x ---*/ + for (auto level = nLevels; level > 0;) { + --level; + const auto begin = precond_level_ptr[level]; + const auto size = precond_level_ptr[level + 1] - begin; + if (size == 0) continue; + LU_SGS_BackwardKernel<<>>(d_precond_level_idx, begin, size, nPointDomain, nVar, M, d_q_blocks.u, d_q_scale.u, d_invM, d_prod, quantized_mode); + } + + gpuErrChk(cudaStreamEndCapture(aux_stream, &graph)); + gpuErrChk(cudaGraphInstantiate(&precond_bwd_graph_exec, graph, nullptr, nullptr, 0)); + gpuErrChk(cudaGraphDestroy(graph)); + precond_bwd_graph_prod = d_prod; + + } + + gpuErrChk(cudaGraphLaunch(precond_bwd_graph_exec, aux_stream)); + gpuErrChk(cudaStreamSynchronize(aux_stream)); + gpuErrChk(cudaGetLastError()); + +} + template void CSysMatrix::HtDTransfer(bool trigger) const { SU2_ZONE_SCOPED @@ -765,12 +1021,16 @@ template void CSysMatrix::MatrixVectorProductGPU(const CSysVector& v template void CSysMatrix::QuantizeDiagonalBlocksGPU(); \ template void CSysMatrix::BuildJacobiPreconditionerGPU(); \ template void CSysMatrix::BuildILUPreconditionerGPU(); \ +template void CSysMatrix::BuildLU_SGSPreconditionerGPU(); \ template void CSysMatrix::ComputeILUPreconditionerGPU(const CSysVector& vec, \ CSysVector& prod) const; \ template void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector& vec, \ CSysVector& prod, \ CGeometry* geometry, \ - const CConfig* config) const; + const CConfig* config) const;\ +template void CSysMatrix::ComputeLU_SGSForwardGPU(const CSysVector& vec, \ + CSysVector& prod) const; \ +template void CSysMatrix::ComputeLU_SGSBackwardGPU(CSysVector& prod) const; INSTANTIATE_MATRIX(su2mixedfloat) #if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION)