From 7b5ce77941a72d87d115f966dda1bdb2e394524e Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 10 Aug 2026 21:58:30 -0700 Subject: [PATCH 01/30] message format --- Common/src/linear_algebra/CSysMatrix.cpp | 4 ++-- SU2_CFD/include/solvers/CFVMFlowSolverBase.inl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index c9668aa0e88..fff62a363be 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -350,9 +350,9 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi static bool printed = false; if (rank == MASTER_NODE && !printed) { - cout << "GPU ILU scheduling (worst rank): " << nColorsMax << " colors for the factorization (~" + cout << "GPU ILU scheduling (worst rank):\n " << nColorsMax << " colors for the factorization (~" << minAvgColorSize << " points/color on average),\n" - << " " << nLevelsMax << " levels for the triangular solves (~" << minAvgLevelSize + << " " << nLevelsMax << " levels for the triangular solves (~" << minAvgLevelSize << " points/level on average)." << endl; printed = true; } diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 30dc89f0a43..5b1c49b75d9 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -362,9 +362,9 @@ void CFVMFlowSolverBase::HybridParallelInitialization(const CConfig& confi if (SU2_MPI::GetRank() == MASTER_NODE && numRanksUsingReducer != SU2_MPI::GetSize()) { cout << "Among the ranks that use edge coloring,\n" - << " the minimum efficiency is " << minColoredParallelEff << ",\n" - << " the maximum number of colors is " << maxColoredNumColors << ",\n" - << " the minimum edge color group size is " << minColoredEdgeColorGroupSize << "." << endl; + << " the minimum efficiency is " << minColoredParallelEff << ",\n" + << " the maximum number of colors is " << maxColoredNumColors << ",\n" + << " the minimum edge color group size is " << minColoredEdgeColorGroupSize << "." << endl; } } From 4cb8044a606c494fcfcc1068dfdfd60a700956ab Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 11 Aug 2026 19:38:39 -0700 Subject: [PATCH 02/30] Add quantized Jacobi and identity preconditioners Extends the Q_LU_SGS quantization scheme to JACOBI and IDENTITY: the matrix-vector product used by the Krylov solver gets int8-quantized off-diagonal storage, while each preconditioner's own apply step stays at full precision (Jacobi's inverse diagonal is always computed from the full-precision diagonal; identity is a no-op either way). Co-Authored-By: Claude Sonnet 5 --- .../linear_algebra/CPreconditioner.hpp | 68 +++++++++++++++++++ Common/include/linear_algebra/CSysMatrix.hpp | 8 ++- Common/include/option_structure.hpp | 7 ++ Common/src/CConfig.cpp | 2 + Common/src/linear_algebra/CSysMatrix.cpp | 20 +++--- Common/src/linear_algebra/CSysSolve.cpp | 7 ++ config_template.cfg | 12 ++-- 7 files changed, 108 insertions(+), 16 deletions(-) diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index e23f5de381f..e1c9a27718d 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -116,6 +116,30 @@ class CIdentityPreconditioner final : public CPreconditioner { inline bool IsIdentity() const override { return true; } }; +/*! + * \class CQuantizedIdentityPreconditioner + * \brief No-op preconditioner, same as CIdentityPreconditioner, but requests quantized (int8) + * off-diagonal storage for the matrix-vector product shared with the Krylov solver, as Q_LU_SGS/Q_JACOBI. + */ +template +class CQuantizedIdentityPreconditioner final : public CPreconditioner { + private: + CSysMatrix& sparse_matrix; + + public: + inline explicit CQuantizedIdentityPreconditioner(CSysMatrix& matrix_ref) : sparse_matrix(matrix_ref) {} + + CQuantizedIdentityPreconditioner() = delete; + + inline void operator()(const CSysVector& u, CSysVector& v) const override { v = u; } + + inline bool IsIdentity() const override { return true; } + + /*! \brief Quantize the diagonal blocks (off diagonals are quantized on the fly), needed by the + * quantized matrix-vector product even though this preconditioner itself is a no-op. */ + inline void Build() override { sparse_matrix.QuantizeDiagonalBlocks(); } +}; + /*! * \class CJacobiPreconditioner * \brief Specialization of preconditioner that uses CSysMatrix class. @@ -162,6 +186,44 @@ class CJacobiPreconditioner final : public CPreconditioner { inline void Build() override { sparse_matrix.BuildJacobiPreconditioner(); } }; +/*! + * \class CQuantizedJacobiPreconditioner + * \brief Specialization of preconditioner that uses CSysMatrix class. + * \note The preconditioner operation itself is identical to CJacobiPreconditioner (it only ever + * applies the full precision inverse diagonal); quantization here only affects the off-diagonal + * blocks used by the matrix-vector product shared with the Krylov solver, exactly as Q_LU_SGS. + */ +template +class CQuantizedJacobiPreconditioner final : public CPreconditioner { + private: + CSysMatrix& sparse_matrix; + CGeometry* geometry; + const CConfig* config; + + public: + inline CQuantizedJacobiPreconditioner(CSysMatrix& matrix_ref, CGeometry* geometry_ref, + const CConfig* config_ref) + : sparse_matrix(matrix_ref) { + if ((geometry_ref == nullptr) || (config_ref == nullptr)) + SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION); + geometry = geometry_ref; + config = config_ref; + } + + CQuantizedJacobiPreconditioner() = delete; + + inline void operator()(const CSysVector& u, CSysVector& v) const override { + sparse_matrix.ComputeJacobiPreconditioner(u, v, geometry, config); + } + + /*! \brief Build the (full precision) inverse diagonal, then quantize the diagonal blocks so the + * quantized matrix-vector product (off diagonals are quantized on the fly) is consistent. */ + inline void Build() override { + sparse_matrix.BuildJacobiPreconditioner(); + sparse_matrix.QuantizeDiagonalBlocks(); + } +}; + /*! * \class CILUPreconditioner * \brief Specialization of preconditioner that uses CSysMatrix class @@ -406,9 +468,15 @@ CPreconditioner* CPreconditioner::Create(ENUM_LINEAR_SOL case IDENTITY: prec = new CIdentityPreconditioner(); break; + case Q_IDENTITY: + prec = new CQuantizedIdentityPreconditioner(jacobian); + break; case JACOBI: prec = new CJacobiPreconditioner(jacobian, geometry, config); break; + case Q_JACOBI: + prec = new CQuantizedJacobiPreconditioner(jacobian, geometry, config); + break; case LINELET: prec = new CLineletPreconditioner(jacobian, geometry, config); break; diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 38f3c902214..9133aaabbb8 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -256,9 +256,13 @@ class CSysMatrix { /*--- Quantized off-diagonal storage (used when quantized_mode == true). ---*/ using QuantType = int8_t; - /*! \brief Set by Initialize() when preconditioner == Q_LU_SGS. + /*! \brief Set by Initialize() when preconditioner == Q_LU_SGS, Q_JACOBI or Q_IDENTITY. * mat.l and mat.u are NOT allocated; off-diagonal blocks live in the - * q_* arrays below. */ + * q_* arrays below. Only the matrix-vector product (used by the Krylov + * solver and, for Q_LU_SGS, by the sweeps) reads the quantized blocks; + * the Jacobi preconditioner never touches them since it only applies the + * (full precision) inverse diagonal, and the identity preconditioner does + * not touch the matrix at all. */ #ifndef CODI_REVERSE_TYPE bool quantized_mode = false; #else diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 2dcd241f9de..c6bddc15aa6 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2529,17 +2529,24 @@ enum ENUM_LINEAR_SOLVER_PREC { LINELET, /*!< \brief Line implicit preconditioner. */ ILU, /*!< \brief ILU(k) preconditioner. */ Q_LU_SGS, /*!< \brief LU-SGS with quantized (int8) off-diagonal storage; L/U are never allocated as ScalarType. */ + Q_JACOBI, /*!< \brief Jacobi with quantized (int8) off-diagonal storage; same matvec quantization as Q_LU_SGS, + the diagonal inverse is still computed and applied at full precision. */ + Q_IDENTITY, /*!< \brief No preconditioner, but the matrix-vector product still uses quantized (int8) + off-diagonal storage, same matvec quantization as Q_LU_SGS/Q_JACOBI. */ PASTIX_ILU=10, /*!< \brief PaStiX ILU(k) preconditioner. */ PASTIX_LU_P, /*!< \brief PaStiX LU as preconditioner. */ PASTIX_LDLT_P, /*!< \brief PaStiX LDLT as preconditioner. */ }; static const MapType Linear_Solver_Prec_Map = { MakePair("NONE", IDENTITY) + MakePair("IDENTITY", IDENTITY) MakePair("JACOBI", JACOBI) MakePair("LU_SGS", LU_SGS) MakePair("LINELET", LINELET) MakePair("ILU", ILU) MakePair("Q_LU_SGS", Q_LU_SGS) + MakePair("Q_JACOBI", Q_JACOBI) + MakePair("Q_IDENTITY", Q_IDENTITY) MakePair("PASTIX_ILU", PASTIX_ILU) MakePair("PASTIX_LU", PASTIX_LU_P) MakePair("PASTIX_LDLT", PASTIX_LDLT_P) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index e075375001b..e3786124458 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -7473,6 +7473,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { case LU_SGS: cout << "Using LU-SGS preconditioning."<< endl; break; case Q_LU_SGS: cout << "Using LU-SGS preconditioning with matrix quantization."<< endl; break; case JACOBI: cout << "Using Jacobi preconditioning."<< endl; break; + case Q_JACOBI: cout << "Using Jacobi preconditioning with matrix quantization."<< endl; break; } break; case SMOOTHER: @@ -7482,6 +7483,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { case LU_SGS: cout << "A LU-SGS"; break; case Q_LU_SGS: cout << "A quantized LU-SGS"; break; case JACOBI: cout << "A Jacobi"; break; + case Q_JACOBI: cout << "A quantized Jacobi"; break; } cout << " method is used for smoothing the linear system." << endl; break; diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index fff62a363be..df6fa449540 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -207,16 +207,17 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi useCuda = config->GetCUDA(); const bool ilu_needed = (prec == ILU); - const bool diag_needed = (prec == JACOBI) || (prec == LINELET); + const bool diag_needed = (prec == JACOBI) || (prec == Q_JACOBI) || (prec == LINELET); /*--- Linelet also builds the Jacobi preconditioner but reads the inverse diagonal blocks on - * the host, so only plain Jacobi can keep them exclusively on the device. ---*/ - jacobi_on_device = useCuda && (prec == JACOBI); + * the host, so only plain (or quantized) Jacobi can keep them exclusively on the device. ---*/ + jacobi_on_device = useCuda && (prec == JACOBI || prec == Q_JACOBI); #ifndef CODI_REVERSE_TYPE - const bool q_lus_needed = allow_quant && !useCuda && (prec == Q_LU_SGS); + const bool quantized_offdiag_needed = + allow_quant && !useCuda && (prec == Q_LU_SGS || prec == Q_JACOBI || prec == Q_IDENTITY); #else /*--- No quantization in adjoint mode for now because TransposeInPlace would get complicated. ---*/ - const bool q_lus_needed = false; + const bool quantized_offdiag_needed = false; #endif /*--- Basic dimensions. ---*/ @@ -242,9 +243,9 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi } allocAndInit(mat.d, nPoint * nVar * nEqn); - if (q_lus_needed) { - /*--- Q_LU_SGS: no full-precision L/U; off-diagonal blocks live in quantized storage. - * L/U are quantized on-the-fly during assembly; diagonal is quantized in Build step. ---*/ + if (quantized_offdiag_needed) { + /*--- Q_LU_SGS / Q_JACOBI / Q_IDENTITY: no full-precision L/U; off-diagonal blocks live in quantized storage. + * L/U are quantized on-the-fly during assembly; diagonal is quantized in the Build step. ---*/ #ifndef CODI_REVERSE_TYPE quantized_mode = true; #endif @@ -732,7 +733,8 @@ void CSysMatrix::QuantizeDiagonalBlocks() { SU2_ZONE_SCOPED if (quantized_mode) { - /*--- Q_LU_SGS: L/U were quantized during assembly; only the diagonal needs quantization now. ---*/ + /*--- Q_LU_SGS / Q_JACOBI / Q_IDENTITY: L/U were quantized during assembly; only the diagonal needs quantization + * now. ---*/ SU2_OMP_FOR_DYN(omp_heavy_size) for (auto i = 0ul; i < nPointDomain; ++i) QuantizeBlock(&mat.d[i * nVar * nVar], &q_scale_d[i * nVar], &q_blocks_d[i * nVar * nVar]); diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index dcd88c452c5..a008a503169 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -1568,6 +1568,13 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con case LINELET: if (RequiresTranspose) Jacobian.BuildJacobiPreconditioner(); break; + case Q_JACOBI: + /*--- Rebuild the full-precision inverse diagonal and re-quantize it, mirroring Build(). ---*/ + if (RequiresTranspose) { + Jacobian.BuildJacobiPreconditioner(); + Jacobian.QuantizeDiagonalBlocks(); + } + break; case LU_SGS: case Q_LU_SGS: /*--- Nothing to build (transpose path not supported for Q_LU_SGS, see CSysMatrix::Initialize). ---*/ diff --git a/config_template.cfg b/config_template.cfg index 922636bc6b4..ff33a439456 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1648,11 +1648,13 @@ DISCADJ_LIN_SOLVER= FGMRES % Use CUDA GPU Acceleration for FGMRES Linear Solver Only ENABLE_CUDA=NO % -% Preconditioner of the Krylov linear solver or type of smoother (ILU, LU_SGS, Q_LU_SGS, LINELET, JACOBI) -% Q_LU_SGS triggers the use of quantization to reduce the size of the sparse matrix (for compressible flow -% the matrix becomes 3x smaller relative to mixed-precision mode). This is only used by the compressible -% and incompressible solvers, others fallback silently to LU_SGS. A suitable nondimensionalization mode -% MUST be used otherwise the solver is very likely to diverge. +% Preconditioner of the Krylov linear solver or type of smoother (ILU, LU_SGS, Q_LU_SGS, LINELET, JACOBI, Q_JACOBI) +% Q_LU_SGS and Q_JACOBI trigger the use of quantization to reduce the size of the sparse matrix (for +% compressible flow the matrix becomes 3x smaller relative to mixed-precision mode). Q_IDENTITY requests the +% same quantized matrix-vector product without any actual preconditioning (equivalent to NONE otherwise). +% This is only used by the compressible and incompressible solvers, others fallback silently to +% LU_SGS/JACOBI/NONE. A suitable nondimensionalization mode MUST be used otherwise the solver is very +% likely to diverge. LINEAR_SOLVER_PREC= ILU % % Same for discrete adjoint (JACOBI or ILU), replaces LINEAR_SOLVER_PREC in SU2_*_AD codes. From 31b6b8ab47128e9b92990d4ae9c6e78d1ab812b3 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 11 Aug 2026 21:09:32 -0700 Subject: [PATCH 03/30] Enable quantized matrix-vector product on GPU for Jacobi/identity Q_JACOBI and Q_IDENTITY no longer force quantized_mode off under CUDA: their shared matrix-vector product has no sequential dependency (unlike Q_LU_SGS's sweeps, which stay host-only), so it can run on the device. Adds device mirrors of the quantized off-diagonal storage (refreshed by HtDTransfer) and a new QuantizedBlockLDU_SpMV_kernel mirroring the host QuantizedRowProduct. The quantized diagonal is populated by a new QuantizeDiagonalBlocksGPU() kernel that reads directly from gpu.d rather than transferring the host-quantized diagonal, since diagonal quantization only happens in Build(), which runs after HtDTransfer has already uploaded the matrix for the solve. Validated against the host path on a NACA0012 inviscid case (RTX 4070 Ti): GPU quantized-Jacobi results match CPU quantized-Jacobi exactly, and GPU plain-Jacobi is unaffected. Co-Authored-By: Claude Sonnet 5 --- Common/include/linear_algebra/CSysMatrix.hpp | 26 +++- Common/src/linear_algebra/CSysMatrix.cpp | 62 +++++++-- Common/src/linear_algebra/CSysMatrixGPU.cu | 137 ++++++++++++++++++- 3 files changed, 209 insertions(+), 16 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 9133aaabbb8..9dd557c392a 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -273,9 +273,25 @@ class CSysMatrix { QuantType* q_scale_u; /*!< \brief Same as q_scale_l for the upper entries. */ QuantType* q_blocks_u; /*!< \brief Same as q_blocks_l for the upper entries. */ QuantType* q_scale_d; /*!< \brief Same as q_scale_l for the diagonal entries, [nPoint * nVar]. - * Populated by QuantizeDiagonalBlocks(). */ + * Populated by QuantizeDiagonalBlocks() (host only, see below). */ QuantType* q_blocks_d; /*!< \brief Same as q_blocks_l for the diagonal entries. */ + /*--- Device mirrors of the quantized off-diagonal storage, only allocated when + * quantized_mode && useCuda (currently only reachable for Q_JACOBI/Q_IDENTITY, Q_LU_SGS stays + * host-only). d_q_scale_l/d_q_blocks_l/d_q_scale_u/d_q_blocks_u are plain device-side copies of + * q_scale_l/q_blocks_l/q_scale_u/q_blocks_u, refreshed by HtDTransfer(). The diagonal mirrors + * d_q_scale_d/d_q_blocks_d are populated directly on the device from gpu.d by + * QuantizeDiagonalBlocksGPU() instead: q_scale_d/q_blocks_d (host) are only ever quantized from + * mat.d inside QuantizeDiagonalBlocks() itself, which for Q_JACOBI/Q_IDENTITY runs after + * HtDTransfer() has already copied the (unquantized) diagonal, so transferring the host + * quantized diagonal would race the one point in the solve where it is actually computed. ---*/ + QuantType* d_q_scale_l = nullptr; + QuantType* d_q_blocks_l = nullptr; + QuantType* d_q_scale_u = nullptr; + QuantType* d_q_blocks_u = nullptr; + QuantType* d_q_scale_d = nullptr; + QuantType* d_q_blocks_d = nullptr; + bool useCuda = false; /*!< \brief Whether CUDA is enabled. */ /*!< \brief Whether the inverse diagonal blocks are only needed on the device. False for the @@ -558,6 +574,14 @@ class CSysMatrix { void MatrixVectorProductGPU(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const; + /*! + * \brief Quantize the diagonal blocks directly on the device, from gpu.d into + * d_q_scale_d/d_q_blocks_d. Device counterpart of the host branch of + * QuantizeDiagonalBlocks(), only reachable when quantized_mode && useCuda. + * \note Requires the device matrix to be up to date, see HtDTransfer. + */ + void QuantizeDiagonalBlocksGPU(); + /*! * \brief Build the Jacobi preconditioner on the device, from the device copy of the matrix. * \note Requires the device matrix to be up to date, see HtDTransfer. diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index df6fa449540..6de7de4178a 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -152,6 +152,12 @@ CSysMatrix::~CSysMatrix() { GPUMemoryAllocation::gpu_free(d_invM); GPUMemoryAllocation::gpu_free(d_ilu_color_idx); GPUMemoryAllocation::gpu_free(d_ilu_level_idx); + GPUMemoryAllocation::gpu_free(d_q_scale_l); + GPUMemoryAllocation::gpu_free(d_q_blocks_l); + GPUMemoryAllocation::gpu_free(d_q_scale_u); + GPUMemoryAllocation::gpu_free(d_q_blocks_u); + GPUMemoryAllocation::gpu_free(d_q_scale_d); + GPUMemoryAllocation::gpu_free(d_q_blocks_d); #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); @@ -213,8 +219,12 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi * the host, so only plain (or quantized) Jacobi can keep them exclusively on the device. ---*/ jacobi_on_device = useCuda && (prec == JACOBI || prec == Q_JACOBI); #ifndef CODI_REVERSE_TYPE + /*--- Q_LU_SGS stays host-only: its forward/backward sweeps are inherently sequential and + * already bracketed to the host via ApplyPreconditionerOnHost, so there is nothing to gain + * from a device SpMV there. Q_JACOBI and Q_IDENTITY have no such sequential preconditioner + * step, so their (shared) quantized matrix-vector product is allowed on the device too. ---*/ const bool quantized_offdiag_needed = - allow_quant && !useCuda && (prec == Q_LU_SGS || prec == Q_JACOBI || prec == Q_IDENTITY); + allow_quant && (prec == Q_JACOBI || prec == Q_IDENTITY || (prec == Q_LU_SGS && !useCuda)); #else /*--- No quantization in adjoint mode for now because TransposeInPlace would get complicated. ---*/ const bool quantized_offdiag_needed = false; @@ -275,12 +285,29 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi SU2_MPI::Error("CUDA CSysMatrix block-LDU SpMV requires square blocks.", CURRENT_FUNCTION); } GPUAllocAndInit(gpu.d, nPoint * nVar * nEqn); - GPUAllocAndInit(gpu.l, mat.nnz_l * nVar * nEqn); - GPUAllocAndInit(gpu.u, mat.nnz_u * nVar * nEqn); GPUAllocAndCopy(gpu.row_ptr_l, mat.row_ptr_l, nPointDomain + 1); GPUAllocAndCopy(gpu.col_ind_l, mat.col_ind_l, mat.nnz_l); GPUAllocAndCopy(gpu.row_ptr_u, mat.row_ptr_u, nPointDomain + 1); GPUAllocAndCopy(gpu.col_ind_u, mat.col_ind_u, mat.nnz_u); + + if (quantized_mode) { + /*--- Device mirrors of the host quantized off-diagonal storage; gpu.l/gpu.u are not + * allocated (nothing would ever read them). The diagonal is quantized straight from + * gpu.d by QuantizeDiagonalBlocksGPU(), so there is no d_q_scale_d/d_q_blocks_d transfer + * to arrange, only the allocation here. ---*/ + auto GPUAllocQ = [](QuantType*& ptr, unsigned long n) { + ptr = GPUMemoryAllocation::gpu_alloc(n * sizeof(QuantType)); + }; + GPUAllocQ(d_q_scale_l, mat.nnz_l * nVar); + GPUAllocQ(d_q_blocks_l, mat.nnz_l * nVar * nEqn); + GPUAllocQ(d_q_scale_u, mat.nnz_u * nVar); + GPUAllocQ(d_q_blocks_u, mat.nnz_u * nVar * nEqn); + GPUAllocQ(d_q_scale_d, nPoint * nVar); + GPUAllocQ(d_q_blocks_d, nPoint * nVar * nEqn); + } else { + GPUAllocAndInit(gpu.l, mat.nnz_l * nVar * nEqn); + GPUAllocAndInit(gpu.u, mat.nnz_u * nVar * nEqn); + } } if (type == ConnectivityType::FiniteVolume) { @@ -732,14 +759,29 @@ template void CSysMatrix::QuantizeDiagonalBlocks() { SU2_ZONE_SCOPED - if (quantized_mode) { - /*--- Q_LU_SGS / Q_JACOBI / Q_IDENTITY: L/U were quantized during assembly; only the diagonal needs quantization - * now. ---*/ - SU2_OMP_FOR_DYN(omp_heavy_size) - for (auto i = 0ul; i < nPointDomain; ++i) - QuantizeBlock(&mat.d[i * nVar * nVar], &q_scale_d[i * nVar], &q_blocks_d[i * nVar * nVar]); - END_SU2_OMP_FOR + if (!quantized_mode) return; + + if (useCuda) { +#ifdef SU2_ENABLE_CUDA_KERNELS + if constexpr (su2_gpu_capable_v) { + /*--- Quantize straight from gpu.d, see the d_q_scale_d/d_q_blocks_d comment in + * CSysMatrix.hpp for why this does not go through the host q_scale_d/q_blocks_d. ---*/ + SU2_DEVICE_REGION(QuantizeDiagonalBlocksGPU();) + return; + } else { + GPUNotAvailable(CURRENT_FUNCTION); + } +#else + GPUNotAvailable(CURRENT_FUNCTION); +#endif } + + /*--- Q_LU_SGS / Q_JACOBI / Q_IDENTITY: L/U were quantized during assembly; only the diagonal needs quantization + * now. ---*/ + SU2_OMP_FOR_DYN(omp_heavy_size) + for (auto i = 0ul; i < nPointDomain; ++i) + QuantizeBlock(&mat.d[i * nVar * nVar], &q_scale_d[i * nVar], &q_blocks_d[i * nVar * nVar]); + END_SU2_OMP_FOR } template diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 1cd5de51d5c..af0bc3ed91b 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -353,6 +353,96 @@ __global__ void BlockLDU_SpMV_kernel(unsigned long nRows, unsigned long nVar, y[iRow * nVar + iVar] = sum; } +/*! + * \brief Device counterpart of DecodeQuantScale (CSysMatrix.hpp), bit-identical: reconstructs + * the float row-scale 2^e from its packed int8 binary exponent. + */ +__device__ FORCEINLINE float DecodeQuantScaleDevice(int8_t e) { + const uint32_t bits = static_cast(max(0, static_cast(e) + 127)) << 23; + return __uint_as_float(bits); +} + +/*! + * \brief Device counterpart of EncodeQuantBlock (CSysMatrix.hpp), bit-identical: quantizes one + * nVar x nVar block (row-major, pointed to by \p blk) into per-row int8 storage. One + * thread handles one whole block (i.e. one matrix row). + */ +template +__device__ void EncodeQuantBlockDevice(const ScalarType* __restrict__ blk, int8_t* __restrict__ qs, + int8_t* __restrict__ qv, unsigned long nVar) { + for (auto r = 0ul; r < nVar; ++r) { + constexpr uint32_t eps_bits = 0x34000000u; + uint32_t max_abs_bits = eps_bits; + for (auto c = 0ul; c < nVar; ++c) { + const float fv = static_cast(blk[r * nVar + c]); + max_abs_bits = max(max_abs_bits, __float_as_uint(fv) & 0x7FFFFFFFu); + } + const int e = min(127, max(-128, static_cast(max_abs_bits >> 23) - 133)); + qs[r] = static_cast(e); + const float inv_rscale = __uint_as_float(static_cast(127 - e) << 23); + for (auto c = 0ul; c < nVar; ++c) { + qv[r * nVar + c] = + static_cast(max(-128.f, min(127.f, roundf(static_cast(blk[r * nVar + c]) * inv_rscale)))); + } + } +} + +/*! + * \brief Quantize the diagonal blocks straight from the device diagonal (gpu.d), device version + * of CSysMatrix::QuantizeBlock applied row-by-row (host QuantizeDiagonalBlocks). + * One thread per row. + */ +template +__global__ void QuantizeDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, + const ScalarType* __restrict__ mat_d, int8_t* __restrict__ q_scale_d, + int8_t* __restrict__ q_blocks_d) { + const auto iRow = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (iRow >= nRows) return; + + EncodeQuantBlockDevice(mat_d + iRow * nVar * nVar, q_scale_d + iRow * nVar, q_blocks_d + iRow * nVar * nVar, nVar); +} + +/*! + * \brief Quantized block-LDU SpMV kernel: y[iRow] = (L + D + U) * x per block-row, reading int8 + * row-scaled quantized blocks instead of full precision ones. Device version of + * QuantizedRowProduct/QuantizedMatVecAdd (CSysMatrix.inl). Same launch/thread layout as + * BlockLDU_SpMV_kernel: one CUDA block per block-row, threadIdx.x indexes output variable. + */ +template +__global__ void QuantizedBlockLDU_SpMV_kernel( + unsigned long nRows, unsigned long nVar, const su2uint* __restrict__ row_ptr_l, + const su2uint* __restrict__ col_ind_l, const int8_t* __restrict__ q_scale_l, + const int8_t* __restrict__ q_blocks_l, const int8_t* __restrict__ q_scale_d, + const int8_t* __restrict__ q_blocks_d, const su2uint* __restrict__ row_ptr_u, + const su2uint* __restrict__ col_ind_u, const int8_t* __restrict__ q_scale_u, + const int8_t* __restrict__ q_blocks_u, const ScalarType* __restrict__ x, ScalarType* __restrict__ y) { + const unsigned long iRow = blockIdx.x; + const unsigned long iVar = threadIdx.x; + if (iRow >= nRows || iVar >= nVar) return; + + auto addBlock = [&](const int8_t* __restrict__ qs, const int8_t* __restrict__ qv, const ScalarType* __restrict__ xk) { + const float row_scale = DecodeQuantScaleDevice(qs[iVar]); + ScalarType partial = 0; + for (unsigned long jVar = 0; jVar < nVar; ++jVar) partial += qv[iVar * nVar + jVar] * xk[jVar]; + return static_cast(row_scale) * partial; + }; + + ScalarType sum = 0; + /* Lower */ + for (auto k = row_ptr_l[iRow]; k < row_ptr_l[iRow + 1]; ++k) { + const auto col = col_ind_l[k]; + sum += addBlock(q_scale_l + k * nVar, q_blocks_l + k * nVar * nVar, x + col * nVar); + } + /* Diagonal */ + sum += addBlock(q_scale_d + iRow * nVar, q_blocks_d + iRow * nVar * nVar, x + iRow * nVar); + /* Upper */ + for (auto k = row_ptr_u[iRow]; k < row_ptr_u[iRow + 1]; ++k) { + const auto col = col_ind_u[k]; + sum += addBlock(q_scale_u + k * nVar, q_blocks_u + k * nVar * nVar, x + col * nVar); + } + y[iRow * nVar + iVar] = sum; +} + } // namespace template @@ -397,6 +487,23 @@ void CSysMatrix::BuildJacobiPreconditionerGPU() { gpuErrChk(cudaGetLastError()); } +template +void CSysMatrix::QuantizeDiagonalBlocksGPU() { + SU2_ZONE_SCOPED + + 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. ---*/ + constexpr unsigned threadsPerBlock = 128; + const auto blocks = static_cast((nPointDomain + threadsPerBlock - 1) / threadsPerBlock); + QuantizeDiagonalBlocksKernel + <<>>(nPointDomain, nVar, gpu.d, d_q_scale_d, d_q_blocks_d); + /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ + gpuErrChk(cudaStreamSynchronize(nullptr)); + gpuErrChk(cudaGetLastError()); +} + template void CSysMatrix::BuildILUPreconditionerGPU() { SU2_ZONE_SCOPED @@ -535,8 +642,21 @@ void CSysMatrix::HtDTransfer(bool trigger) const { SU2_ZONE_SCOPED if (!trigger) return; gpuErrChk(cudaMemcpy(gpu.d, mat.d, sizeof(ScalarType) * nPoint * nVar * nEqn, cudaMemcpyHostToDevice)); - gpuErrChk(cudaMemcpy(gpu.l, mat.l, sizeof(ScalarType) * mat.nnz_l * nVar * nEqn, cudaMemcpyHostToDevice)); - gpuErrChk(cudaMemcpy(gpu.u, mat.u, sizeof(ScalarType) * mat.nnz_u * nVar * nEqn, cudaMemcpyHostToDevice)); + if (quantized_mode) { + /*--- No gpu.l/gpu.u to transfer (never allocated); mirror the host quantized off-diagonal + * storage instead. The diagonal mirrors (d_q_scale_d/d_q_blocks_d) are NOT transferred here, + * they are populated straight from gpu.d by QuantizeDiagonalBlocksGPU(), see the comment on + * those members in CSysMatrix.hpp for why. ---*/ + gpuErrChk(cudaMemcpy(d_q_scale_l, q_scale_l, sizeof(QuantType) * mat.nnz_l * nVar, cudaMemcpyHostToDevice)); + gpuErrChk( + cudaMemcpy(d_q_blocks_l, q_blocks_l, sizeof(QuantType) * mat.nnz_l * nVar * nEqn, cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpy(d_q_scale_u, q_scale_u, sizeof(QuantType) * mat.nnz_u * nVar, cudaMemcpyHostToDevice)); + gpuErrChk( + cudaMemcpy(d_q_blocks_u, q_blocks_u, sizeof(QuantType) * mat.nnz_u * nVar * nEqn, cudaMemcpyHostToDevice)); + } else { + gpuErrChk(cudaMemcpy(gpu.l, mat.l, sizeof(ScalarType) * mat.nnz_l * nVar * nEqn, cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpy(gpu.u, mat.u, sizeof(ScalarType) * mat.nnz_u * nVar * nEqn, cudaMemcpyHostToDevice)); + } } template @@ -549,9 +669,15 @@ void CSysMatrix::MatrixVectorProductGPU(const CSysVector 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); + if (quantized_mode) { + QuantizedBlockLDU_SpMV_kernel<<>>( + nPointDomain, nVar, gpu.row_ptr_l, gpu.col_ind_l, d_q_scale_l, d_q_blocks_l, d_q_scale_d, d_q_blocks_d, + gpu.row_ptr_u, gpu.col_ind_u, d_q_scale_u, d_q_blocks_u, d_vec, d_prod); + } else { + 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); + } /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ gpuErrChk(cudaStreamSynchronize(nullptr)); gpuErrChk(cudaGetLastError()); @@ -563,6 +689,7 @@ template void CSysMatrix::MatrixVectorProductGPU(const CSysVector& v CSysVector& prod, \ CGeometry* geometry, \ const CConfig* config) const; \ +template void CSysMatrix::QuantizeDiagonalBlocksGPU(); \ template void CSysMatrix::BuildJacobiPreconditionerGPU(); \ template void CSysMatrix::BuildILUPreconditionerGPU(); \ template void CSysMatrix::ComputeILUPreconditionerGPU(const CSysVector& vec, \ From e16cad76466a021852a77f92e139acec4823e4c9 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 11 Aug 2026 21:40:30 -0700 Subject: [PATCH 04/30] Group quantized storage into LDU and overlap diagonal quantize with L/U upload Templatize the existing LDU struct on its value type and reuse it for the quantized storage: q_scale/q_blocks (host) and d_q_scale/d_q_blocks (device) replace the twelve loose QuantType* members, cutting the constructor/destructor/Initialize boilerplate down to a couple of generic-lambda calls per group. Also change how the quantized diagonal reaches the device: instead of a GPU kernel that quantized straight from gpu.d, the diagonal is always quantized on the host in QuantizeDiagonalBlocks() and then uploaded. HtDTransfer() kicks off the (larger) L/U quantized transfer via cudaMemcpyAsync before Build() is reached, so it can overlap with that host-side diagonal quantization instead of blocking on it first; same default stream ordering keeps any later kernel launch correct without an extra explicit sync. Re-validated against the pre-refactor build on the NACA0012 case (RTX 4070 Ti): CPU/GPU x plain/quantized Jacobi and GPU quantized-identity all reproduce identical CL/CD to the prior commit. Co-Authored-By: Claude Sonnet 5 --- Common/include/linear_algebra/CSysMatrix.hpp | 102 +++++++++---------- Common/include/linear_algebra/CSysMatrix.inl | 16 +-- Common/src/linear_algebra/CSysMatrix.cpp | 95 ++++++++--------- Common/src/linear_algebra/CSysMatrixGPU.cu | 81 ++++----------- 4 files changed, 118 insertions(+), 176 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 9dd557c392a..0ac1d2e7bb9 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -231,14 +231,19 @@ class CSysMatrix { /*! * \brief Aggregates value arrays and sparse-structure pointers for an LDU-partitioned matrix. - * Each CSysMatrix holds three LDU instances: the host matrix (mat), its device copy (gpu), - * and the ILU factorization (ilu). Ownership of the value arrays (d/l/u) and whether - * the pointers address host or device memory is managed by CSysMatrix. + * Each CSysMatrix holds three LDU instances: the host matrix (mat), its + * device copy (gpu), and the ILU factorization (ilu). Ownership of the value arrays + * (d/l/u) and whether the pointers address host or device memory is managed by + * CSysMatrix. Also reused with T = QuantType to group the quantized scale/blocks + * storage (q_scale, q_blocks, d_q_scale, d_q_blocks) the same way; for those the pattern + * fields (row_ptr_l, col_ind_l, row_ptr_u, col_ind_u, nnz_l, nnz_u) are simply left + * unused, since the sparsity pattern is already available from mat/gpu. */ + template struct LDU { - ScalarType* d = nullptr; /*!< \brief Diagonal block values. */ - ScalarType* l = nullptr; /*!< \brief Strictly-lower block values. */ - ScalarType* u = nullptr; /*!< \brief Strictly-upper block values. */ + T* d = nullptr; /*!< \brief Diagonal block values. */ + T* l = nullptr; /*!< \brief Strictly-lower block values. */ + T* u = nullptr; /*!< \brief Strictly-upper block values. */ const su2uint* row_ptr_l = nullptr; /*!< \brief Row pointers for L (geometry-owned or GPU copy). */ const su2uint* col_ind_l = nullptr; /*!< \brief Column indices for L. */ const su2uint* row_ptr_u = nullptr; /*!< \brief Row pointers for U. */ @@ -247,50 +252,40 @@ class CSysMatrix { unsigned long nnz_u = 0; /*!< \brief Number of U nonzeros. */ }; - LDU mat; /*!< \brief Host matrix (values owned via aligned_alloc; pattern from geometry). */ - 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). */ + LDU mat; /*!< \brief Host matrix (values owned via aligned_alloc; pattern from geometry). */ + 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. */ /*--- Quantized off-diagonal storage (used when quantized_mode == true). ---*/ using QuantType = int8_t; /*! \brief Set by Initialize() when preconditioner == Q_LU_SGS, Q_JACOBI or Q_IDENTITY. - * mat.l and mat.u are NOT allocated; off-diagonal blocks live in the - * q_* arrays below. Only the matrix-vector product (used by the Krylov - * solver and, for Q_LU_SGS, by the sweeps) reads the quantized blocks; - * the Jacobi preconditioner never touches them since it only applies the - * (full precision) inverse diagonal, and the identity preconditioner does - * not touch the matrix at all. */ + * mat.l and mat.u are NOT allocated; off-diagonal blocks live in q_scale/q_blocks + * below. Only the matrix-vector product (used by the Krylov solver and, for Q_LU_SGS, + * by the sweeps) reads the quantized blocks; the Jacobi preconditioner never touches + * them since it only applies the (full precision) inverse diagonal, and the identity + * preconditioner does not touch the matrix at all. */ #ifndef CODI_REVERSE_TYPE bool quantized_mode = false; #else static constexpr bool quantized_mode = false; #endif - QuantType* q_scale_l; /*!< \brief Per-row exponent for L blocks, [nnz_l * nVar]. */ - QuantType* q_blocks_l; /*!< \brief Quantized L block entries, [nnz_l * nVar * nEqn]. */ - QuantType* q_scale_u; /*!< \brief Same as q_scale_l for the upper entries. */ - QuantType* q_blocks_u; /*!< \brief Same as q_blocks_l for the upper entries. */ - QuantType* q_scale_d; /*!< \brief Same as q_scale_l for the diagonal entries, [nPoint * nVar]. - * Populated by QuantizeDiagonalBlocks() (host only, see below). */ - QuantType* q_blocks_d; /*!< \brief Same as q_blocks_l for the diagonal entries. */ - - /*--- Device mirrors of the quantized off-diagonal storage, only allocated when - * quantized_mode && useCuda (currently only reachable for Q_JACOBI/Q_IDENTITY, Q_LU_SGS stays - * host-only). d_q_scale_l/d_q_blocks_l/d_q_scale_u/d_q_blocks_u are plain device-side copies of - * q_scale_l/q_blocks_l/q_scale_u/q_blocks_u, refreshed by HtDTransfer(). The diagonal mirrors - * d_q_scale_d/d_q_blocks_d are populated directly on the device from gpu.d by - * QuantizeDiagonalBlocksGPU() instead: q_scale_d/q_blocks_d (host) are only ever quantized from - * mat.d inside QuantizeDiagonalBlocks() itself, which for Q_JACOBI/Q_IDENTITY runs after - * HtDTransfer() has already copied the (unquantized) diagonal, so transferring the host - * quantized diagonal would race the one point in the solve where it is actually computed. ---*/ - QuantType* d_q_scale_l = nullptr; - QuantType* d_q_blocks_l = nullptr; - QuantType* d_q_scale_u = nullptr; - QuantType* d_q_blocks_u = nullptr; - QuantType* d_q_scale_d = nullptr; - QuantType* d_q_blocks_d = nullptr; + /*!< \brief Per-row exponents; .l/.u sized [nnz_l/u * nVar], .d [nPoint * nVar] + * (populated by QuantizeDiagonalBlocks(), always on the host, see below). */ + LDU q_scale; + /*!< \brief Quantized block entries; .l/.u sized [nnz_l/u * nVar * nEqn], .d [nPoint * nVar * nEqn]. */ + LDU q_blocks; + + /*!< \brief Device mirrors of the quantized storage, only allocated when quantized_mode && + * useCuda (currently only reachable for Q_JACOBI/Q_IDENTITY, Q_LU_SGS stays host-only). + * d_q_scale.l/.u and d_q_blocks.l/.u are plain device-side copies of q_scale.l/.u and + * q_blocks.l/.u, uploaded *asynchronously* by HtDTransfer() so that transfer can overlap with + * the host quantizing the diagonal (see QuantizeDiagonalBlocks()); d_q_scale.d/d_q_blocks.d are + * that host result, uploaded once ready by QuantizeDiagonalBlocksGPU(). */ + LDU d_q_scale; + LDU d_q_blocks; bool useCuda = false; /*!< \brief Whether CUDA is enabled. */ @@ -562,7 +557,7 @@ class CSysMatrix { /*! \brief Diagonal product using quantized D (Q_LU_SGS backward sweep). */ inline void QuantizedDiagonalProduct(const CSysVector& vec, unsigned long row_i, ScalarType* prod) const; - /*! \brief Gauss elimination on the quantized diagonal block: decodes q_blocks_d into a local + /*! \brief Gauss elimination on the quantized diagonal block: decodes q_blocks.d into a local * ScalarType buffer and delegates to the scalar GaussElimination overload. */ inline void QuantizedGaussElimination(unsigned long block_i, ScalarType* rhs) const; @@ -575,10 +570,9 @@ class CSysMatrix { const CConfig* config) const; /*! - * \brief Quantize the diagonal blocks directly on the device, from gpu.d into - * d_q_scale_d/d_q_blocks_d. Device counterpart of the host branch of - * QuantizeDiagonalBlocks(), only reachable when quantized_mode && useCuda. - * \note Requires the device matrix to be up to date, see HtDTransfer. + * \brief Upload the already host-quantized diagonal (q_scale.d/q_blocks.d, computed by + * QuantizeDiagonalBlocks() just before calling this) into d_q_scale.d/d_q_blocks.d. + * Only reachable when quantized_mode && useCuda. */ void QuantizeDiagonalBlocksGPU(); @@ -695,10 +689,10 @@ class CSysMatrix { } \ if (block_j < block_i) { \ for (auto k = mat.row_ptr_l[block_i]; k < mat.row_ptr_l[block_i + 1]; ++k) \ - if (mat.col_ind_l[k] == block_j) return {nullptr, &q_scale_l[k * nVar], &q_blocks_l[k * nVar * nVar], nVar}; \ + if (mat.col_ind_l[k] == block_j) return {nullptr, &q_scale.l[k * nVar], &q_blocks.l[k * nVar * nVar], nVar}; \ } else { \ for (auto k = mat.row_ptr_u[block_i]; k < mat.row_ptr_u[block_i + 1]; ++k) \ - if (mat.col_ind_u[k] == block_j) return {nullptr, &q_scale_u[k * nVar], &q_blocks_u[k * nVar * nVar], nVar}; \ + if (mat.col_ind_u[k] == block_j) return {nullptr, &q_scale.u[k * nVar], &q_blocks.u[k * nVar * nVar], nVar}; \ } \ return {} GET_BLOCK_VIEW_IMPL; @@ -825,9 +819,9 @@ class CSysMatrix { bij_buf[offset] = PassiveAssign(block_j[iVar][jVar] * scale); bji_buf[offset] = -PassiveAssign(block_i[iVar][jVar] * scale); } - QuantizeBlock(bij_buf, &q_scale_u[iEdge * nVar], &q_blocks_u[iEdge * blkSz]); + QuantizeBlock(bij_buf, &q_scale.u[iEdge * nVar], &q_blocks.u[iEdge * blkSz]); const auto k_l = edge_ptr_l[iEdge]; - QuantizeBlock(bji_buf, &q_scale_l[k_l * nVar], &q_blocks_l[k_l * blkSz]); + QuantizeBlock(bji_buf, &q_scale.l[k_l * nVar], &q_blocks.l[k_l * blkSz]); return; } @@ -895,9 +889,9 @@ class CSysMatrix { bii[i] -= blk_i[k][i]; bjj[i] -= blk_j[k][i]; } - QuantizeBlock(blk_j[k], &q_scale_u[iEdge[k] * nVar], &q_blocks_u[iEdge[k] * blkSz]); + QuantizeBlock(blk_j[k], &q_scale.u[iEdge[k] * nVar], &q_blocks.u[iEdge[k] * blkSz]); const auto k_l = edge_ptr_l[iEdge[k]]; - QuantizeBlock(blk_i[k], &q_scale_l[k_l * nVar], &q_blocks_l[k_l * blkSz]); + QuantizeBlock(blk_i[k], &q_scale.l[k_l * nVar], &q_blocks.l[k_l * blkSz]); } else { auto bij = &mat.u[iEdge[k] * blkSz]; auto bji = &mat.l[edge_ptr_l[iEdge[k]] * blkSz]; @@ -938,9 +932,9 @@ class CSysMatrix { bij_buf[offset] = PassiveAssign(block_j[iVar][jVar] * scale); bji_buf[offset] = -PassiveAssign(block_i[iVar][jVar] * scale); } - QuantizeBlock(bij_buf, &q_scale_u[iEdge * nVar], &q_blocks_u[iEdge * blkSz]); + QuantizeBlock(bij_buf, &q_scale.u[iEdge * nVar], &q_blocks.u[iEdge * blkSz]); const auto k_l = edge_ptr_l[iEdge]; - QuantizeBlock(bji_buf, &q_scale_l[k_l * nVar], &q_blocks_l[k_l * blkSz]); + QuantizeBlock(bji_buf, &q_scale.l[k_l * nVar], &q_blocks.l[k_l * blkSz]); return; } @@ -1001,9 +995,9 @@ class CSysMatrix { if (mask[k] == 0) continue; if (quantized_mode) { - QuantizeBlock(blk_j[k], &q_scale_u[iEdge[k] * nVar], &q_blocks_u[iEdge[k] * blkSz]); + QuantizeBlock(blk_j[k], &q_scale.u[iEdge[k] * nVar], &q_blocks.u[iEdge[k] * blkSz]); const auto k_l = edge_ptr_l[iEdge[k]]; - QuantizeBlock(blk_i[k], &q_scale_l[k_l * nVar], &q_blocks_l[k_l * blkSz]); + QuantizeBlock(blk_i[k], &q_scale.l[k_l * nVar], &q_blocks.l[k_l * blkSz]); } else { ScalarType* bij = &mat.u[iEdge[k] * blkSz]; ScalarType* bji = &mat.l[edge_ptr_l[iEdge[k]] * blkSz]; diff --git a/Common/include/linear_algebra/CSysMatrix.inl b/Common/include/linear_algebra/CSysMatrix.inl index b5411016af6..c0ac196c961 100644 --- a/Common/include/linear_algebra/CSysMatrix.inl +++ b/Common/include/linear_algebra/CSysMatrix.inl @@ -147,8 +147,8 @@ FORCEINLINE void CSysMatrix::GaussElimination(unsigned long block_i, template FORCEINLINE void CSysMatrix::QuantizedGaussElimination(unsigned long block_i, ScalarType* rhs) const { ScalarType block[MAXNVAR * MAXNVAR]; - const QuantType* __restrict qs = &q_scale_d[block_i * nVar]; - const QuantType* __restrict qv = &q_blocks_d[block_i * nVar * nVar]; + const QuantType* __restrict qs = &q_scale.d[block_i * nVar]; + const QuantType* __restrict qv = &q_blocks.d[block_i * nVar * nVar]; for (auto r = 0ul; r < nVar; ++r) { const float row_scale = DecodeQuantScale(qs[r]); for (auto c = 0ul; c < nVar; ++c) block[r * nVar + c] = static_cast(qv[r * nVar + c] * row_scale); @@ -242,12 +242,12 @@ FORCEINLINE void CSysMatrix::QuantizedRowProduct(const CSysVector::QuantizedUpperProduct(const CSysVector< for (auto index = mat.row_ptr_u[row_i]; index < mat.row_ptr_u[row_i + 1]; index++) { auto col_j = mat.col_ind_u[index]; if (col_j < col_ub || col_j >= nPointDomain) { - QuantizedMatVecAdd(&q_scale_u[index * nVar], &q_blocks_u[index * nVar * nEqn], &vec[col_j * nEqn], prod); + QuantizedMatVecAdd(&q_scale.u[index * nVar], &q_blocks.u[index * nVar * nEqn], &vec[col_j * nEqn], prod); } } } @@ -272,7 +272,7 @@ FORCEINLINE void CSysMatrix::QuantizedLowerProduct(const CSysVector< for (auto index = mat.row_ptr_l[row_i]; index < mat.row_ptr_l[row_i + 1]; index++) { auto col_j = mat.col_ind_l[index]; if (col_j >= col_lb) { - QuantizedMatVecAdd(&q_scale_l[index * nVar], &q_blocks_l[index * nVar * nEqn], &vec[col_j * nEqn], prod); + QuantizedMatVecAdd(&q_scale.l[index * nVar], &q_blocks.l[index * nVar * nEqn], &vec[col_j * nEqn], prod); } } } @@ -281,5 +281,5 @@ template FORCEINLINE void CSysMatrix::QuantizedDiagonalProduct(const CSysVector& vec, unsigned long row_i, ScalarType* prod) const { for (auto iVar = 0ul; iVar < nVar; iVar++) prod[iVar] = 0.0; - QuantizedMatVecAdd(&q_scale_d[row_i * nVar], &q_blocks_d[row_i * nVar * nEqn], &vec[row_i * nEqn], prod); + QuantizedMatVecAdd(&q_scale.d[row_i * nVar], &q_blocks.d[row_i * nVar * nEqn], &vec[row_i * nEqn], prod); } diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 6de7de4178a..b8e55afd9e4 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -99,12 +99,8 @@ CSysMatrix::CSysMatrix() : rank(SU2_MPI::GetRank()), size(SU2_MPI::G ilu.d = nullptr; ilu.u = nullptr; - q_scale_l = nullptr; - q_blocks_l = nullptr; - q_scale_u = nullptr; - q_blocks_u = nullptr; - q_scale_d = nullptr; - q_blocks_d = nullptr; + q_scale = {}; + q_blocks = {}; invM = nullptr; d_invM = nullptr; @@ -122,23 +118,19 @@ CSysMatrix::~CSysMatrix() { SU2_ZONE_SCOPED delete[] omp_partitions; - auto freeHostLDU = [](LDU& m) { + auto freeHostLDU = [](auto& m) { MemoryAllocation::aligned_free(m.d); MemoryAllocation::aligned_free(m.l); MemoryAllocation::aligned_free(m.u); }; freeHostLDU(mat); freeHostLDU(ilu); + freeHostLDU(q_scale); + freeHostLDU(q_blocks); MemoryAllocation::aligned_free(invM); - MemoryAllocation::aligned_free(q_scale_l); - MemoryAllocation::aligned_free(q_blocks_l); - MemoryAllocation::aligned_free(q_scale_u); - MemoryAllocation::aligned_free(q_blocks_u); - MemoryAllocation::aligned_free(q_scale_d); - MemoryAllocation::aligned_free(q_blocks_d); if (useCuda) { - auto freeLDU = [](LDU& m) { + auto freeLDU = [](auto& m) { GPUMemoryAllocation::gpu_free(m.d); GPUMemoryAllocation::gpu_free(m.l); GPUMemoryAllocation::gpu_free(m.u); @@ -149,15 +141,11 @@ CSysMatrix::~CSysMatrix() { }; freeLDU(gpu); freeLDU(gpu_ilu); + freeLDU(d_q_scale); + 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_q_scale_l); - GPUMemoryAllocation::gpu_free(d_q_blocks_l); - GPUMemoryAllocation::gpu_free(d_q_scale_u); - GPUMemoryAllocation::gpu_free(d_q_blocks_u); - GPUMemoryAllocation::gpu_free(d_q_scale_d); - GPUMemoryAllocation::gpu_free(d_q_blocks_d); #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); @@ -262,12 +250,12 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi auto allocQ = [](QuantType*& ptr, unsigned long n) { ptr = MemoryAllocation::aligned_alloc(64, n * sizeof(QuantType)); }; - allocQ(q_scale_l, mat.nnz_l * nVar); - allocQ(q_blocks_l, mat.nnz_l * nVar * nEqn); - allocQ(q_scale_u, mat.nnz_u * nVar); - allocQ(q_blocks_u, mat.nnz_u * nVar * nEqn); - allocQ(q_scale_d, nPoint * nVar); - allocQ(q_blocks_d, nPoint * nVar * nEqn); + allocQ(q_scale.l, mat.nnz_l * nVar); + allocQ(q_blocks.l, mat.nnz_l * nVar * nEqn); + allocQ(q_scale.u, mat.nnz_u * nVar); + allocQ(q_blocks.u, mat.nnz_u * nVar * nEqn); + allocQ(q_scale.d, nPoint * nVar); + allocQ(q_blocks.d, nPoint * nVar * nEqn); } else { allocAndInit(mat.l, mat.nnz_l * nVar * nEqn); allocAndInit(mat.u, mat.nnz_u * nVar * nEqn); @@ -291,19 +279,18 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi GPUAllocAndCopy(gpu.col_ind_u, mat.col_ind_u, mat.nnz_u); if (quantized_mode) { - /*--- Device mirrors of the host quantized off-diagonal storage; gpu.l/gpu.u are not - * allocated (nothing would ever read them). The diagonal is quantized straight from - * gpu.d by QuantizeDiagonalBlocksGPU(), so there is no d_q_scale_d/d_q_blocks_d transfer - * to arrange, only the allocation here. ---*/ + /*--- Device mirrors of the host quantized storage; gpu.l/gpu.u are not allocated (nothing + * would ever read them). d_q_scale.d/d_q_blocks.d are uploaded from the host result once + * QuantizeDiagonalBlocks() has computed it, see the comment on those members. ---*/ auto GPUAllocQ = [](QuantType*& ptr, unsigned long n) { ptr = GPUMemoryAllocation::gpu_alloc(n * sizeof(QuantType)); }; - GPUAllocQ(d_q_scale_l, mat.nnz_l * nVar); - GPUAllocQ(d_q_blocks_l, mat.nnz_l * nVar * nEqn); - GPUAllocQ(d_q_scale_u, mat.nnz_u * nVar); - GPUAllocQ(d_q_blocks_u, mat.nnz_u * nVar * nEqn); - GPUAllocQ(d_q_scale_d, nPoint * nVar); - GPUAllocQ(d_q_blocks_d, nPoint * nVar * nEqn); + GPUAllocQ(d_q_scale.l, mat.nnz_l * nVar); + GPUAllocQ(d_q_blocks.l, mat.nnz_l * nVar * nEqn); + GPUAllocQ(d_q_scale.u, mat.nnz_u * nVar); + GPUAllocQ(d_q_blocks.u, mat.nnz_u * nVar * nEqn); + GPUAllocQ(d_q_scale.d, nPoint * nVar); + GPUAllocQ(d_q_blocks.d, nPoint * nVar * nEqn); } else { GPUAllocAndInit(gpu.l, mat.nnz_l * nVar * nEqn); GPUAllocAndInit(gpu.u, mat.nnz_u * nVar * nEqn); @@ -761,13 +748,20 @@ void CSysMatrix::QuantizeDiagonalBlocks() { if (!quantized_mode) return; + /*--- Q_LU_SGS / Q_JACOBI / Q_IDENTITY: L/U were quantized during assembly; only the diagonal needs quantization + * now. Always done on the host, even under CUDA: HtDTransfer() already kicked off the (larger) + * L/U quantized transfer asynchronously before Build() reached this point, so quantizing the + * diagonal here on the CPU overlaps with that transfer instead of waiting on it first. ---*/ + SU2_OMP_FOR_DYN(omp_heavy_size) + for (auto i = 0ul; i < nPointDomain; ++i) + QuantizeBlock(&mat.d[i * nVar * nVar], &q_scale.d[i * nVar], &q_blocks.d[i * nVar * nVar]); + END_SU2_OMP_FOR + if (useCuda) { #ifdef SU2_ENABLE_CUDA_KERNELS if constexpr (su2_gpu_capable_v) { - /*--- Quantize straight from gpu.d, see the d_q_scale_d/d_q_blocks_d comment in - * CSysMatrix.hpp for why this does not go through the host q_scale_d/q_blocks_d. ---*/ + /*--- Upload the host result just computed above. ---*/ SU2_DEVICE_REGION(QuantizeDiagonalBlocksGPU();) - return; } else { GPUNotAvailable(CURRENT_FUNCTION); } @@ -775,13 +769,6 @@ void CSysMatrix::QuantizeDiagonalBlocks() { GPUNotAvailable(CURRENT_FUNCTION); #endif } - - /*--- Q_LU_SGS / Q_JACOBI / Q_IDENTITY: L/U were quantized during assembly; only the diagonal needs quantization - * now. ---*/ - SU2_OMP_FOR_DYN(omp_heavy_size) - for (auto i = 0ul; i < nPointDomain; ++i) - QuantizeBlock(&mat.d[i * nVar * nVar], &q_scale_d[i * nVar], &q_blocks_d[i * nVar * nVar]); - END_SU2_OMP_FOR } template @@ -802,10 +789,10 @@ void CSysMatrix::SetValZero() { zeroChunk(mat.l, mat.nnz_l * nVar * nEqn); zeroChunk(mat.u, mat.nnz_u * nVar * nEqn); } else { - zeroChunk(q_scale_l, mat.nnz_l * nVar); - zeroChunk(q_scale_u, mat.nnz_l * nVar); - zeroChunk(q_blocks_l, mat.nnz_l * nVar * nEqn); - zeroChunk(q_blocks_u, mat.nnz_u * nVar * nEqn); + zeroChunk(q_scale.l, mat.nnz_l * nVar); + zeroChunk(q_scale.u, mat.nnz_l * nVar); + zeroChunk(q_blocks.l, mat.nnz_l * nVar * nEqn); + zeroChunk(q_blocks.u, mat.nnz_u * nVar * nEqn); } SU2_OMP_BARRIER } @@ -907,10 +894,10 @@ void CSysMatrix::DeleteValsRowi(unsigned long block_i, unsigned long if (quantized_mode) { for (auto k = mat.row_ptr_l[block_i]; k < mat.row_ptr_l[block_i + 1]; ++k) { - for (auto iVar = 0u; iVar < nEqn; iVar++) q_blocks_l[k * blkSz + row * nEqn + iVar] = 0; + for (auto iVar = 0u; iVar < nEqn; iVar++) q_blocks.l[k * blkSz + row * nEqn + iVar] = 0; } for (auto k = mat.row_ptr_u[block_i]; k < mat.row_ptr_u[block_i + 1]; ++k) { - for (auto iVar = 0u; iVar < nEqn; iVar++) q_blocks_u[k * blkSz + row * nEqn + iVar] = 0; + for (auto iVar = 0u; iVar < nEqn; iVar++) q_blocks.u[k * blkSz + row * nEqn + iVar] = 0; } } else { for (auto k = mat.row_ptr_l[block_i]; k < mat.row_ptr_l[block_i + 1]; ++k) { @@ -1605,9 +1592,9 @@ void CSysMatrix::SetDiagonalAsColumnSum() { for (auto j = 0ul; j < nEqn; ++j) d_i[i * nEqn + j] -= view(i, j); }; for (auto k_l = mat.row_ptr_l[iPoint]; k_l < mat.row_ptr_l[iPoint + 1]; ++k_l) - subtractTransp(l_to_u_transp[k_l], q_scale_u, q_blocks_u); + subtractTransp(l_to_u_transp[k_l], q_scale.u, q_blocks.u); for (auto k_u = mat.row_ptr_u[iPoint]; k_u < mat.row_ptr_u[iPoint + 1]; ++k_u) - subtractTransp(u_to_l_transp[k_u], q_scale_l, q_blocks_l); + subtractTransp(u_to_l_transp[k_u], q_scale.l, q_blocks.l); } } END_SU2_OMP_FOR diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index af0bc3ed91b..286aa5999f7 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -362,46 +362,6 @@ __device__ FORCEINLINE float DecodeQuantScaleDevice(int8_t e) { return __uint_as_float(bits); } -/*! - * \brief Device counterpart of EncodeQuantBlock (CSysMatrix.hpp), bit-identical: quantizes one - * nVar x nVar block (row-major, pointed to by \p blk) into per-row int8 storage. One - * thread handles one whole block (i.e. one matrix row). - */ -template -__device__ void EncodeQuantBlockDevice(const ScalarType* __restrict__ blk, int8_t* __restrict__ qs, - int8_t* __restrict__ qv, unsigned long nVar) { - for (auto r = 0ul; r < nVar; ++r) { - constexpr uint32_t eps_bits = 0x34000000u; - uint32_t max_abs_bits = eps_bits; - for (auto c = 0ul; c < nVar; ++c) { - const float fv = static_cast(blk[r * nVar + c]); - max_abs_bits = max(max_abs_bits, __float_as_uint(fv) & 0x7FFFFFFFu); - } - const int e = min(127, max(-128, static_cast(max_abs_bits >> 23) - 133)); - qs[r] = static_cast(e); - const float inv_rscale = __uint_as_float(static_cast(127 - e) << 23); - for (auto c = 0ul; c < nVar; ++c) { - qv[r * nVar + c] = - static_cast(max(-128.f, min(127.f, roundf(static_cast(blk[r * nVar + c]) * inv_rscale)))); - } - } -} - -/*! - * \brief Quantize the diagonal blocks straight from the device diagonal (gpu.d), device version - * of CSysMatrix::QuantizeBlock applied row-by-row (host QuantizeDiagonalBlocks). - * One thread per row. - */ -template -__global__ void QuantizeDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, - const ScalarType* __restrict__ mat_d, int8_t* __restrict__ q_scale_d, - int8_t* __restrict__ q_blocks_d) { - const auto iRow = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; - if (iRow >= nRows) return; - - EncodeQuantBlockDevice(mat_d + iRow * nVar * nVar, q_scale_d + iRow * nVar, q_blocks_d + iRow * nVar * nVar, nVar); -} - /*! * \brief Quantized block-LDU SpMV kernel: y[iRow] = (L + D + U) * x per block-row, reading int8 * row-scaled quantized blocks instead of full precision ones. Device version of @@ -493,15 +453,12 @@ void CSysMatrix::QuantizeDiagonalBlocksGPU() { 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. ---*/ - constexpr unsigned threadsPerBlock = 128; - const auto blocks = static_cast((nPointDomain + threadsPerBlock - 1) / threadsPerBlock); - QuantizeDiagonalBlocksKernel - <<>>(nPointDomain, nVar, gpu.d, d_q_scale_d, d_q_blocks_d); - /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ - gpuErrChk(cudaStreamSynchronize(nullptr)); - gpuErrChk(cudaGetLastError()); + /*--- The diagonal is quantized on the host, in QuantizeDiagonalBlocks(), right before this is + * called; this just uploads the result. Doing it on the host lets it run while the (larger) + * async L/U transfer kicked off earlier by HtDTransfer() is still in flight on the device. ---*/ + gpuErrChk(cudaMemcpy(d_q_scale.d, q_scale.d, sizeof(QuantType) * nPointDomain * nVar, cudaMemcpyHostToDevice)); + gpuErrChk( + cudaMemcpy(d_q_blocks.d, q_blocks.d, sizeof(QuantType) * nPointDomain * nVar * nVar, cudaMemcpyHostToDevice)); } template @@ -644,15 +601,19 @@ void CSysMatrix::HtDTransfer(bool trigger) const { gpuErrChk(cudaMemcpy(gpu.d, mat.d, sizeof(ScalarType) * nPoint * nVar * nEqn, cudaMemcpyHostToDevice)); if (quantized_mode) { /*--- No gpu.l/gpu.u to transfer (never allocated); mirror the host quantized off-diagonal - * storage instead. The diagonal mirrors (d_q_scale_d/d_q_blocks_d) are NOT transferred here, - * they are populated straight from gpu.d by QuantizeDiagonalBlocksGPU(), see the comment on - * those members in CSysMatrix.hpp for why. ---*/ - gpuErrChk(cudaMemcpy(d_q_scale_l, q_scale_l, sizeof(QuantType) * mat.nnz_l * nVar, cudaMemcpyHostToDevice)); - gpuErrChk( - cudaMemcpy(d_q_blocks_l, q_blocks_l, sizeof(QuantType) * mat.nnz_l * nVar * nEqn, cudaMemcpyHostToDevice)); - gpuErrChk(cudaMemcpy(d_q_scale_u, q_scale_u, sizeof(QuantType) * mat.nnz_u * nVar, cudaMemcpyHostToDevice)); - gpuErrChk( - cudaMemcpy(d_q_blocks_u, q_blocks_u, sizeof(QuantType) * mat.nnz_u * nVar * nEqn, cudaMemcpyHostToDevice)); + * storage instead. Issued as async copies on the default stream, then left in flight: the + * caller (CSysMatrixVectorProduct's constructor) returns right after this, and QuantizeDiag- + * onalBlocks() -> Build() quantizes the diagonal on the host next, before anything is + * launched on the device again. Any later kernel that reads d_q_scale/d_q_blocks (also issued + * on the default stream) still waits for these correctly, by stream ordering, without an + * explicit sync here; the diagonal mirrors (d_q_scale.d/d_q_blocks.d) are uploaded once that + * host quantization is done, by QuantizeDiagonalBlocksGPU(). ---*/ + gpuErrChk(cudaMemcpyAsync(d_q_scale.l, q_scale.l, sizeof(QuantType) * mat.nnz_l * nVar, cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpyAsync(d_q_blocks.l, q_blocks.l, sizeof(QuantType) * mat.nnz_l * nVar * nEqn, + cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpyAsync(d_q_scale.u, q_scale.u, sizeof(QuantType) * mat.nnz_u * nVar, cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpyAsync(d_q_blocks.u, q_blocks.u, sizeof(QuantType) * mat.nnz_u * nVar * nEqn, + cudaMemcpyHostToDevice)); } else { gpuErrChk(cudaMemcpy(gpu.l, mat.l, sizeof(ScalarType) * mat.nnz_l * nVar * nEqn, cudaMemcpyHostToDevice)); gpuErrChk(cudaMemcpy(gpu.u, mat.u, sizeof(ScalarType) * mat.nnz_u * nVar * nEqn, cudaMemcpyHostToDevice)); @@ -671,8 +632,8 @@ void CSysMatrix::MatrixVectorProductGPU(const CSysVector dim3 gridDim(static_cast(nPointDomain), 1, 1); if (quantized_mode) { QuantizedBlockLDU_SpMV_kernel<<>>( - nPointDomain, nVar, gpu.row_ptr_l, gpu.col_ind_l, d_q_scale_l, d_q_blocks_l, d_q_scale_d, d_q_blocks_d, - gpu.row_ptr_u, gpu.col_ind_u, d_q_scale_u, d_q_blocks_u, d_vec, d_prod); + nPointDomain, nVar, gpu.row_ptr_l, gpu.col_ind_l, d_q_scale.l, d_q_blocks.l, d_q_scale.d, d_q_blocks.d, + gpu.row_ptr_u, gpu.col_ind_u, d_q_scale.u, d_q_blocks.u, d_vec, d_prod); } else { BlockLDU_SpMV_kernel<<>>( nPointDomain, nVar, gpu.row_ptr_l, gpu.col_ind_l, gpu.l, gpu.d, From 9bdd89c200916e17d4053b746e643f77a184429f Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 11 Aug 2026 21:51:45 -0700 Subject: [PATCH 05/30] Fold quantized diagonal build into the base preconditioner wrappers CQuantizedLUSGSPreconditioner, CQuantizedJacobiPreconditioner and CQuantizedIdentityPreconditioner each differed from their plain counterpart only by an extra QuantizeDiagonalBlocks() call in Build(). Since that call is already a no-op unless the matrix was actually set up for quantization, push it into the plain wrappers instead: - CLU_SGSPreconditioner/CIdentityPreconditioner now call QuantizeDiagonalBlocks() directly in Build(). - CSysMatrix::BuildJacobiPreconditioner() calls QuantizeDiagonalBlocks() itself (before the jacobi_on_device branch, since diagonal quantization is independent of invM and jacobi_on_device returns early), so CJacobiPreconditioner needs no change at all. This removes three duplicate classes; the factory now maps Q_LU_SGS/ Q_JACOBI/Q_IDENTITY to the same wrapper as their plain counterpart, and CSysSolve's discrete-adjoint transpose-rebuild switch merges Q_JACOBI into the JACOBI/LINELET case. Re-validated on the NACA0012 case: LU_SGS/Q_LU_SGS (host) and CPU/GPU x plain/quantized Jacobi and GPU quantized-identity all reproduce the same CL/CD as before this refactor. Co-Authored-By: Claude Sonnet 5 --- .../linear_algebra/CPreconditioner.hpp | 111 +++--------------- Common/src/linear_algebra/CSysMatrix.cpp | 5 + Common/src/linear_algebra/CSysSolve.cpp | 9 +- 3 files changed, 24 insertions(+), 101 deletions(-) diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index e1c9a27718d..98e252749ff 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -107,36 +107,26 @@ CPreconditioner::~CPreconditioner() {} /*! * \class CIdentityPreconditioner * \brief No-op preconditioner used when Krylov solvers run without preconditioning. + * \note Also serves Q_IDENTITY: Build() requests quantization of the diagonal blocks, needed by + * the matrix-vector product shared with the Krylov solver (off diagonals are quantized on the + * fly during assembly), even though this preconditioner's own operation is a no-op either way; + * CSysMatrix::QuantizeDiagonalBlocks() itself no-ops unless the matrix was actually set up for + * quantization (Q_IDENTITY/Q_JACOBI/Q_LU_SGS), so this is free for plain IDENTITY. */ template class CIdentityPreconditioner final : public CPreconditioner { - public: - inline void operator()(const CSysVector& u, CSysVector& v) const override { v = u; } - - inline bool IsIdentity() const override { return true; } -}; - -/*! - * \class CQuantizedIdentityPreconditioner - * \brief No-op preconditioner, same as CIdentityPreconditioner, but requests quantized (int8) - * off-diagonal storage for the matrix-vector product shared with the Krylov solver, as Q_LU_SGS/Q_JACOBI. - */ -template -class CQuantizedIdentityPreconditioner final : public CPreconditioner { private: CSysMatrix& sparse_matrix; public: - inline explicit CQuantizedIdentityPreconditioner(CSysMatrix& matrix_ref) : sparse_matrix(matrix_ref) {} + inline explicit CIdentityPreconditioner(CSysMatrix& matrix_ref) : sparse_matrix(matrix_ref) {} - CQuantizedIdentityPreconditioner() = delete; + CIdentityPreconditioner() = delete; inline void operator()(const CSysVector& u, CSysVector& v) const override { v = u; } inline bool IsIdentity() const override { return true; } - /*! \brief Quantize the diagonal blocks (off diagonals are quantized on the fly), needed by the - * quantized matrix-vector product even though this preconditioner itself is a no-op. */ inline void Build() override { sparse_matrix.QuantizeDiagonalBlocks(); } }; @@ -181,49 +171,13 @@ class CJacobiPreconditioner final : public CPreconditioner { } /*! - * \note Request the associated matrix to build the preconditioner. + * \note Request the associated matrix to build the preconditioner. Also serves Q_JACOBI: + * BuildJacobiPreconditioner() quantizes the diagonal blocks itself when the matrix was + * set up for it, so there is nothing extra to do here for the quantized case. */ inline void Build() override { sparse_matrix.BuildJacobiPreconditioner(); } }; -/*! - * \class CQuantizedJacobiPreconditioner - * \brief Specialization of preconditioner that uses CSysMatrix class. - * \note The preconditioner operation itself is identical to CJacobiPreconditioner (it only ever - * applies the full precision inverse diagonal); quantization here only affects the off-diagonal - * blocks used by the matrix-vector product shared with the Krylov solver, exactly as Q_LU_SGS. - */ -template -class CQuantizedJacobiPreconditioner final : public CPreconditioner { - private: - CSysMatrix& sparse_matrix; - CGeometry* geometry; - const CConfig* config; - - public: - inline CQuantizedJacobiPreconditioner(CSysMatrix& matrix_ref, CGeometry* geometry_ref, - const CConfig* config_ref) - : sparse_matrix(matrix_ref) { - if ((geometry_ref == nullptr) || (config_ref == nullptr)) - SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION); - geometry = geometry_ref; - config = config_ref; - } - - CQuantizedJacobiPreconditioner() = delete; - - inline void operator()(const CSysVector& u, CSysVector& v) const override { - sparse_matrix.ComputeJacobiPreconditioner(u, v, geometry, config); - } - - /*! \brief Build the (full precision) inverse diagonal, then quantize the diagonal blocks so the - * quantized matrix-vector product (off diagonals are quantized on the fly) is consistent. */ - inline void Build() override { - sparse_matrix.BuildJacobiPreconditioner(); - sparse_matrix.QuantizeDiagonalBlocks(); - } -}; - /*! * \class CILUPreconditioner * \brief Specialization of preconditioner that uses CSysMatrix class @@ -310,36 +264,11 @@ class CLU_SGSPreconditioner final : public CPreconditioner { inline void operator()(const CSysVector& u, CSysVector& v) const override { ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputeLU_SGSPreconditioner(u, v, geometry, config); }); } -}; - -/*! - * \class CQuantizedLUSGSPreconditioner - * \brief Specialization of preconditioner that uses CSysMatrix class. - */ -template -class CQuantizedLUSGSPreconditioner final : public CPreconditioner { - private: - CSysMatrix& sparse_matrix; - CGeometry* geometry; - const CConfig* config; - public: - inline CQuantizedLUSGSPreconditioner(CSysMatrix& matrix_ref, CGeometry* geometry_ref, - const CConfig* config_ref) - : sparse_matrix(matrix_ref) { - if ((geometry_ref == nullptr) || (config_ref == nullptr)) - SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION); - geometry = geometry_ref; - config = config_ref; - } - - CQuantizedLUSGSPreconditioner() = delete; - - inline void operator()(const CSysVector& u, CSysVector& v) const override { - ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputeLU_SGSPreconditioner(u, v, geometry, config); }); - } - - /*! \brief Quantize the diagonal blocks (off diagonals are quantized on the fly). */ + /*! + * \note Also serves Q_LU_SGS: quantizes the diagonal blocks (off diagonals are quantized on + * the fly during assembly); a no-op for plain LU_SGS. + */ inline void Build() override { sparse_matrix.QuantizeDiagonalBlocks(); } }; @@ -466,25 +395,19 @@ CPreconditioner* CPreconditioner::Create(ENUM_LINEAR_SOL switch (kind) { case IDENTITY: - prec = new CIdentityPreconditioner(); - break; case Q_IDENTITY: - prec = new CQuantizedIdentityPreconditioner(jacobian); + prec = new CIdentityPreconditioner(jacobian); break; case JACOBI: - prec = new CJacobiPreconditioner(jacobian, geometry, config); - break; case Q_JACOBI: - prec = new CQuantizedJacobiPreconditioner(jacobian, geometry, config); + prec = new CJacobiPreconditioner(jacobian, geometry, config); break; case LINELET: prec = new CLineletPreconditioner(jacobian, geometry, config); break; case LU_SGS: - prec = new CLU_SGSPreconditioner(jacobian, geometry, config); - break; case Q_LU_SGS: - prec = new CQuantizedLUSGSPreconditioner(jacobian, geometry, config); + prec = new CLU_SGSPreconditioner(jacobian, geometry, config); break; case ILU: prec = new CILUPreconditioner(jacobian, geometry, config); diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index b8e55afd9e4..625f190bdd8 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -969,6 +969,11 @@ template void CSysMatrix::BuildJacobiPreconditioner() { SU2_ZONE_SCOPED + /*--- Independent of invM (reads/quantizes mat.d, a no-op unless quantized_mode); done first, + * unconditionally, so it runs whichever branch below builds invM (in particular the + * jacobi_on_device one, which returns early). ---*/ + QuantizeDiagonalBlocks(); + if (jacobi_on_device) { #ifdef SU2_ENABLE_CUDA_KERNELS if constexpr (su2_gpu_capable_v) { diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index a008a503169..3c9176d8716 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -1566,14 +1566,9 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con break; case JACOBI: case LINELET: - if (RequiresTranspose) Jacobian.BuildJacobiPreconditioner(); - break; case Q_JACOBI: - /*--- Rebuild the full-precision inverse diagonal and re-quantize it, mirroring Build(). ---*/ - if (RequiresTranspose) { - Jacobian.BuildJacobiPreconditioner(); - Jacobian.QuantizeDiagonalBlocks(); - } + /*--- BuildJacobiPreconditioner() quantizes the diagonal itself when needed. ---*/ + if (RequiresTranspose) Jacobian.BuildJacobiPreconditioner(); break; case LU_SGS: case Q_LU_SGS: From 76c27f22ba3b22a6c31ba2d1ec2415dd0ae17c32 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 11 Aug 2026 22:00:47 -0700 Subject: [PATCH 06/30] Inline the diagonal quantization upload instead of a separate GPU hook QuantizeDiagonalBlocksGPU() had shrunk to two cudaMemcpy calls with no kernel involved, so it doesn't need the usual CSysMatrixGPU.cu hook indirection: cudaMemcpy (unlike a kernel launch) is a plain CUDA runtime call, already available in CSysMatrix.cpp transitively via allocation_toolbox.hpp -> GPUComms.cuh, the same way Initialize() already calls cudaMalloc/cudaMemcpy directly for the other GPU buffers. Inlines the upload at the end of QuantizeDiagonalBlocks() and removes the now-empty hook (declaration, .cu implementation, and its INSTANTIATE_MATRIX entry). Re-validated on the NACA0012 case: GPU quantized-Jacobi and quantized- identity reproduce the same CL/CD as before this change. Co-Authored-By: Claude Sonnet 5 --- Common/include/linear_algebra/CSysMatrix.hpp | 11 ++--------- Common/src/linear_algebra/CSysMatrix.cpp | 10 ++++++++-- Common/src/linear_algebra/CSysMatrixGPU.cu | 17 +---------------- 3 files changed, 11 insertions(+), 27 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 0ac1d2e7bb9..f22e59a3fcb 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -282,8 +282,8 @@ class CSysMatrix { * useCuda (currently only reachable for Q_JACOBI/Q_IDENTITY, Q_LU_SGS stays host-only). * d_q_scale.l/.u and d_q_blocks.l/.u are plain device-side copies of q_scale.l/.u and * q_blocks.l/.u, uploaded *asynchronously* by HtDTransfer() so that transfer can overlap with - * the host quantizing the diagonal (see QuantizeDiagonalBlocks()); d_q_scale.d/d_q_blocks.d are - * that host result, uploaded once ready by QuantizeDiagonalBlocksGPU(). */ + * the host quantizing the diagonal; d_q_scale.d/d_q_blocks.d are that host result, uploaded (a + * plain cudaMemcpy, not a kernel) once ready, at the end of QuantizeDiagonalBlocks(). */ LDU d_q_scale; LDU d_q_blocks; @@ -569,13 +569,6 @@ class CSysMatrix { void MatrixVectorProductGPU(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const; - /*! - * \brief Upload the already host-quantized diagonal (q_scale.d/q_blocks.d, computed by - * QuantizeDiagonalBlocks() just before calling this) into d_q_scale.d/d_q_blocks.d. - * Only reachable when quantized_mode && useCuda. - */ - void QuantizeDiagonalBlocksGPU(); - /*! * \brief Build the Jacobi preconditioner on the device, from the device copy of the matrix. * \note Requires the device matrix to be up to date, see HtDTransfer. diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 625f190bdd8..15ac1f0c433 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -760,8 +760,14 @@ void CSysMatrix::QuantizeDiagonalBlocks() { if (useCuda) { #ifdef SU2_ENABLE_CUDA_KERNELS if constexpr (su2_gpu_capable_v) { - /*--- Upload the host result just computed above. ---*/ - SU2_DEVICE_REGION(QuantizeDiagonalBlocksGPU();) + /*--- Just an upload of the host result above, no computation, so a plain CUDA runtime call + * (available here via GPUComms.cuh, transitively included through allocation_toolbox.hpp) + * rather than a kernel dispatched through a CSysMatrixGPU.cu hook. ---*/ + BEGIN_SU2_DEVICE_REGION + gpuErrChk(cudaMemcpy(d_q_scale.d, q_scale.d, sizeof(QuantType) * nPointDomain * nVar, cudaMemcpyHostToDevice)); + gpuErrChk( + cudaMemcpy(d_q_blocks.d, q_blocks.d, sizeof(QuantType) * nPointDomain * nVar * nVar, cudaMemcpyHostToDevice)); + END_SU2_DEVICE_REGION } else { GPUNotAvailable(CURRENT_FUNCTION); } diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 286aa5999f7..513993dcba1 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -447,20 +447,6 @@ void CSysMatrix::BuildJacobiPreconditionerGPU() { gpuErrChk(cudaGetLastError()); } -template -void CSysMatrix::QuantizeDiagonalBlocksGPU() { - SU2_ZONE_SCOPED - - if (nPointDomain == 0) return; - - /*--- The diagonal is quantized on the host, in QuantizeDiagonalBlocks(), right before this is - * called; this just uploads the result. Doing it on the host lets it run while the (larger) - * async L/U transfer kicked off earlier by HtDTransfer() is still in flight on the device. ---*/ - gpuErrChk(cudaMemcpy(d_q_scale.d, q_scale.d, sizeof(QuantType) * nPointDomain * nVar, cudaMemcpyHostToDevice)); - gpuErrChk( - cudaMemcpy(d_q_blocks.d, q_blocks.d, sizeof(QuantType) * nPointDomain * nVar * nVar, cudaMemcpyHostToDevice)); -} - template void CSysMatrix::BuildILUPreconditionerGPU() { SU2_ZONE_SCOPED @@ -607,7 +593,7 @@ void CSysMatrix::HtDTransfer(bool trigger) const { * launched on the device again. Any later kernel that reads d_q_scale/d_q_blocks (also issued * on the default stream) still waits for these correctly, by stream ordering, without an * explicit sync here; the diagonal mirrors (d_q_scale.d/d_q_blocks.d) are uploaded once that - * host quantization is done, by QuantizeDiagonalBlocksGPU(). ---*/ + * host quantization is done, by a plain cudaMemcpy at the end of QuantizeDiagonalBlocks(). ---*/ gpuErrChk(cudaMemcpyAsync(d_q_scale.l, q_scale.l, sizeof(QuantType) * mat.nnz_l * nVar, cudaMemcpyHostToDevice)); gpuErrChk(cudaMemcpyAsync(d_q_blocks.l, q_blocks.l, sizeof(QuantType) * mat.nnz_l * nVar * nEqn, cudaMemcpyHostToDevice)); @@ -650,7 +636,6 @@ template void CSysMatrix::MatrixVectorProductGPU(const CSysVector& v CSysVector& prod, \ CGeometry* geometry, \ const CConfig* config) const; \ -template void CSysMatrix::QuantizeDiagonalBlocksGPU(); \ template void CSysMatrix::BuildJacobiPreconditionerGPU(); \ template void CSysMatrix::BuildILUPreconditionerGPU(); \ template void CSysMatrix::ComputeILUPreconditionerGPU(const CSysVector& vec, \ From c62a413f4a27a0dea64af04378617177a414cba7 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 14 Aug 2026 22:13:59 -0700 Subject: [PATCH 07/30] Pin the quantized L/U/D host arrays so the async H2D transfer actually overlaps cudaMemcpyAsync is only genuinely asynchronous (non-blocking on the host) when the host-side buffer is page-locked; from regular pageable memory (what q_scale/q_blocks were using) the driver silently falls back to a synchronous staged copy, so the async transfer added earlier was never actually overlapping with the host quantizing the diagonal. Adds GPUMemoryAllocation::pinned_alloc/pinned_free (cudaMallocHost/ cudaFreeHost) and switches q_scale/q_blocks to use them instead of aligned_alloc whenever useCuda, matched by conditional frees in the destructor. Also makes the diagonal upload async (now pinned too) and adds an explicit cudaStreamSynchronize after it, matching the pattern every other GPU-touching function in this file already follows, rather than relying solely on implicit default-stream ordering. Verified with nsys (gputrace) that the transfers now show up with SrcMemKd=Pinned. Also chased a small (~0.01%) CL/CD difference this introduced for the GPU quantized-identity case: confirmed via a CPU- only, non-CUDA run of the same case that varying OMP_NUM_THREADS alone produces a difference of the same order, so this is pre-existing floating-point non-associativity sensitivity in that (deliberately unpreconditioned) test case, not a bug in the new transfer code - plain/quantized Jacobi (which is actually preconditioned) stay bit-identical throughout. Co-Authored-By: Claude Sonnet 5 --- Common/include/linear_algebra/CSysMatrix.hpp | 10 +++-- .../include/toolboxes/allocation_toolbox.hpp | 35 +++++++++++++++ Common/src/linear_algebra/CSysMatrix.cpp | 44 ++++++++++++++++--- Common/src/linear_algebra/CSysMatrixGPU.cu | 13 ++++-- 4 files changed, 88 insertions(+), 14 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index f22e59a3fcb..95db350def8 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -273,7 +273,10 @@ class CSysMatrix { static constexpr bool quantized_mode = false; #endif /*!< \brief Per-row exponents; .l/.u sized [nnz_l/u * nVar], .d [nPoint * nVar] - * (populated by QuantizeDiagonalBlocks(), always on the host, see below). */ + * (populated by QuantizeDiagonalBlocks(), always on the host, see below). + * Pinned (cudaMallocHost) rather than aligned_alloc when useCuda, see Initialize(), + * so the async uploads below are genuinely asynchronous instead of silently + * blocking (cudaMemcpyAsync only overlaps with the host from pinned memory). */ LDU q_scale; /*!< \brief Quantized block entries; .l/.u sized [nnz_l/u * nVar * nEqn], .d [nPoint * nVar * nEqn]. */ LDU q_blocks; @@ -282,8 +285,9 @@ class CSysMatrix { * useCuda (currently only reachable for Q_JACOBI/Q_IDENTITY, Q_LU_SGS stays host-only). * d_q_scale.l/.u and d_q_blocks.l/.u are plain device-side copies of q_scale.l/.u and * q_blocks.l/.u, uploaded *asynchronously* by HtDTransfer() so that transfer can overlap with - * the host quantizing the diagonal; d_q_scale.d/d_q_blocks.d are that host result, uploaded (a - * plain cudaMemcpy, not a kernel) once ready, at the end of QuantizeDiagonalBlocks(). */ + * the host quantizing the diagonal; d_q_scale.d/d_q_blocks.d are that host result, likewise + * uploaded asynchronously (a plain cudaMemcpyAsync, not a kernel), at the end of + * QuantizeDiagonalBlocks(). */ LDU d_q_scale; LDU d_q_blocks; diff --git a/Common/include/toolboxes/allocation_toolbox.hpp b/Common/include/toolboxes/allocation_toolbox.hpp index 9c357405b75..66db0a04035 100644 --- a/Common/include/toolboxes/allocation_toolbox.hpp +++ b/Common/include/toolboxes/allocation_toolbox.hpp @@ -144,4 +144,39 @@ inline T* gpu_alloc_cpy(const T* src_ptr, size_t size) noexcept { return static_cast(ptr); } + +/*! + * \brief Page-locked ("pinned") host memory allocation. + * \note Unlike regular (pageable) host memory, cudaMemcpyAsync from/to a pinned buffer is + * actually asynchronous with respect to the host thread; from pageable memory the driver + * silently falls back to a synchronous staged copy. Only worth it for host buffers that + * are the source/destination of an async transfer meant to overlap with other host work. + * \param[in] size in bytes. + * \tparam ZeroInit, initialize memory to 0. + * \return Pointer to memory, always use pinned_free to deallocate. + */ +template +inline T* pinned_alloc(size_t size) noexcept { + void* ptr = nullptr; + +#if defined(HAVE_CUDA) + gpuErrChk(cudaMallocHost((void**)(&ptr), size)); + if (ZeroInit) memset(ptr, 0, size); +#else + return 0; +#endif + + return static_cast(ptr); +} + +/*! + * \brief Free memory allocated with pinned_alloc. + * \param[in] ptr, pointer to memory we want to release. + */ +template +inline void pinned_free(T* ptr) noexcept { +#ifdef HAVE_CUDA + gpuErrChk(cudaFreeHost((void*)ptr)); +#endif +} } // namespace GPUMemoryAllocation diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 15ac1f0c433..a4c52d1e134 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -125,10 +125,24 @@ CSysMatrix::~CSysMatrix() { }; freeHostLDU(mat); freeHostLDU(ilu); - freeHostLDU(q_scale); - freeHostLDU(q_blocks); MemoryAllocation::aligned_free(invM); + /*--- q_scale/q_blocks are pinned (cudaMallocHost) rather than aligned_alloc when useCuda, see + * the comment in Initialize(); free with the matching deallocator either way. ---*/ + auto freeQuantLDU = [this](auto& m) { + if (useCuda) { + GPUMemoryAllocation::pinned_free(m.d); + GPUMemoryAllocation::pinned_free(m.l); + GPUMemoryAllocation::pinned_free(m.u); + } else { + MemoryAllocation::aligned_free(m.d); + MemoryAllocation::aligned_free(m.l); + MemoryAllocation::aligned_free(m.u); + } + }; + freeQuantLDU(q_scale); + freeQuantLDU(q_blocks); + if (useCuda) { auto freeLDU = [](auto& m) { GPUMemoryAllocation::gpu_free(m.d); @@ -247,8 +261,16 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi #ifndef CODI_REVERSE_TYPE quantized_mode = true; #endif - auto allocQ = [](QuantType*& ptr, unsigned long n) { - ptr = MemoryAllocation::aligned_alloc(64, n * sizeof(QuantType)); + /*--- Pinned (page-locked) when useCuda: HtDTransfer()/QuantizeDiagonalBlocks() upload these + * with cudaMemcpyAsync, which is only genuinely asynchronous from pinned host memory (see + * GPUMemoryAllocation::pinned_alloc); from regular pageable memory it silently degrades to a + * blocking copy, defeating the overlap with host-side diagonal quantization. ---*/ + auto allocQ = [useCuda = this->useCuda](QuantType*& ptr, unsigned long n) { + if (useCuda) { + ptr = GPUMemoryAllocation::pinned_alloc(n * sizeof(QuantType)); + } else { + ptr = MemoryAllocation::aligned_alloc(64, n * sizeof(QuantType)); + } }; allocQ(q_scale.l, mat.nnz_l * nVar); allocQ(q_blocks.l, mat.nnz_l * nVar * nEqn); @@ -762,11 +784,19 @@ void CSysMatrix::QuantizeDiagonalBlocks() { if constexpr (su2_gpu_capable_v) { /*--- Just an upload of the host result above, no computation, so a plain CUDA runtime call * (available here via GPUComms.cuh, transitively included through allocation_toolbox.hpp) - * rather than a kernel dispatched through a CSysMatrixGPU.cu hook. ---*/ + * rather than a kernel dispatched through a CSysMatrixGPU.cu hook. Async (q_scale.d/ + * q_blocks.d are pinned, see Initialize()) so issuing it does not block the host; the sync + * right after does not undo that overlap; by this point the host has already spent the + * whole diagonal-quantization loop above letting the earlier (larger) L/U transfer, kicked + * off by HtDTransfer(), drain in the background; it only makes sure that transfer and this + * one are actually finished before Build() returns, matching every other GPU-touching + * function in this file (they all sync at the end, see the comment on those calls). ---*/ BEGIN_SU2_DEVICE_REGION - gpuErrChk(cudaMemcpy(d_q_scale.d, q_scale.d, sizeof(QuantType) * nPointDomain * nVar, cudaMemcpyHostToDevice)); gpuErrChk( - cudaMemcpy(d_q_blocks.d, q_blocks.d, sizeof(QuantType) * nPointDomain * nVar * nVar, cudaMemcpyHostToDevice)); + cudaMemcpyAsync(d_q_scale.d, q_scale.d, sizeof(QuantType) * nPointDomain * nVar, cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpyAsync(d_q_blocks.d, q_blocks.d, sizeof(QuantType) * nPointDomain * nVar * nVar, + cudaMemcpyHostToDevice)); + gpuErrChk(cudaStreamSynchronize(nullptr)); END_SU2_DEVICE_REGION } else { GPUNotAvailable(CURRENT_FUNCTION); diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 513993dcba1..9c2d0fb0639 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -590,10 +590,15 @@ void CSysMatrix::HtDTransfer(bool trigger) const { * storage instead. Issued as async copies on the default stream, then left in flight: the * caller (CSysMatrixVectorProduct's constructor) returns right after this, and QuantizeDiag- * onalBlocks() -> Build() quantizes the diagonal on the host next, before anything is - * launched on the device again. Any later kernel that reads d_q_scale/d_q_blocks (also issued - * on the default stream) still waits for these correctly, by stream ordering, without an - * explicit sync here; the diagonal mirrors (d_q_scale.d/d_q_blocks.d) are uploaded once that - * host quantization is done, by a plain cudaMemcpy at the end of QuantizeDiagonalBlocks(). ---*/ + * launched on the device again. This is only genuinely asynchronous (i.e. the host thread + * does not block here waiting for the copy) because q_scale.l/q_blocks.l/q_scale.u/ + * q_blocks.u are pinned host memory, see the comment on those members / Initialize() - + * cudaMemcpyAsync silently degrades to a blocking copy from regular pageable memory. Any + * later kernel that reads d_q_scale/d_q_blocks is issued on the same default stream too, so + * strictly it would not need to wait for these explicitly; QuantizeDiagonalBlocks() still + * syncs right after uploading the diagonal mirrors (d_q_scale.d/d_q_blocks.d, the same way, + * async from pinned memory) to be safe and consistent with every other GPU-touching function + * in this file, all of which sync before returning. ---*/ gpuErrChk(cudaMemcpyAsync(d_q_scale.l, q_scale.l, sizeof(QuantType) * mat.nnz_l * nVar, cudaMemcpyHostToDevice)); gpuErrChk(cudaMemcpyAsync(d_q_blocks.l, q_blocks.l, sizeof(QuantType) * mat.nnz_l * nVar * nEqn, cudaMemcpyHostToDevice)); From a1aaabf322dbc3a5311249e2113f384a56eb19a5 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 14 Aug 2026 22:43:27 -0700 Subject: [PATCH 08/30] Batch rows per block and vectorize int8 loads in the SpMV kernels Both BlockLDU_SpMV_kernel and QuantizedBlockLDU_SpMV_kernel launched one CUDA block per matrix row with blockDim.x == nVar (~4-6 threads), leaving most of a warp idle and capping occupancy at a few resident, mostly-empty warps per SM - latency-bound well before DRAM bandwidth becomes the limit, which is why the quantized kernel (4x fewer matrix bytes) wasn't measurably faster than the float one: shrinking the matrix doesn't help a kernel that isn't bandwidth-bound in the first place, and int8 element loads cost the same one-load-per-element as the float32 ones did. Rows are now batched into blocks of ~128 threads (threadIdx.x maps to (row-within-block, output variable) via divmod by nVar), the same batching for both kernels since it's the float kernel's launch shape that needs fixing first. The quantized kernel additionally reads each row of quantized mantissas 4 bytes at a time (memcpy into a uint32_t, row bytes aren't generally 4-byte aligned) instead of one byte at a time; true dp4a doesn't apply since only the matrix is quantized, not the vector, so accumulation stays one scalar FMA per element in the same order as before (bit-identical results) - the win is fewer, wider load instructions. Verified bit-identical results against the pre-change build on the NACA0012 case (nVar=4, full vectorized path) and a quick 5-iteration ONERAM6 3D run (nVar=5, exercises the vectorized-load remainder loop), GPU quantized-Jacobi matching the CPU reference exactly in both cases. Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysMatrixGPU.cu | 54 ++++++++++++++++------ 1 file changed, 41 insertions(+), 13 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 9c2d0fb0639..ddb024df7c6 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -26,6 +26,7 @@ */ #include +#include #include "../../include/linear_algebra/CMatrixInverse.hpp" #include "../../include/linear_algebra/CSysMatrix.inl" @@ -315,8 +316,12 @@ __global__ void IluBackwardKernel(const su2uint* __restrict__ level_idx, unsigne } /*! - * \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). + * \brief Block-LDU SpMV kernel: y[iRow] = (L + D + U) * x per block-row. Several rows are + * batched into one CUDA block (blockDim.x / nVar of them, see MatrixVectorProductGPU) + * instead of one row per block: nVar is typically ~4-6, so one-row-per-block leaves most + * of a warp's lanes permanently idle and caps occupancy at a few resident (mostly-empty) + * warps per SM, well before DRAM bandwidth is the limit. threadIdx.x indexes + * (row-within-block, output variable) as (threadIdx.x / nVar, threadIdx.x % nVar). */ template __global__ void BlockLDU_SpMV_kernel(unsigned long nRows, unsigned long nVar, @@ -328,9 +333,10 @@ __global__ void BlockLDU_SpMV_kernel(unsigned long nRows, unsigned long nVar, const su2uint* __restrict__ col_ind_u, const ScalarType* __restrict__ mat_u, const ScalarType* __restrict__ x, ScalarType* __restrict__ y) { - const unsigned long iRow = blockIdx.x; - const unsigned long iVar = threadIdx.x; - if (iRow >= nRows || iVar >= nVar) return; + const unsigned long rowsPerBlock = blockDim.x / nVar; + const unsigned long iVar = threadIdx.x % nVar; + const unsigned long iRow = static_cast(blockIdx.x) * rowsPerBlock + threadIdx.x / nVar; + if (iRow >= nRows) return; ScalarType sum = 0; /* Lower */ @@ -365,8 +371,12 @@ __device__ FORCEINLINE float DecodeQuantScaleDevice(int8_t e) { /*! * \brief Quantized block-LDU SpMV kernel: y[iRow] = (L + D + U) * x per block-row, reading int8 * row-scaled quantized blocks instead of full precision ones. Device version of - * QuantizedRowProduct/QuantizedMatVecAdd (CSysMatrix.inl). Same launch/thread layout as - * BlockLDU_SpMV_kernel: one CUDA block per block-row, threadIdx.x indexes output variable. + * QuantizedRowProduct/QuantizedMatVecAdd (CSysMatrix.inl). Rows are batched per block the + * same way as BlockLDU_SpMV_kernel (see its comment). Each row of quantized mantissas is + * also read 4 bytes at a time (one 32-bit load instead of four 8-bit ones): dp4a does not + * apply here since only the matrix is quantized, not x, so the accumulation itself is + * still one scalar FMA per element in the exact same order as the scalar loop (bit- + * identical result) - the win is fewer, wider load instructions, not fewer FMAs. */ template __global__ void QuantizedBlockLDU_SpMV_kernel( @@ -376,14 +386,27 @@ __global__ void QuantizedBlockLDU_SpMV_kernel( const int8_t* __restrict__ q_blocks_d, const su2uint* __restrict__ row_ptr_u, const su2uint* __restrict__ col_ind_u, const int8_t* __restrict__ q_scale_u, const int8_t* __restrict__ q_blocks_u, const ScalarType* __restrict__ x, ScalarType* __restrict__ y) { - const unsigned long iRow = blockIdx.x; - const unsigned long iVar = threadIdx.x; - if (iRow >= nRows || iVar >= nVar) return; + const unsigned long rowsPerBlock = blockDim.x / nVar; + const unsigned long iVar = threadIdx.x % nVar; + const unsigned long iRow = static_cast(blockIdx.x) * rowsPerBlock + threadIdx.x / nVar; + if (iRow >= nRows) return; auto addBlock = [&](const int8_t* __restrict__ qs, const int8_t* __restrict__ qv, const ScalarType* __restrict__ xk) { const float row_scale = DecodeQuantScaleDevice(qs[iVar]); + const int8_t* __restrict__ row = qv + iVar * nVar; ScalarType partial = 0; - for (unsigned long jVar = 0; jVar < nVar; ++jVar) partial += qv[iVar * nVar + jVar] * xk[jVar]; + unsigned long jVar = 0; + for (; jVar + 4 <= nVar; jVar += 4) { + /*--- Row bytes are not generally 4-byte aligned (nVar*nVar need not be a multiple of 4), + * so this must go through memcpy rather than a reinterpret_cast deref. ---*/ + uint32_t packed; + memcpy(&packed, row + jVar, sizeof(packed)); + partial += static_cast(packed) * xk[jVar]; + partial += static_cast(packed >> 8) * xk[jVar + 1]; + partial += static_cast(packed >> 16) * xk[jVar + 2]; + partial += static_cast(packed >> 24) * xk[jVar + 3]; + } + for (; jVar < nVar; ++jVar) partial += row[jVar] * xk[jVar]; return static_cast(row_scale) * partial; }; @@ -619,8 +642,13 @@ void CSysMatrix::MatrixVectorProductGPU(const CSysVector ScalarType* d_vec = vec.GetDevicePointer(); ScalarType* d_prod = prod.GetDevicePointer(); - dim3 blockDim(static_cast(nVar), 1, 1); - dim3 gridDim(static_cast(nPointDomain), 1, 1); + /*--- Batch several rows per block (see BlockLDU_SpMV_kernel's comment): nVar is small + * (typically ~4-6), so one row per block would leave most of a warp idle. Aim for ~128 + * threads/block, the largest whole number of rows that fits. ---*/ + constexpr unsigned long targetThreadsPerBlock = 128; + const auto rowsPerBlock = std::max(1, targetThreadsPerBlock / nVar); + dim3 blockDim(static_cast(rowsPerBlock * nVar), 1, 1); + dim3 gridDim(static_cast((nPointDomain + rowsPerBlock - 1) / rowsPerBlock), 1, 1); if (quantized_mode) { QuantizedBlockLDU_SpMV_kernel<<>>( nPointDomain, nVar, gpu.row_ptr_l, gpu.col_ind_l, d_q_scale.l, d_q_blocks.l, d_q_scale.d, d_q_blocks.d, From a113c70b32a14c76e4ddf22f9f7caea3a4d5a8d0 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 14 Aug 2026 23:00:44 -0700 Subject: [PATCH 09/30] Thread-per-row layout for the Jacobi preconditioner apply kernel ApplyJacobiPreconditionerKernel already had good occupancy (128 threads/block, one thread per point), so the SpMV row-batching fix doesn't apply here as-is. Its actual problem is coalescing: each thread read a whole nVar x nVar invM block serially, so at any given load instruction consecutive threads (consecutive points) were nVar^2 elements apart - e.g. 200 bytes for nVar=5, double - with essentially no sharing of the 32-byte memory sectors. Switches to the same thread mapping as the SpMV kernels (threadIdx.x -> (point-within-block, output variable) via divmod by nVar), so each thread computes one output row via its own dot-product loop rather than the whole block. This shrinks the inter-thread stride to nVar elements (a real, if partial, coalescing improvement - not the full stride-1 coalescing a one-thread-per-block-entry-plus-reduction version would get) while adding no shared memory or synchronization, since occupancy was never the problem here. Verified bit-identical to the pre-change kernel: NACA0012 (nVar=4) GPU Jacobi/Q_Jacobi runs reproduce their prior CL/CD exactly, and a 3D ONERAM6 (nVar=5) GPU Jacobi run matches its CPU reference exactly. Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysMatrixGPU.cu | 39 ++++++++++++++-------- 1 file changed, 25 insertions(+), 14 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index ddb024df7c6..9e121beefe6 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -34,23 +34,32 @@ namespace { +/*! + * \brief Apply the Jacobi preconditioner: prod = invM * vec (block-diagonal). Points are batched + * into blocks of ~128 threads (see ComputeJacobiPreconditionerGPU), threadIdx.x mapping + * to (point-within-block, output variable) via divmod by nVar - the same layout as + * BlockLDU_SpMV_kernel. Occupancy was already fine here (one thread per point, full + * warps), but consecutive threads used to land nVar^2 elements apart in invM (one point's + * whole dense block per thread); with this mapping they land nVar elements apart instead, + * a real (if partial, since it is still not stride-1) coalescing win that needs no shared + * memory or synchronization, unlike a fully-coalesced one-thread-per-block-entry version + * would. + */ template -__global__ void ApplyJacobiPreconditionerKernel(const ScalarType* invM, const ScalarType* vec, ScalarType* prod, +__global__ void ApplyJacobiPreconditionerKernel(const ScalarType* __restrict__ invM, + const ScalarType* __restrict__ vec, ScalarType* __restrict__ prod, unsigned long nPointDomain, unsigned long nVar) { - const auto iPoint = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + const unsigned long rowsPerBlock = blockDim.x / nVar; + const unsigned long iVar = threadIdx.x % nVar; + const unsigned long iPoint = static_cast(blockIdx.x) * rowsPerBlock + threadIdx.x / nVar; if (iPoint >= nPointDomain) return; - const auto block = &invM[iPoint * nVar * nVar]; - const auto rhs = &vec[iPoint * nVar]; - auto out = &prod[iPoint * nVar]; + const auto* block = &invM[iPoint * nVar * nVar + iVar * nVar]; + const auto* rhs = &vec[iPoint * nVar]; - for (auto iVar = 0ul; iVar < nVar; ++iVar) { - auto sum = ScalarType(0); - for (auto jVar = 0ul; jVar < nVar; ++jVar) { - sum += block[iVar * nVar + jVar] * rhs[jVar]; - } - out[iVar] = sum; - } + auto sum = ScalarType(0); + for (unsigned long jVar = 0; jVar < nVar; ++jVar) sum += block[jVar] * rhs[jVar]; + prod[iPoint * nVar + iVar] = sum; } /*--- ILU. The factorization is scheduled by coloring: colors are true independent sets of the @@ -441,8 +450,10 @@ void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector((nPointDomain + threadsPerBlock - 1) / threadsPerBlock); + constexpr unsigned long targetThreadsPerBlock = 128; + const auto rowsPerBlock = std::max(1, targetThreadsPerBlock / nVar); + const auto threadsPerBlock = static_cast(rowsPerBlock * nVar); + const auto blocks = static_cast((nPointDomain + rowsPerBlock - 1) / rowsPerBlock); ApplyJacobiPreconditionerKernel<<>>(d_invM, vec.GetDevicePointer(), prod.GetDevicePointer(), nPointDomain, nVar); /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ From 56fe054d9e89beafe799f2bbea45ddb4b7026092 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 08:10:41 -0700 Subject: [PATCH 10/30] Try thread-per-block-entry layout for the Jacobi apply kernel For completeness, compare against the thread-per-row version just committed: threadIdx.x now maps to (point-within-block, iVar, jVar), so consecutive threads read consecutive elements of invM (stride 1, fully coalesced across the whole block, not just within one point, since points are stored back to back) instead of landing nVar elements apart. Each thread computes one product into shared memory; the jVar==0 thread of each row then sums its nVar partials in the same order as the scalar loop it replaces (bit-identical result). Points are still batched per block for occupancy. The cost is a shared-memory write/read and one __syncthreads() per call, which the thread-per-row version didn't need - not obviously a win at this problem size without measuring, hence trying both. Verified bit-identical to the prior kernel: NACA0012 (nVar=4) GPU Jacobi/Q_Jacobi reproduce their known CL/CD exactly, and a 3D ONERAM6 (nVar=5) GPU Jacobi run matches its CPU reference exactly. Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysMatrixGPU.cu | 61 ++++++++++++++-------- 1 file changed, 38 insertions(+), 23 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 9e121beefe6..d21ea106376 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -35,31 +35,45 @@ namespace { /*! - * \brief Apply the Jacobi preconditioner: prod = invM * vec (block-diagonal). Points are batched - * into blocks of ~128 threads (see ComputeJacobiPreconditionerGPU), threadIdx.x mapping - * to (point-within-block, output variable) via divmod by nVar - the same layout as - * BlockLDU_SpMV_kernel. Occupancy was already fine here (one thread per point, full - * warps), but consecutive threads used to land nVar^2 elements apart in invM (one point's - * whole dense block per thread); with this mapping they land nVar elements apart instead, - * a real (if partial, since it is still not stride-1) coalescing win that needs no shared - * memory or synchronization, unlike a fully-coalesced one-thread-per-block-entry version - * would. + * \brief Apply the Jacobi preconditioner: prod = invM * vec (block-diagonal). One thread per + * block *entry* rather than per output row: threadIdx.x maps to (point-within-block, + * iVar, jVar) via divmod by nVar*nVar then nVar, so consecutive threads land at + * consecutive elements of invM (stride 1, fully coalesced - contiguous even across a + * point boundary, since points are stored back to back). Each thread computes one + * product and stages it in shared memory; the nVar threads of a row (fixed point, iVar) + * then get reduced by the jVar==0 thread, in the same j=0..nVar-1 order as the scalar + * loop this replaces, so the result is bit-identical. Points are still batched per block + * (see ComputeJacobiPreconditionerGPU) for occupancy, same as the row-per-thread version + * this replaces, but here the added shared-memory write/read and one __syncthreads() per + * call are the price paid for going from partial to full coalescing. */ template __global__ void ApplyJacobiPreconditionerKernel(const ScalarType* __restrict__ invM, const ScalarType* __restrict__ vec, ScalarType* __restrict__ prod, unsigned long nPointDomain, unsigned long nVar) { - const unsigned long rowsPerBlock = blockDim.x / nVar; - const unsigned long iVar = threadIdx.x % nVar; - const unsigned long iPoint = static_cast(blockIdx.x) * rowsPerBlock + threadIdx.x / nVar; - if (iPoint >= nPointDomain) return; + const unsigned long blockSize = nVar * nVar; + const unsigned long pointsPerBlock = blockDim.x / blockSize; + const unsigned long pointInBlock = threadIdx.x / blockSize; + const unsigned long localTid = threadIdx.x % blockSize; + const unsigned long iVar = localTid / nVar; + const unsigned long jVar = localTid % nVar; + const unsigned long iPoint = static_cast(blockIdx.x) * pointsPerBlock + pointInBlock; + + extern __shared__ __align__(sizeof(double)) char smem[]; + auto* partial = reinterpret_cast(smem); - const auto* block = &invM[iPoint * nVar * nVar + iVar * nVar]; - const auto* rhs = &vec[iPoint * nVar]; + ScalarType contribution = 0; + if (iPoint < nPointDomain) { + contribution = invM[iPoint * blockSize + localTid] * vec[iPoint * nVar + jVar]; + } + partial[threadIdx.x] = contribution; + __syncthreads(); - auto sum = ScalarType(0); - for (unsigned long jVar = 0; jVar < nVar; ++jVar) sum += block[jVar] * rhs[jVar]; - prod[iPoint * nVar + iVar] = sum; + if (iPoint < nPointDomain && jVar == 0) { + ScalarType sum = 0; + for (unsigned long j = 0; j < nVar; ++j) sum += partial[pointInBlock * blockSize + iVar * nVar + j]; + prod[iPoint * nVar + iVar] = sum; + } } /*--- ILU. The factorization is scheduled by coloring: colors are true independent sets of the @@ -451,11 +465,12 @@ void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector(1, targetThreadsPerBlock / nVar); - const auto threadsPerBlock = static_cast(rowsPerBlock * nVar); - const auto blocks = static_cast((nPointDomain + rowsPerBlock - 1) / rowsPerBlock); - ApplyJacobiPreconditionerKernel<<>>(d_invM, vec.GetDevicePointer(), prod.GetDevicePointer(), - nPointDomain, nVar); + const unsigned long blockSize = nVar * nVar; + const auto pointsPerBlock = std::max(1, targetThreadsPerBlock / blockSize); + const auto threadsPerBlock = static_cast(pointsPerBlock * blockSize); + const auto blocks = static_cast((nPointDomain + pointsPerBlock - 1) / pointsPerBlock); + ApplyJacobiPreconditionerKernel<<>>( + d_invM, vec.GetDevicePointer(), prod.GetDevicePointer(), nPointDomain, nVar); /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ gpuErrChk(cudaStreamSynchronize(nullptr)); gpuErrChk(cudaGetLastError()); From feeae49fcfc2f07e17ae32b084bee3def53b9aa3 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 08:58:47 -0700 Subject: [PATCH 11/30] Revert "Try thread-per-block-entry layout for the Jacobi apply kernel" This reverts commit 56fe054d9e89beafe799f2bbea45ddb4b7026092. --- Common/src/linear_algebra/CSysMatrixGPU.cu | 61 ++++++++-------------- 1 file changed, 23 insertions(+), 38 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index d21ea106376..9e121beefe6 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -35,45 +35,31 @@ namespace { /*! - * \brief Apply the Jacobi preconditioner: prod = invM * vec (block-diagonal). One thread per - * block *entry* rather than per output row: threadIdx.x maps to (point-within-block, - * iVar, jVar) via divmod by nVar*nVar then nVar, so consecutive threads land at - * consecutive elements of invM (stride 1, fully coalesced - contiguous even across a - * point boundary, since points are stored back to back). Each thread computes one - * product and stages it in shared memory; the nVar threads of a row (fixed point, iVar) - * then get reduced by the jVar==0 thread, in the same j=0..nVar-1 order as the scalar - * loop this replaces, so the result is bit-identical. Points are still batched per block - * (see ComputeJacobiPreconditionerGPU) for occupancy, same as the row-per-thread version - * this replaces, but here the added shared-memory write/read and one __syncthreads() per - * call are the price paid for going from partial to full coalescing. + * \brief Apply the Jacobi preconditioner: prod = invM * vec (block-diagonal). Points are batched + * into blocks of ~128 threads (see ComputeJacobiPreconditionerGPU), threadIdx.x mapping + * to (point-within-block, output variable) via divmod by nVar - the same layout as + * BlockLDU_SpMV_kernel. Occupancy was already fine here (one thread per point, full + * warps), but consecutive threads used to land nVar^2 elements apart in invM (one point's + * whole dense block per thread); with this mapping they land nVar elements apart instead, + * a real (if partial, since it is still not stride-1) coalescing win that needs no shared + * memory or synchronization, unlike a fully-coalesced one-thread-per-block-entry version + * would. */ template __global__ void ApplyJacobiPreconditionerKernel(const ScalarType* __restrict__ invM, const ScalarType* __restrict__ vec, ScalarType* __restrict__ prod, unsigned long nPointDomain, unsigned long nVar) { - const unsigned long blockSize = nVar * nVar; - const unsigned long pointsPerBlock = blockDim.x / blockSize; - const unsigned long pointInBlock = threadIdx.x / blockSize; - const unsigned long localTid = threadIdx.x % blockSize; - const unsigned long iVar = localTid / nVar; - const unsigned long jVar = localTid % nVar; - const unsigned long iPoint = static_cast(blockIdx.x) * pointsPerBlock + pointInBlock; - - extern __shared__ __align__(sizeof(double)) char smem[]; - auto* partial = reinterpret_cast(smem); + const unsigned long rowsPerBlock = blockDim.x / nVar; + const unsigned long iVar = threadIdx.x % nVar; + const unsigned long iPoint = static_cast(blockIdx.x) * rowsPerBlock + threadIdx.x / nVar; + if (iPoint >= nPointDomain) return; - ScalarType contribution = 0; - if (iPoint < nPointDomain) { - contribution = invM[iPoint * blockSize + localTid] * vec[iPoint * nVar + jVar]; - } - partial[threadIdx.x] = contribution; - __syncthreads(); + const auto* block = &invM[iPoint * nVar * nVar + iVar * nVar]; + const auto* rhs = &vec[iPoint * nVar]; - if (iPoint < nPointDomain && jVar == 0) { - ScalarType sum = 0; - for (unsigned long j = 0; j < nVar; ++j) sum += partial[pointInBlock * blockSize + iVar * nVar + j]; - prod[iPoint * nVar + iVar] = sum; - } + auto sum = ScalarType(0); + for (unsigned long jVar = 0; jVar < nVar; ++jVar) sum += block[jVar] * rhs[jVar]; + prod[iPoint * nVar + iVar] = sum; } /*--- ILU. The factorization is scheduled by coloring: colors are true independent sets of the @@ -465,12 +451,11 @@ void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector(1, targetThreadsPerBlock / blockSize); - const auto threadsPerBlock = static_cast(pointsPerBlock * blockSize); - const auto blocks = static_cast((nPointDomain + pointsPerBlock - 1) / pointsPerBlock); - ApplyJacobiPreconditionerKernel<<>>( - d_invM, vec.GetDevicePointer(), prod.GetDevicePointer(), nPointDomain, nVar); + const auto rowsPerBlock = std::max(1, targetThreadsPerBlock / nVar); + const auto threadsPerBlock = static_cast(rowsPerBlock * nVar); + const auto blocks = static_cast((nPointDomain + rowsPerBlock - 1) / rowsPerBlock); + ApplyJacobiPreconditionerKernel<<>>(d_invM, vec.GetDevicePointer(), prod.GetDevicePointer(), + nPointDomain, nVar); /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ gpuErrChk(cudaStreamSynchronize(nullptr)); gpuErrChk(cudaGetLastError()); From e65bfa8eb89f48ef7f891faaad2c232b2583eccf Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 09:04:35 -0700 Subject: [PATCH 12/30] Batch rows per block in the Jacobi diagonal-inversion kernel InvertDiagonalBlocksKernel launched one block per row with blockDim.x == nVar*nVar (16-36 threads for nVar=4-6), under a full warp, wasting lanes the same way the SpMV/apply kernels did before their fix. Batches ~128 threads worth of rows per block instead (threadIdx.x -> (row-within-block, block entry) via divmod by nVar*nVar), each row getting its own slice of the dynamic shared memory buffer. The inversion itself (SU2_LinAlg::MatrixInverse, shared host/device code) stays serial on one thread per row: it's a row-oriented Gauss-Jordan with a genuine sequential dependency chain of ~nVar^2/2 elimination substeps, so spreading it across threads would trade cheap serial FLOPs (~nVar^3, tiny at this size) for that many __syncthreads() barriers - likely a net loss, and not worth it anyway since this kernel only runs once per Build(), not per Krylov iteration. Batching just fixes the launch-shape waste, same as the other kernels. Verified bit-identical: NACA0012 (nVar=4) GPU Jacobi/Q_Jacobi reproduce their known CL/CD exactly, and a 3D ONERAM6 (nVar=5) GPU Jacobi run matches its CPU reference exactly. Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysMatrixGPU.cu | 44 ++++++++++++++-------- 1 file changed, 29 insertions(+), 15 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 9e121beefe6..aec6575b837 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -109,26 +109,37 @@ __device__ FORCEINLINE ScalarType* GetBlockILU(const DeviceLDU& M, u /*! * \brief Invert the diagonal blocks of the matrix, device version of InverseDiagonalBlock. - * \note Grid: one block per row, blockDim.x == nVar*nVar (one thread per block entry, they - * stage the input in shared memory). Dynamic shared memory: nVar*nVar scalars. + * Rows (points) are batched into blocks of ~128 threads (see + * BuildJacobiPreconditionerGPU), same reasoning as the SpMV/Jacobi-apply kernels: for + * small nVar a block of just nVar*nVar threads (16-36 for nVar=4-6) is less than a full + * warp, so one-row-per-block wastes lanes. The inversion itself (SU2_LinAlg::MatrixInverse, + * shared __host__ __device__ code) stays serial per row - it is a row-oriented Gauss- + * Jordan with a genuine sequential dependency chain across ~nVar^2/2 elimination + * substeps, so parallelizing it across threads would trade cheap serial FLOPs (~nVar^3, + * tiny for nVar<=~6) for that many __syncthreads() barriers, likely a net loss at this + * size - see the loader threads below, which only cooperate on the (fully parallel) load, + * not the inversion. threadIdx.x maps to (row-within-block, block entry) via divmod by + * nVar*nVar. Dynamic shared memory: blockDim.x scalars, one nVar*nVar work buffer per row + * in the block. */ template __global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, const ScalarType* __restrict__ mat_d, ScalarType* __restrict__ invM) { - const unsigned long iRow = blockIdx.x; - if (iRow >= nRows) return; - - const auto blockSize = nVar * nVar; - const unsigned long tid = threadIdx.x; - - /*--- The inversion destroys its input, so it cannot work on the matrix itself. ---*/ + const unsigned long blockSize = nVar * nVar; + const unsigned long rowsPerBlock = blockDim.x / blockSize; + const unsigned long rowInBlock = threadIdx.x / blockSize; + const unsigned long localTid = threadIdx.x % blockSize; + const unsigned long iRow = static_cast(blockIdx.x) * rowsPerBlock + rowInBlock; + + /*--- The inversion destroys its input, so it cannot work on the matrix itself; each row in the + * batch gets its own blockSize-sized slice of shared memory. ---*/ extern __shared__ __align__(sizeof(double)) char smem[]; - auto* work = reinterpret_cast(smem); + auto* work = reinterpret_cast(smem) + rowInBlock * blockSize; - work[tid] = mat_d[iRow * blockSize + tid]; + if (iRow < nRows) work[localTid] = mat_d[iRow * blockSize + localTid]; __syncthreads(); - if (tid == 0) SU2_LinAlg::MatrixInverse(nVar, work, invM + iRow * blockSize); + if (iRow < nRows && localTid == 0) SU2_LinAlg::MatrixInverse(nVar, work, invM + iRow * blockSize); } /*! @@ -472,10 +483,13 @@ void CSysMatrix::BuildJacobiPreconditionerGPU() { /*--- 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); + constexpr unsigned long targetThreadsPerBlock = 128; + const unsigned long blockSize = nVar * nVar; + const auto rowsPerBlock = std::max(1, targetThreadsPerBlock / blockSize); + const auto threadsPerBlock = static_cast(rowsPerBlock * blockSize); + const auto blocks = static_cast((nPointDomain + rowsPerBlock - 1) / rowsPerBlock); InvertDiagonalBlocksKernel - <<(nPointDomain), blockSize, blockSize * sizeof(ScalarType)>>>(nPointDomain, nVar, gpu.d, - d_invM); + <<>>(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()); From 7c308e6be8af932ecf197582201fbfda3c9969f8 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 09:11:51 -0700 Subject: [PATCH 13/30] Revert "Batch rows per block in the Jacobi diagonal-inversion kernel" This reverts commit e65bfa8eb89f48ef7f891faaad2c232b2583eccf. --- Common/src/linear_algebra/CSysMatrixGPU.cu | 44 ++++++++-------------- 1 file changed, 15 insertions(+), 29 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index aec6575b837..9e121beefe6 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -109,37 +109,26 @@ __device__ FORCEINLINE ScalarType* GetBlockILU(const DeviceLDU& M, u /*! * \brief Invert the diagonal blocks of the matrix, device version of InverseDiagonalBlock. - * Rows (points) are batched into blocks of ~128 threads (see - * BuildJacobiPreconditionerGPU), same reasoning as the SpMV/Jacobi-apply kernels: for - * small nVar a block of just nVar*nVar threads (16-36 for nVar=4-6) is less than a full - * warp, so one-row-per-block wastes lanes. The inversion itself (SU2_LinAlg::MatrixInverse, - * shared __host__ __device__ code) stays serial per row - it is a row-oriented Gauss- - * Jordan with a genuine sequential dependency chain across ~nVar^2/2 elimination - * substeps, so parallelizing it across threads would trade cheap serial FLOPs (~nVar^3, - * tiny for nVar<=~6) for that many __syncthreads() barriers, likely a net loss at this - * size - see the loader threads below, which only cooperate on the (fully parallel) load, - * not the inversion. threadIdx.x maps to (row-within-block, block entry) via divmod by - * nVar*nVar. Dynamic shared memory: blockDim.x scalars, one nVar*nVar work buffer per row - * in the block. + * \note Grid: one block per row, blockDim.x == nVar*nVar (one thread per block entry, they + * stage the input in shared memory). Dynamic shared memory: nVar*nVar scalars. */ template __global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, const ScalarType* __restrict__ mat_d, ScalarType* __restrict__ invM) { - const unsigned long blockSize = nVar * nVar; - const unsigned long rowsPerBlock = blockDim.x / blockSize; - const unsigned long rowInBlock = threadIdx.x / blockSize; - const unsigned long localTid = threadIdx.x % blockSize; - const unsigned long iRow = static_cast(blockIdx.x) * rowsPerBlock + rowInBlock; - - /*--- The inversion destroys its input, so it cannot work on the matrix itself; each row in the - * batch gets its own blockSize-sized slice of shared memory. ---*/ + const unsigned long iRow = blockIdx.x; + if (iRow >= nRows) return; + + const auto blockSize = nVar * nVar; + const unsigned long tid = threadIdx.x; + + /*--- The inversion destroys its input, so it cannot work on the matrix itself. ---*/ extern __shared__ __align__(sizeof(double)) char smem[]; - auto* work = reinterpret_cast(smem) + rowInBlock * blockSize; + auto* work = reinterpret_cast(smem); - if (iRow < nRows) work[localTid] = mat_d[iRow * blockSize + localTid]; + work[tid] = mat_d[iRow * blockSize + tid]; __syncthreads(); - if (iRow < nRows && localTid == 0) SU2_LinAlg::MatrixInverse(nVar, work, invM + iRow * blockSize); + if (tid == 0) SU2_LinAlg::MatrixInverse(nVar, work, invM + iRow * blockSize); } /*! @@ -483,13 +472,10 @@ void CSysMatrix::BuildJacobiPreconditionerGPU() { /*--- 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. ---*/ - constexpr unsigned long targetThreadsPerBlock = 128; - const unsigned long blockSize = nVar * nVar; - const auto rowsPerBlock = std::max(1, targetThreadsPerBlock / blockSize); - const auto threadsPerBlock = static_cast(rowsPerBlock * blockSize); - const auto blocks = static_cast((nPointDomain + rowsPerBlock - 1) / rowsPerBlock); + const auto blockSize = static_cast(nVar * nVar); InvertDiagonalBlocksKernel - <<>>(nPointDomain, nVar, gpu.d, d_invM); + <<(nPointDomain), blockSize, 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()); From 40007df65ce36af7c6a5c3051c82e040ed4f1c6a Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 09:44:25 -0700 Subject: [PATCH 14/30] Quantize the diagonal on the GPU straight from gpu.d, sharing the encoding logic gpu.d is already uploaded unconditionally by HtDTransfer() (Jacobi's own build needs the full precision diagonal regardless of quantization), so quantize it there directly under CUDA instead of computing on the host and uploading the int8 result - one less host/device round trip for the diagonal specifically (L/U still have to be quantized on the host, during assembly, since only the host ever touches the matrix as it is being built). To do this without duplicating the encoding algorithm (as previously happened with a hand-rolled DecodeQuantScaleDevice), DecodeQuantScale and EncodeQuantBlock (CSysMatrix.hpp) are now genuinely shared __host__ __device__ code, following the same SU2_CUDA_HOST_DEVICE pattern CMatrixInverse.hpp already established for MatrixInverse. That macro was previously defined locally (with a matching #undef at EOF) in three separate files - CMatrixInverse.hpp, CSysVector.hpp and vector_expressions.hpp - which do not compose: whichever of them is included last silently undefines it for anything appearing after, in the same translation unit. Centralized the single definition in code_config.hpp instead. Getting the shared quantization code to actually behave correctly under nvcc took two rounds of fixes, found by bisecting with a local, device-only reimplementation as a known-good reference: - memcpy-based float/uint32_t reinterpretation, used for the packed scale exponent, silently misinterprets the bits when compiled for the device (no diagnostic, just a wrong answer) - now branches on __CUDA_ARCH__ to use __float_as_uint/__uint_as_float on device. - std::max/std::min are not actually callable from device code without --expt-relaxed-constexpr (not used in this build); nvcc only warns ("calling a constexpr __host__ function ... is not allowed") rather than erroring, and silently emits broken device code. Added QuantMax/QuantMin, which branch the same way (CUDA's own unqualified max/min on device, std::max/std::min on host). Both bugs were fully silent (clean compile, no crash, no NaN) and produced a wrong-but-plausible-looking answer, so they were only caught by comparing against the known-correct CPU/pre-refactor GPU results - a compute-sanitizer racecheck run would have been the more direct tool but the toolchain in this environment (2022.4.1) could not attach to this MPI-launched binary. Verified against CPU/pre-refactor references: NACA0012 (nVar=4) GPU Jacobi/Q_Jacobi/Q_Identity all reproduce their known CL/CD exactly, and a 3D ONERAM6 (nVar=5) GPU Q_Jacobi run matches its CPU counterpart exactly. Co-Authored-By: Claude Sonnet 5 --- Common/include/code_config.hpp | 10 ++ .../include/linear_algebra/CMatrixInverse.hpp | 8 +- Common/include/linear_algebra/CSysMatrix.hpp | 99 +++++++++++++++---- Common/include/linear_algebra/CSysVector.hpp | 5 +- .../linear_algebra/vector_expressions.hpp | 8 +- Common/src/linear_algebra/CSysMatrix.cpp | 71 +++++++------ Common/src/linear_algebra/CSysMatrixGPU.cu | 71 ++++++++----- 7 files changed, 177 insertions(+), 95 deletions(-) diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp index 41d3c747cf8..93c82eaad88 100644 --- a/Common/include/code_config.hpp +++ b/Common/include/code_config.hpp @@ -50,6 +50,16 @@ #define NEVERINLINE inline #endif +/*--- Marks a function callable from both host and device code (e.g. a small numerical kernel + * shared between a CSysMatrix host routine and its CSysMatrixGPU.cu counterpart), a no-op outside + * nvcc. Centralized here (rather than defined locally, as e.g. CMatrixInverse.hpp used to) so + * every shared routine uses the exact same macro instead of independent copies drifting apart. ---*/ +#ifdef __CUDACC__ +#define SU2_CUDA_HOST_DEVICE __host__ __device__ +#else +#define SU2_CUDA_HOST_DEVICE +#endif + #if defined(__INTEL_COMPILER) /*--- Disable warnings related to inline attributes. ---*/ #pragma warning disable 2196 diff --git a/Common/include/linear_algebra/CMatrixInverse.hpp b/Common/include/linear_algebra/CMatrixInverse.hpp index f3d64630349..6d6dfa73959 100644 --- a/Common/include/linear_algebra/CMatrixInverse.hpp +++ b/Common/include/linear_algebra/CMatrixInverse.hpp @@ -30,11 +30,7 @@ #include -#ifdef __CUDACC__ -#define SU2_CUDA_HOST_DEVICE __host__ __device__ -#else -#define SU2_CUDA_HOST_DEVICE -#endif +#include "../code_config.hpp" namespace SU2_LinAlg { @@ -95,5 +91,3 @@ SU2_CUDA_HOST_DEVICE inline void MatrixInverse(unsigned long nVar, ScalarType* m } } // namespace SU2_LinAlg - -#undef SU2_CUDA_HOST_DEVICE diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 95db350def8..75df44b94a9 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -29,6 +29,7 @@ #pragma once #include "../CConfig.hpp" +#include "../code_config.hpp" #include "CSysVector.hpp" #include "CPastixWrapper.hpp" #include "../toolboxes/graph_toolbox.hpp" @@ -109,44 +110,94 @@ struct CSysMatrixComms { MPI_QUANTITIES commType = MPI_QUANTITIES::SOLUTION_MATRIX); }; +/*! + * \brief std::max/std::min, usable from both host and device code. Calling std::max/std::min + * directly from a SU2_CUDA_HOST_DEVICE function compiles without error but is not + * actually valid without --expt-relaxed-constexpr (not used in this build) - nvcc only + * warns ("calling a constexpr __host__ function ... is not allowed"), then silently + * emits device code that does not do what it looks like it does. CUDA's device compiler + * provides its own (unqualified) max/min built-ins instead, which these two forward to; + * on the host they just forward to std::max/std::min. + */ +template +SU2_CUDA_HOST_DEVICE inline T QuantMax(T a, T b) noexcept { +#ifdef __CUDA_ARCH__ + return max(a, b); +#else + return std::max(a, b); +#endif +} +template +SU2_CUDA_HOST_DEVICE inline T QuantMin(T a, T b) noexcept { +#ifdef __CUDA_ARCH__ + return min(a, b); +#else + return std::min(a, b); +#endif +} + /*! * \brief Reconstruct the float row-scale from a stored int8 binary exponent. * The exponent \p e was packed as (e + 127) into the IEEE 754 biased-exponent field * with a zero mantissa, giving an exact power of two: 2^e. * This is the inverse of the encoding in EncodeQuantBlock. + * \note Shared verbatim with the device (the quantized SpMV kernel in CSysMatrixGPU.cu decodes + * through this same function, not a separate copy), see SU2_CUDA_HOST_DEVICE. The bit + * reinterpretation itself branches on __CUDA_ARCH__ (mirroring RegularizePivot in + * CMatrixInverse.hpp): __uint_as_float on device, memcpy on host - plain memcpy compiles + * for the device too, but was observed to silently misinterpret the bits there rather than + * reinterpreting them. */ -FORCEINLINE float DecodeQuantScale(int8_t e) noexcept { - const uint32_t bits = static_cast(std::max(0, static_cast(e) + 127)) << 23; +SU2_CUDA_HOST_DEVICE inline float DecodeQuantScale(int8_t e) noexcept { + const uint32_t bits = static_cast(QuantMax(0, static_cast(e) + 127)) << 23; +#ifdef __CUDA_ARCH__ + return __uint_as_float(bits); +#else float scale; memcpy(&scale, &bits, sizeof(bits)); return scale; +#endif } /*! * \brief Encode one nVar×nVar block into per-row int8 quantized storage. * \p f(r,c) is called twice per entry (max-abs scan then encoding); it should be cheap. * Stores a per-row scale exponent in \p qs and clamped int8 values in \p qv. + * \note Shared verbatim with the device (CSysMatrixGPU.cu's diagonal-quantization kernel encodes + * through this same function, not a separate copy), see SU2_CUDA_HOST_DEVICE. \p f must + * already return a passive (non-AD) value: every caller does, since quantization is + * compile-time disabled whenever ScalarType could be AD-active (see quantized_mode), so + * there is no SU2_TYPE::PassiveValue call here to keep this device-callable. The bit + * reinterpretations branch on __CUDA_ARCH__ the same way DecodeQuantScale does, see there. */ template -FORCEINLINE void EncodeQuantBlock(const F& f, int8_t* __restrict qs, int8_t* __restrict qv, - unsigned long nVar) noexcept { +SU2_CUDA_HOST_DEVICE inline void EncodeQuantBlock(const F& f, int8_t* __restrict qs, int8_t* __restrict qv, + unsigned long nVar) noexcept { for (auto r = 0ul; r < nVar; ++r) { constexpr uint32_t eps_bits = 0x34000000u; uint32_t max_abs_bits = eps_bits; for (auto c = 0ul; c < nVar; ++c) { - const float fv = SU2_TYPE::PassiveValue(f(r, c)); + const float fv = static_cast(f(r, c)); +#ifdef __CUDA_ARCH__ + const uint32_t fb = __float_as_uint(fv); +#else uint32_t fb; memcpy(&fb, &fv, sizeof(fb)); - max_abs_bits = std::max(max_abs_bits, fb & 0x7FFFFFFFu); +#endif + max_abs_bits = QuantMax(max_abs_bits, fb & 0x7FFFFFFFu); } - const int e = std::min(127, std::max(-128, static_cast(max_abs_bits >> 23) - 133)); + const int e = QuantMin(127, QuantMax(-128, static_cast(max_abs_bits >> 23) - 133)); qs[r] = static_cast(e); const uint32_t inv_bits = static_cast(127 - e) << 23; +#ifdef __CUDA_ARCH__ + const float inv_rscale = __uint_as_float(inv_bits); +#else float inv_rscale; memcpy(&inv_rscale, &inv_bits, sizeof(inv_rscale)); +#endif for (auto c = 0ul; c < nVar; ++c) { qv[r * nVar + c] = - static_cast(std::max(-128.f, std::min(127.f, roundf(SU2_TYPE::PassiveValue(f(r, c)) * inv_rscale)))); + static_cast(QuantMax(-128.f, QuantMin(127.f, roundf(static_cast(f(r, c)) * inv_rscale)))); } } } @@ -272,22 +323,27 @@ class CSysMatrix { #else static constexpr bool quantized_mode = false; #endif - /*!< \brief Per-row exponents; .l/.u sized [nnz_l/u * nVar], .d [nPoint * nVar] - * (populated by QuantizeDiagonalBlocks(), always on the host, see below). - * Pinned (cudaMallocHost) rather than aligned_alloc when useCuda, see Initialize(), - * so the async uploads below are genuinely asynchronous instead of silently - * blocking (cudaMemcpyAsync only overlaps with the host from pinned memory). */ + /*!< \brief Per-row exponents; .l/.u sized [nnz_l/u * nVar], .d [nPoint * nVar]. .l/.u are + * populated during assembly (quantized on the fly); .d is populated by + * QuantizeDiagonalBlocks() on the host, but only when !useCuda - under CUDA the + * diagonal is quantized straight from gpu.d instead (QuantizeDiagonalBlocksGPU()), + * so q_scale.d/q_blocks.d are simply unused there. .l/.u are pinned (cudaMallocHost) + * rather than aligned_alloc when useCuda, see Initialize(), so HtDTransfer()'s async + * uploads of them are genuinely asynchronous instead of silently blocking + * (cudaMemcpyAsync only overlaps with the host from pinned memory); .d is never + * uploaded, so it is always plain aligned_alloc regardless of useCuda. */ LDU q_scale; /*!< \brief Quantized block entries; .l/.u sized [nnz_l/u * nVar * nEqn], .d [nPoint * nVar * nEqn]. */ LDU q_blocks; /*!< \brief Device mirrors of the quantized storage, only allocated when quantized_mode && * useCuda (currently only reachable for Q_JACOBI/Q_IDENTITY, Q_LU_SGS stays host-only). - * d_q_scale.l/.u and d_q_blocks.l/.u are plain device-side copies of q_scale.l/.u and - * q_blocks.l/.u, uploaded *asynchronously* by HtDTransfer() so that transfer can overlap with - * the host quantizing the diagonal; d_q_scale.d/d_q_blocks.d are that host result, likewise - * uploaded asynchronously (a plain cudaMemcpyAsync, not a kernel), at the end of - * QuantizeDiagonalBlocks(). */ + * d_q_scale.l/.u and d_q_blocks.l/.u are device-side copies of q_scale.l/.u and q_blocks.l/.u, + * uploaded *asynchronously* by HtDTransfer(). d_q_scale.d/d_q_blocks.d are populated directly + * on the device from gpu.d by QuantizeDiagonalBlocksGPU() instead - gpu.d is uploaded + * unconditionally by HtDTransfer() anyway (Jacobi's own build needs the full precision + * diagonal regardless of quantization), so quantizing it there avoids a host quantize + + * upload round trip for the diagonal specifically. */ LDU d_q_scale; LDU d_q_blocks; @@ -573,6 +629,13 @@ class CSysMatrix { void MatrixVectorProductGPU(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const; + /*! + * \brief Quantize the diagonal blocks directly on the device, from gpu.d into + * d_q_scale.d/d_q_blocks.d (EncodeQuantBlock, shared with the host path). + * \note Requires the device matrix to be up to date, see HtDTransfer. + */ + void QuantizeDiagonalBlocksGPU(); + /*! * \brief Build the Jacobi preconditioner on the device, from the device copy of the matrix. * \note Requires the device matrix to be up to date, see HtDTransfer. diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 1b57a6ffe94..87367f5f151 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -31,6 +31,7 @@ #include #include +#include "../code_config.hpp" #include "../parallelization/mpi_structure.hpp" #include "../parallelization/omp_structure.hpp" #include "../parallelization/vectorization.hpp" @@ -39,9 +40,6 @@ #ifdef __CUDACC__ #include "GPUComms.cuh" -#define SU2_CUDA_HOST_DEVICE __host__ __device__ -#else -#define SU2_CUDA_HOST_DEVICE #endif template @@ -709,4 +707,3 @@ CVectorView::CVectorView(const CSysVector& vector) #undef CSYSVEC_PARFOR #undef END_CSYSVEC_PARFOR -#undef SU2_CUDA_HOST_DEVICE diff --git a/Common/include/linear_algebra/vector_expressions.hpp b/Common/include/linear_algebra/vector_expressions.hpp index f9b941a0bdb..d59dfdf51ae 100644 --- a/Common/include/linear_algebra/vector_expressions.hpp +++ b/Common/include/linear_algebra/vector_expressions.hpp @@ -28,6 +28,7 @@ #pragma once #include "../basic_types/datatype_structure.hpp" +#include "../code_config.hpp" #include #include #include @@ -39,12 +40,6 @@ namespace VecExpr { /// \addtogroup VecExpr /// @{ -#ifdef __CUDACC__ -#define SU2_CUDA_HOST_DEVICE __host__ __device__ -#else -#define SU2_CUDA_HOST_DEVICE -#endif - /*! * \brief Base vector expression class. * \ingroup BLAS @@ -251,5 +246,4 @@ MAKE_BINARY_FUN(operator>, gt_, gt_impl) #undef MAKE_BINARY_FUN /// @} -#undef SU2_CUDA_HOST_DEVICE } // namespace VecExpr diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index a4c52d1e134..7b3e87ebdaf 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -127,15 +127,15 @@ CSysMatrix::~CSysMatrix() { freeHostLDU(ilu); MemoryAllocation::aligned_free(invM); - /*--- q_scale/q_blocks are pinned (cudaMallocHost) rather than aligned_alloc when useCuda, see - * the comment in Initialize(); free with the matching deallocator either way. ---*/ + /*--- q_scale/q_blocks' .l/.u are pinned (cudaMallocHost) rather than aligned_alloc when + * useCuda, .d never is; see the comment in Initialize(). Free each with its matching + * deallocator. ---*/ auto freeQuantLDU = [this](auto& m) { + MemoryAllocation::aligned_free(m.d); if (useCuda) { - GPUMemoryAllocation::pinned_free(m.d); GPUMemoryAllocation::pinned_free(m.l); GPUMemoryAllocation::pinned_free(m.u); } else { - MemoryAllocation::aligned_free(m.d); MemoryAllocation::aligned_free(m.l); MemoryAllocation::aligned_free(m.u); } @@ -261,21 +261,26 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi #ifndef CODI_REVERSE_TYPE quantized_mode = true; #endif - /*--- Pinned (page-locked) when useCuda: HtDTransfer()/QuantizeDiagonalBlocks() upload these - * with cudaMemcpyAsync, which is only genuinely asynchronous from pinned host memory (see + /*--- .l/.u are pinned (page-locked) when useCuda: HtDTransfer() uploads them with + * cudaMemcpyAsync, which is only genuinely asynchronous from pinned host memory (see * GPUMemoryAllocation::pinned_alloc); from regular pageable memory it silently degrades to a - * blocking copy, defeating the overlap with host-side diagonal quantization. ---*/ - auto allocQ = [useCuda = this->useCuda](QuantType*& ptr, unsigned long n) { + * blocking copy. The diagonal (.d) is never uploaded - under CUDA it is quantized straight + * from gpu.d on the device instead (QuantizeDiagonalBlocksGPU()) - so it stays plain + * aligned_alloc regardless of useCuda. ---*/ + auto allocQ = [](QuantType*& ptr, unsigned long n) { + ptr = MemoryAllocation::aligned_alloc(64, n * sizeof(QuantType)); + }; + auto allocPinnedIfCuda = [useCuda = this->useCuda](QuantType*& ptr, unsigned long n) { if (useCuda) { ptr = GPUMemoryAllocation::pinned_alloc(n * sizeof(QuantType)); } else { ptr = MemoryAllocation::aligned_alloc(64, n * sizeof(QuantType)); } }; - allocQ(q_scale.l, mat.nnz_l * nVar); - allocQ(q_blocks.l, mat.nnz_l * nVar * nEqn); - allocQ(q_scale.u, mat.nnz_u * nVar); - allocQ(q_blocks.u, mat.nnz_u * nVar * nEqn); + allocPinnedIfCuda(q_scale.l, mat.nnz_l * nVar); + allocPinnedIfCuda(q_blocks.l, mat.nnz_l * nVar * nEqn); + allocPinnedIfCuda(q_scale.u, mat.nnz_u * nVar); + allocPinnedIfCuda(q_blocks.u, mat.nnz_u * nVar * nEqn); allocQ(q_scale.d, nPoint * nVar); allocQ(q_blocks.d, nPoint * nVar * nEqn); } else { @@ -770,34 +775,19 @@ void CSysMatrix::QuantizeDiagonalBlocks() { if (!quantized_mode) return; - /*--- Q_LU_SGS / Q_JACOBI / Q_IDENTITY: L/U were quantized during assembly; only the diagonal needs quantization - * now. Always done on the host, even under CUDA: HtDTransfer() already kicked off the (larger) - * L/U quantized transfer asynchronously before Build() reached this point, so quantizing the - * diagonal here on the CPU overlaps with that transfer instead of waiting on it first. ---*/ - SU2_OMP_FOR_DYN(omp_heavy_size) - for (auto i = 0ul; i < nPointDomain; ++i) - QuantizeBlock(&mat.d[i * nVar * nVar], &q_scale.d[i * nVar], &q_blocks.d[i * nVar * nVar]); - END_SU2_OMP_FOR - if (useCuda) { #ifdef SU2_ENABLE_CUDA_KERNELS if constexpr (su2_gpu_capable_v) { - /*--- Just an upload of the host result above, no computation, so a plain CUDA runtime call - * (available here via GPUComms.cuh, transitively included through allocation_toolbox.hpp) - * rather than a kernel dispatched through a CSysMatrixGPU.cu hook. Async (q_scale.d/ - * q_blocks.d are pinned, see Initialize()) so issuing it does not block the host; the sync - * right after does not undo that overlap; by this point the host has already spent the - * whole diagonal-quantization loop above letting the earlier (larger) L/U transfer, kicked - * off by HtDTransfer(), drain in the background; it only makes sure that transfer and this - * one are actually finished before Build() returns, matching every other GPU-touching - * function in this file (they all sync at the end, see the comment on those calls). ---*/ - BEGIN_SU2_DEVICE_REGION - gpuErrChk( - cudaMemcpyAsync(d_q_scale.d, q_scale.d, sizeof(QuantType) * nPointDomain * nVar, cudaMemcpyHostToDevice)); - gpuErrChk(cudaMemcpyAsync(d_q_blocks.d, q_blocks.d, sizeof(QuantType) * nPointDomain * nVar * nVar, - cudaMemcpyHostToDevice)); - gpuErrChk(cudaStreamSynchronize(nullptr)); - END_SU2_DEVICE_REGION + /*--- gpu.d is already on the device - HtDTransfer() uploads it unconditionally, since + * Jacobi's own build needs the full precision diagonal regardless of quantization - so + * quantize straight from it here instead of quantizing on the host and uploading the + * result: one less host/device round trip for the diagonal specifically (L/U still have + * to be quantized on the host, during assembly, since only the host ever touches the + * matrix as it is being built). This calls the exact same EncodeQuantBlock routine the + * host path below uses (CSysMatrix.hpp, SU2_CUDA_HOST_DEVICE), not a separate device copy. + * ---*/ + SU2_DEVICE_REGION(QuantizeDiagonalBlocksGPU();) + return; } else { GPUNotAvailable(CURRENT_FUNCTION); } @@ -805,6 +795,13 @@ void CSysMatrix::QuantizeDiagonalBlocks() { GPUNotAvailable(CURRENT_FUNCTION); #endif } + + /*--- Q_LU_SGS / Q_JACOBI / Q_IDENTITY: L/U were quantized during assembly; only the diagonal needs quantization + * now. ---*/ + SU2_OMP_FOR_DYN(omp_heavy_size) + for (auto i = 0ul; i < nPointDomain; ++i) + QuantizeBlock(&mat.d[i * nVar * nVar], &q_scale.d[i * nVar], &q_blocks.d[i * nVar * nVar]); + END_SU2_OMP_FOR } template diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 9e121beefe6..2fd81aca048 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -131,6 +131,25 @@ __global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nV if (tid == 0) SU2_LinAlg::MatrixInverse(nVar, work, invM + iRow * blockSize); } +/*! + * \brief Quantize the diagonal blocks straight from the device diagonal (gpu.d), device + * counterpart of CSysMatrix::QuantizeBlock (CSysMatrix.cpp), applied row by row. Calls the + * exact same EncodeQuantBlock encoding routine (CSysMatrix.hpp, SU2_CUDA_HOST_DEVICE) that + * the host path uses, rather than a separate device copy of the encoding logic. One + * thread per row. + */ +template +__global__ void QuantizeDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, + const ScalarType* __restrict__ mat_d, int8_t* __restrict__ q_scale_d, + int8_t* __restrict__ q_blocks_d) { + const auto iRow = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (iRow >= nRows) return; + + const auto* blk = mat_d + iRow * nVar * nVar; + EncodeQuantBlock([&](unsigned long r, unsigned long c) { return blk[r * nVar + c]; }, q_scale_d + iRow * nVar, + q_blocks_d + iRow * nVar * nVar, nVar); +} + /*! * \brief Factorize the rows of one color, one sweep of an iterative (colored Gauss-Seidel) * ILU factorization: same order/pattern as the exact level-scheduled algorithm (this @@ -368,15 +387,6 @@ __global__ void BlockLDU_SpMV_kernel(unsigned long nRows, unsigned long nVar, y[iRow * nVar + iVar] = sum; } -/*! - * \brief Device counterpart of DecodeQuantScale (CSysMatrix.hpp), bit-identical: reconstructs - * the float row-scale 2^e from its packed int8 binary exponent. - */ -__device__ FORCEINLINE float DecodeQuantScaleDevice(int8_t e) { - const uint32_t bits = static_cast(max(0, static_cast(e) + 127)) << 23; - return __uint_as_float(bits); -} - /*! * \brief Quantized block-LDU SpMV kernel: y[iRow] = (L + D + U) * x per block-row, reading int8 * row-scaled quantized blocks instead of full precision ones. Device version of @@ -401,7 +411,7 @@ __global__ void QuantizedBlockLDU_SpMV_kernel( if (iRow >= nRows) return; auto addBlock = [&](const int8_t* __restrict__ qs, const int8_t* __restrict__ qv, const ScalarType* __restrict__ xk) { - const float row_scale = DecodeQuantScaleDevice(qs[iVar]); + const float row_scale = DecodeQuantScale(qs[iVar]); const int8_t* __restrict__ row = qv + iVar * nVar; ScalarType partial = 0; unsigned long jVar = 0; @@ -461,6 +471,23 @@ void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector +void CSysMatrix::QuantizeDiagonalBlocksGPU() { + SU2_ZONE_SCOPED + + 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. ---*/ + constexpr unsigned threadsPerBlock = 128; + const auto blocks = static_cast((nPointDomain + threadsPerBlock - 1) / threadsPerBlock); + QuantizeDiagonalBlocksKernel + <<>>(nPointDomain, nVar, gpu.d, d_q_scale.d, d_q_blocks.d); + /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ + gpuErrChk(cudaStreamSynchronize(nullptr)); + gpuErrChk(cudaGetLastError()); +} + template void CSysMatrix::BuildJacobiPreconditionerGPU() { SU2_ZONE_SCOPED @@ -621,18 +648,17 @@ void CSysMatrix::HtDTransfer(bool trigger) const { gpuErrChk(cudaMemcpy(gpu.d, mat.d, sizeof(ScalarType) * nPoint * nVar * nEqn, cudaMemcpyHostToDevice)); if (quantized_mode) { /*--- No gpu.l/gpu.u to transfer (never allocated); mirror the host quantized off-diagonal - * storage instead. Issued as async copies on the default stream, then left in flight: the - * caller (CSysMatrixVectorProduct's constructor) returns right after this, and QuantizeDiag- - * onalBlocks() -> Build() quantizes the diagonal on the host next, before anything is - * launched on the device again. This is only genuinely asynchronous (i.e. the host thread - * does not block here waiting for the copy) because q_scale.l/q_blocks.l/q_scale.u/ - * q_blocks.u are pinned host memory, see the comment on those members / Initialize() - - * cudaMemcpyAsync silently degrades to a blocking copy from regular pageable memory. Any - * later kernel that reads d_q_scale/d_q_blocks is issued on the same default stream too, so - * strictly it would not need to wait for these explicitly; QuantizeDiagonalBlocks() still - * syncs right after uploading the diagonal mirrors (d_q_scale.d/d_q_blocks.d, the same way, - * async from pinned memory) to be safe and consistent with every other GPU-touching function - * in this file, all of which sync before returning. ---*/ + * storage instead (the diagonal mirrors, d_q_scale.d/d_q_blocks.d, are not touched here at + * all - QuantizeDiagonalBlocksGPU() populates them straight from gpu.d, just uploaded above, + * with no host round trip). Issued as async copies on the default stream, then left in + * flight: the caller (CSysMatrixVectorProduct's constructor) returns right after this, and + * whatever the preconditioner's Build() does next runs concurrently with the transfer still + * draining. This is only genuinely asynchronous (i.e. the host thread does not block here + * waiting for the copy) because q_scale.l/q_blocks.l/q_scale.u/q_blocks.u are pinned host + * memory, see the comment on those members / Initialize() - cudaMemcpyAsync silently + * degrades to a blocking copy from regular pageable memory. Any later kernel that reads + * d_q_scale.l/.u/d_q_blocks.l/.u is issued on the same default stream too, so it correctly + * waits for these by stream ordering without an explicit sync here. ---*/ gpuErrChk(cudaMemcpyAsync(d_q_scale.l, q_scale.l, sizeof(QuantType) * mat.nnz_l * nVar, cudaMemcpyHostToDevice)); gpuErrChk(cudaMemcpyAsync(d_q_blocks.l, q_blocks.l, sizeof(QuantType) * mat.nnz_l * nVar * nEqn, cudaMemcpyHostToDevice)); @@ -680,6 +706,7 @@ template void CSysMatrix::MatrixVectorProductGPU(const CSysVector& v CSysVector& prod, \ CGeometry* geometry, \ const CConfig* config) const; \ +template void CSysMatrix::QuantizeDiagonalBlocksGPU(); \ template void CSysMatrix::BuildJacobiPreconditionerGPU(); \ template void CSysMatrix::BuildILUPreconditionerGPU(); \ template void CSysMatrix::ComputeILUPreconditionerGPU(const CSysVector& vec, \ From 0bfd05baaee1f3e56385e54ef431dee3d945f573 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 15:18:50 -0700 Subject: [PATCH 15/30] Warp-align the row batching in the Jacobi diagonal-inversion kernel The reverted batching attempt (e65bfa8eb8) floor-divided a ~128 thread target by nVar*nVar to pick rows-per-block, which for nVar values that don't divide 32 evenly (e.g. nVar=5, blockSize=25 -> 125 threads) left the block's last warp partly idle. Round rowsPerBlock up instead to the smallest multiple that makes blockDim.x itself a multiple of 32. This does not touch the actual bottleneck (MatrixInverse is a serial per-row Gauss-Jordan, so only one thread per row ever does real work regardless of alignment) - it only tidies up the parallel load-into- shared-memory phase. Kept for the learning value of confirming that, even warp-aligned, batching does not help here. Verified bit-identical: NACA0012 (nVar=4) GPU Jacobi/Q_Jacobi/Q_Identity reproduce their known CL/CD exactly, and Q_Jacobi on ONERAM6 (nVar=5, which now batches 32 rows/800 threads per block instead of 5 rows/125 threads) matches its CPU reference to 5 decimals. Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysMatrixGPU.cu | 56 +++++++++++++++------- 1 file changed, 40 insertions(+), 16 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 2fd81aca048..bcbd05d79ff 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -27,6 +27,7 @@ #include #include +#include #include "../../include/linear_algebra/CMatrixInverse.hpp" #include "../../include/linear_algebra/CSysMatrix.inl" @@ -109,26 +110,39 @@ __device__ FORCEINLINE ScalarType* GetBlockILU(const DeviceLDU& M, u /*! * \brief Invert the diagonal blocks of the matrix, device version of InverseDiagonalBlock. - * \note Grid: one block per row, blockDim.x == nVar*nVar (one thread per block entry, they - * stage the input in shared memory). Dynamic shared memory: nVar*nVar scalars. + * Rows (points) are batched into blocks whose thread count is a whole multiple of the + * warp size (see BuildJacobiPreconditionerGPU), unlike the previous batched attempt + * (reverted: e65bfa8eb8) which just floor-divided a ~128 thread target by nVar*nVar - for + * nVar not dividing 32 evenly (e.g. nVar=5, blockSize=25 -> 125 threads) that left the + * block's last warp partly empty. The inversion itself (SU2_LinAlg::MatrixInverse, shared + * __host__ __device__ code) still runs serially on one thread per row - it is a row- + * oriented Gauss-Jordan with a genuine sequential dependency chain across ~nVar^2/2 + * elimination substeps, so spreading it across threads would trade cheap serial FLOPs + * (~nVar^3, tiny at this size) for that many __syncthreads() barriers. Warp-alignment does + * not change that: it only avoids wasting the last warp's tail lanes on the (fully + * parallel) load into shared memory, the inversion itself still runs one point at a time + * per warp regardless of alignment. threadIdx.x maps to (row-within-block, block entry) + * via divmod by nVar*nVar. Dynamic shared memory: blockDim.x scalars, one nVar*nVar work + * buffer per row in the block. */ template __global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, const ScalarType* __restrict__ mat_d, ScalarType* __restrict__ invM) { - const unsigned long iRow = blockIdx.x; - if (iRow >= nRows) return; - - const auto blockSize = nVar * nVar; - const unsigned long tid = threadIdx.x; - - /*--- The inversion destroys its input, so it cannot work on the matrix itself. ---*/ + const unsigned long blockSize = nVar * nVar; + const unsigned long rowsPerBlock = blockDim.x / blockSize; + const unsigned long rowInBlock = threadIdx.x / blockSize; + const unsigned long localTid = threadIdx.x % blockSize; + const unsigned long iRow = static_cast(blockIdx.x) * rowsPerBlock + rowInBlock; + + /*--- The inversion destroys its input, so it cannot work on the matrix itself; each row in the + * batch gets its own blockSize-sized slice of shared memory. ---*/ extern __shared__ __align__(sizeof(double)) char smem[]; - auto* work = reinterpret_cast(smem); + auto* work = reinterpret_cast(smem) + rowInBlock * blockSize; - work[tid] = mat_d[iRow * blockSize + tid]; + if (iRow < nRows) work[localTid] = mat_d[iRow * blockSize + localTid]; __syncthreads(); - if (tid == 0) SU2_LinAlg::MatrixInverse(nVar, work, invM + iRow * blockSize); + if (iRow < nRows && localTid == 0) SU2_LinAlg::MatrixInverse(nVar, work, invM + iRow * blockSize); } /*! @@ -498,11 +512,21 @@ void CSysMatrix::BuildJacobiPreconditionerGPU() { 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); + * CSysMatrixVectorProduct, which is created before the preconditioner is built. Round the + * number of rows batched per block up to the smallest multiple that makes blockDim.x itself a + * multiple of the warp size (32), so the block's last warp is never partly idle; then pick the + * smallest such multiple at or above the ~128 thread target instead of always taking the very + * first one, so nVar values that need a huge single step to reach alignment (e.g. nVar=5, + * blockSize=25 -> step of 32 rows/800 threads) don't blow the block far past the target. ---*/ + constexpr unsigned long targetThreadsPerBlock = 128; + const unsigned long blockSize = nVar * nVar; + const unsigned long warpAlignStep = 32 / std::gcd(32ul, blockSize); + auto rowsPerBlock = std::max(warpAlignStep, targetThreadsPerBlock / blockSize); + rowsPerBlock = ((rowsPerBlock + warpAlignStep - 1) / warpAlignStep) * warpAlignStep; + const auto threadsPerBlock = static_cast(rowsPerBlock * blockSize); + const auto blocks = static_cast((nPointDomain + rowsPerBlock - 1) / rowsPerBlock); InvertDiagonalBlocksKernel - <<(nPointDomain), blockSize, blockSize * sizeof(ScalarType)>>>(nPointDomain, nVar, gpu.d, - d_invM); + <<>>(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()); From dc95fd84b239f60a5931d22fd404cdec444b05f0 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 15:21:57 -0700 Subject: [PATCH 16/30] Revert "Warp-align the row batching in the Jacobi diagonal-inversion kernel" This reverts commit 0bfd05baaee1f3e56385e54ef431dee3d945f573. --- Common/src/linear_algebra/CSysMatrixGPU.cu | 56 +++++++--------------- 1 file changed, 16 insertions(+), 40 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index bcbd05d79ff..2fd81aca048 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -27,7 +27,6 @@ #include #include -#include #include "../../include/linear_algebra/CMatrixInverse.hpp" #include "../../include/linear_algebra/CSysMatrix.inl" @@ -110,39 +109,26 @@ __device__ FORCEINLINE ScalarType* GetBlockILU(const DeviceLDU& M, u /*! * \brief Invert the diagonal blocks of the matrix, device version of InverseDiagonalBlock. - * Rows (points) are batched into blocks whose thread count is a whole multiple of the - * warp size (see BuildJacobiPreconditionerGPU), unlike the previous batched attempt - * (reverted: e65bfa8eb8) which just floor-divided a ~128 thread target by nVar*nVar - for - * nVar not dividing 32 evenly (e.g. nVar=5, blockSize=25 -> 125 threads) that left the - * block's last warp partly empty. The inversion itself (SU2_LinAlg::MatrixInverse, shared - * __host__ __device__ code) still runs serially on one thread per row - it is a row- - * oriented Gauss-Jordan with a genuine sequential dependency chain across ~nVar^2/2 - * elimination substeps, so spreading it across threads would trade cheap serial FLOPs - * (~nVar^3, tiny at this size) for that many __syncthreads() barriers. Warp-alignment does - * not change that: it only avoids wasting the last warp's tail lanes on the (fully - * parallel) load into shared memory, the inversion itself still runs one point at a time - * per warp regardless of alignment. threadIdx.x maps to (row-within-block, block entry) - * via divmod by nVar*nVar. Dynamic shared memory: blockDim.x scalars, one nVar*nVar work - * buffer per row in the block. + * \note Grid: one block per row, blockDim.x == nVar*nVar (one thread per block entry, they + * stage the input in shared memory). Dynamic shared memory: nVar*nVar scalars. */ template __global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, const ScalarType* __restrict__ mat_d, ScalarType* __restrict__ invM) { - const unsigned long blockSize = nVar * nVar; - const unsigned long rowsPerBlock = blockDim.x / blockSize; - const unsigned long rowInBlock = threadIdx.x / blockSize; - const unsigned long localTid = threadIdx.x % blockSize; - const unsigned long iRow = static_cast(blockIdx.x) * rowsPerBlock + rowInBlock; - - /*--- The inversion destroys its input, so it cannot work on the matrix itself; each row in the - * batch gets its own blockSize-sized slice of shared memory. ---*/ + const unsigned long iRow = blockIdx.x; + if (iRow >= nRows) return; + + const auto blockSize = nVar * nVar; + const unsigned long tid = threadIdx.x; + + /*--- The inversion destroys its input, so it cannot work on the matrix itself. ---*/ extern __shared__ __align__(sizeof(double)) char smem[]; - auto* work = reinterpret_cast(smem) + rowInBlock * blockSize; + auto* work = reinterpret_cast(smem); - if (iRow < nRows) work[localTid] = mat_d[iRow * blockSize + localTid]; + work[tid] = mat_d[iRow * blockSize + tid]; __syncthreads(); - if (iRow < nRows && localTid == 0) SU2_LinAlg::MatrixInverse(nVar, work, invM + iRow * blockSize); + if (tid == 0) SU2_LinAlg::MatrixInverse(nVar, work, invM + iRow * blockSize); } /*! @@ -512,21 +498,11 @@ void CSysMatrix::BuildJacobiPreconditionerGPU() { 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. Round the - * number of rows batched per block up to the smallest multiple that makes blockDim.x itself a - * multiple of the warp size (32), so the block's last warp is never partly idle; then pick the - * smallest such multiple at or above the ~128 thread target instead of always taking the very - * first one, so nVar values that need a huge single step to reach alignment (e.g. nVar=5, - * blockSize=25 -> step of 32 rows/800 threads) don't blow the block far past the target. ---*/ - constexpr unsigned long targetThreadsPerBlock = 128; - const unsigned long blockSize = nVar * nVar; - const unsigned long warpAlignStep = 32 / std::gcd(32ul, blockSize); - auto rowsPerBlock = std::max(warpAlignStep, targetThreadsPerBlock / blockSize); - rowsPerBlock = ((rowsPerBlock + warpAlignStep - 1) / warpAlignStep) * warpAlignStep; - const auto threadsPerBlock = static_cast(rowsPerBlock * blockSize); - const auto blocks = static_cast((nPointDomain + rowsPerBlock - 1) / rowsPerBlock); + * CSysMatrixVectorProduct, which is created before the preconditioner is built. ---*/ + const auto blockSize = static_cast(nVar * nVar); InvertDiagonalBlocksKernel - <<>>(nPointDomain, nVar, gpu.d, d_invM); + <<(nPointDomain), blockSize, 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()); From 249209999908b6c150ba6ae115b07c28651f824d Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 15:29:08 -0700 Subject: [PATCH 17/30] Parallelize the Jacobi diagonal-inversion kernel across all block entries InvertDiagonalBlocksKernel used nVar*nVar threads per row but only thread 0 ever did real work (SU2_LinAlg::MatrixInverse, a serial forward-elimination + back-substitution chain), the rest just cooperated on the parallel load into shared memory. Replace it with full (both-directions) Gauss-Jordan elimination: for each pivot k, every entry update A(i,j) -= A(i,k)*A(k,j) for i != k is independent of every other entry in that same step, so all nVar*nVar threads do real work every step instead of just one. This drops the serial chain from one thread doing ~nVar^3 scalar ops to nVar barrier- separated steps with the full block active. Uses __syncthreads() rather than __syncwarp() so it stays correct once nVar*nVar exceeds one warp (nVar > ~5), at the cost of a full block-wide barrier even when the whole block already fits in one warp. This is a different algorithm from the shared SU2_LinAlg::MatrixInverse (host path, and IluFactorColorKernel, both left untouched), not a device branch of it, since the elimination order genuinely differs - not just a codegen difference from __CUDA_ARCH__. Verified bit-identical to the previous serial-inversion values: NACA0012 (nVar=4) GPU Jacobi/Q_Jacobi/Q_Identity reproduce their known CL/CD exactly, GPU Q_Jacobi matches the CPU reference exactly, and ONERAM6 (nVar=5) GPU Q_Jacobi matches its CPU reference to 5 decimals (same float/double gap as before). Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysMatrixGPU.cu | 56 ++++++++++++++++++---- 1 file changed, 47 insertions(+), 9 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 2fd81aca048..6a3d6a48a01 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -108,9 +108,20 @@ __device__ FORCEINLINE ScalarType* GetBlockILU(const DeviceLDU& M, u } /*! - * \brief Invert the diagonal blocks of the matrix, device version of InverseDiagonalBlock. - * \note Grid: one block per row, blockDim.x == nVar*nVar (one thread per block entry, they - * stage the input in shared memory). Dynamic shared memory: nVar*nVar scalars. + * \brief Invert the diagonal blocks of the matrix, parallel-Gauss-Jordan device version of + * InverseDiagonalBlock. Unlike SU2_LinAlg::MatrixInverse (host, and still used by + * IluFactorColorKernel below) - which does forward elimination on one thread followed by + * a serial back-substitution, an inherently sequential ~nVar^3 chain - this eliminates + * each pivot column from every OTHER row simultaneously (both above and below the pivot, + * no back-substitution needed), so every entry update within one pivot step is independent + * of every other entry in that step: nVar serial steps (each with a couple of barriers), + * instead of one thread serially doing all the work. One thread per block entry (i,j); + * __syncthreads() (not __syncwarp()) is used so this stays correct even when nVar*nVar + * exceeds one warp (nVar > ~5), at the cost of a full block-wide barrier even for the + * common case where the whole block already fits in one warp. + * \note Grid: one block per row, blockDim.x == nVar*nVar. Dynamic shared memory: 2*nVar*nVar + * scalars (the working copy of the block, and the inverse being accumulated in place of + * the old identity-then-eliminate scheme). */ template __global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, @@ -120,15 +131,43 @@ __global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nV const auto blockSize = nVar * nVar; const unsigned long tid = threadIdx.x; + const unsigned long i = tid / nVar; + const unsigned long j = tid % nVar; /*--- The inversion destroys its input, so it cannot work on the matrix itself. ---*/ extern __shared__ __align__(sizeof(double)) char smem[]; - auto* work = reinterpret_cast(smem); + auto* A = reinterpret_cast(smem); + auto* Inv = A + blockSize; - work[tid] = mat_d[iRow * blockSize + tid]; + A[tid] = mat_d[iRow * blockSize + tid]; + Inv[tid] = ScalarType(i == j); __syncthreads(); - if (tid == 0) SU2_LinAlg::MatrixInverse(nVar, work, invM + iRow * blockSize); + for (auto k = 0ul; k < nVar; ++k) { + /*--- Regularize the pivot (shared with the host path, same clamp value). ---*/ + if (i == k && j == k) SU2_LinAlg::RegularizePivot(A[k * nVar + k]); + __syncthreads(); + + /*--- Normalize the pivot row. ---*/ + const ScalarType pivot = A[k * nVar + k]; + if (i == k) { + A[tid] /= pivot; + Inv[tid] /= pivot; + } + __syncthreads(); + + /*--- Eliminate column k from every other row; A(k,*) and Inv(k,*) are already finalized for + * this step (previous barrier), and each thread only ever writes its own (i,j), so this + * needs no further synchronization until the next pivot's regularization reads A(k+1,k+1). ---*/ + if (i != k) { + const ScalarType factor = A[i * nVar + k]; + A[tid] -= factor * A[k * nVar + j]; + Inv[tid] -= factor * Inv[k * nVar + j]; + } + __syncthreads(); + } + + invM[iRow * blockSize + tid] = Inv[tid]; } /*! @@ -500,9 +539,8 @@ void CSysMatrix::BuildJacobiPreconditionerGPU() { /*--- 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, blockSize * sizeof(ScalarType)>>>(nPointDomain, nVar, gpu.d, - d_invM); + 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()); From 1e71a1b6cf99847a9072c1495f1216cd65f1108f Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 15:41:47 -0700 Subject: [PATCH 18/30] Parallelize the diagonal-block quantization kernel across rows QuantizeDiagonalBlocksKernel assigned one thread to a whole point's nVar*nVar block (looping over all nVar rows itself, each with its own max-abs scan and quantize pass), the coarsest-grained kernel in the file - every other kernel here (SpMV, Jacobi apply, diagonal inversion) assigns one thread to at most one block-row. Consecutive threads also landed nVar*nVar elements apart in mat_d, the same de-coalesced pattern ApplyJacobiPreconditionerKernel's row-major thread mapping fixed. Split EncodeQuantBlock's per-row body out into EncodeQuantRow (still SU2_CUDA_HOST_DEVICE, still shared - EncodeQuantBlock now just calls it once per row for the host path), since each block-row's scale and quantization are already fully independent of every other row. The kernel now assigns one thread per (point, row), batched into blocks of ~128 threads the same way as the Jacobi apply/SpMV kernels (threadIdx.x -> (point-within-block, row) via divmod by nVar). Verified bit-identical: NACA0012 (nVar=4) GPU Jacobi/Q_Jacobi/Q_Identity reproduce their known CL/CD exactly, GPU Q_Jacobi matches the CPU reference exactly, and ONERAM6 (nVar=5) GPU Q_Jacobi matches its CPU reference to 5 decimals (same float/double gap as before). Co-Authored-By: Claude Sonnet 5 --- Common/include/linear_algebra/CSysMatrix.hpp | 60 ++++++++++++-------- Common/src/linear_algebra/CSysMatrixGPU.cu | 32 +++++++---- 2 files changed, 57 insertions(+), 35 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 75df44b94a9..780e0817e1a 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -160,48 +160,60 @@ SU2_CUDA_HOST_DEVICE inline float DecodeQuantScale(int8_t e) noexcept { } /*! - * \brief Encode one nVar×nVar block into per-row int8 quantized storage. + * \brief Encode one row of an nVar×nVar block into int8 quantized storage: \p qs receives the + * row's scale exponent, \p qv (nVar entries) the clamped int8 values for row \p r. * \p f(r,c) is called twice per entry (max-abs scan then encoding); it should be cheap. - * Stores a per-row scale exponent in \p qs and clamped int8 values in \p qv. * \note Shared verbatim with the device (CSysMatrixGPU.cu's diagonal-quantization kernel encodes * through this same function, not a separate copy), see SU2_CUDA_HOST_DEVICE. \p f must * already return a passive (non-AD) value: every caller does, since quantization is * compile-time disabled whenever ScalarType could be AD-active (see quantized_mode), so * there is no SU2_TYPE::PassiveValue call here to keep this device-callable. The bit * reinterpretations branch on __CUDA_ARCH__ the same way DecodeQuantScale does, see there. + * Split out from EncodeQuantBlock (which just calls this once per row) so a device kernel + * can assign one thread per row instead of one thread per whole block - each row's scale + * and quantization are already fully independent of every other row. */ template -SU2_CUDA_HOST_DEVICE inline void EncodeQuantBlock(const F& f, int8_t* __restrict qs, int8_t* __restrict qv, - unsigned long nVar) noexcept { - for (auto r = 0ul; r < nVar; ++r) { - constexpr uint32_t eps_bits = 0x34000000u; - uint32_t max_abs_bits = eps_bits; - for (auto c = 0ul; c < nVar; ++c) { - const float fv = static_cast(f(r, c)); +SU2_CUDA_HOST_DEVICE inline void EncodeQuantRow(const F& f, int8_t& qs, int8_t* __restrict qv, unsigned long nVar, + unsigned long r) noexcept { + constexpr uint32_t eps_bits = 0x34000000u; + uint32_t max_abs_bits = eps_bits; + for (auto c = 0ul; c < nVar; ++c) { + const float fv = static_cast(f(r, c)); #ifdef __CUDA_ARCH__ - const uint32_t fb = __float_as_uint(fv); + const uint32_t fb = __float_as_uint(fv); #else - uint32_t fb; - memcpy(&fb, &fv, sizeof(fb)); + uint32_t fb; + memcpy(&fb, &fv, sizeof(fb)); #endif - max_abs_bits = QuantMax(max_abs_bits, fb & 0x7FFFFFFFu); - } - const int e = QuantMin(127, QuantMax(-128, static_cast(max_abs_bits >> 23) - 133)); - qs[r] = static_cast(e); - const uint32_t inv_bits = static_cast(127 - e) << 23; + max_abs_bits = QuantMax(max_abs_bits, fb & 0x7FFFFFFFu); + } + const int e = QuantMin(127, QuantMax(-128, static_cast(max_abs_bits >> 23) - 133)); + qs = static_cast(e); + const uint32_t inv_bits = static_cast(127 - e) << 23; #ifdef __CUDA_ARCH__ - const float inv_rscale = __uint_as_float(inv_bits); + const float inv_rscale = __uint_as_float(inv_bits); #else - float inv_rscale; - memcpy(&inv_rscale, &inv_bits, sizeof(inv_rscale)); + float inv_rscale; + memcpy(&inv_rscale, &inv_bits, sizeof(inv_rscale)); #endif - for (auto c = 0ul; c < nVar; ++c) { - qv[r * nVar + c] = - static_cast(QuantMax(-128.f, QuantMin(127.f, roundf(static_cast(f(r, c)) * inv_rscale)))); - } + for (auto c = 0ul; c < nVar; ++c) { + qv[c] = static_cast(QuantMax(-128.f, QuantMin(127.f, roundf(static_cast(f(r, c)) * inv_rscale)))); } } +/*! + * \brief Encode one nVar×nVar block into per-row int8 quantized storage, see EncodeQuantRow (each + * row's scale/quantization is independent, this just loops over all of them serially for + * the host path). + * \note Shared verbatim with the device, see SU2_CUDA_HOST_DEVICE and EncodeQuantRow. + */ +template +SU2_CUDA_HOST_DEVICE inline void EncodeQuantBlock(const F& f, int8_t* __restrict qs, int8_t* __restrict qv, + unsigned long nVar) noexcept { + for (auto r = 0ul; r < nVar; ++r) EncodeQuantRow(f, qs[r], qv + r * nVar, nVar, r); +} + /*! * \brief View of one matrix block, const-correct via the ScalarType template parameter. * \c CBlockView is read-only; \c CBlockView is mutable diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 6a3d6a48a01..8945db07dfb 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -172,21 +172,29 @@ __global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nV /*! * \brief Quantize the diagonal blocks straight from the device diagonal (gpu.d), device - * counterpart of CSysMatrix::QuantizeBlock (CSysMatrix.cpp), applied row by row. Calls the - * exact same EncodeQuantBlock encoding routine (CSysMatrix.hpp, SU2_CUDA_HOST_DEVICE) that - * the host path uses, rather than a separate device copy of the encoding logic. One - * thread per row. + * counterpart of CSysMatrix::QuantizeBlock (CSysMatrix.cpp). Calls the exact same + * EncodeQuantRow encoding routine (CSysMatrix.hpp, SU2_CUDA_HOST_DEVICE) that the host + * path's EncodeQuantBlock loops over, rather than a separate device copy of the encoding + * logic. One thread per (point, row) - each block-row's scale and quantization are + * already independent of every other row (see EncodeQuantRow), so this is nVar times + * more parallel than one-thread-per-point, and consecutive threads land nVar elements + * apart in mat_d instead of nVar*nVar apart, a partial coalescing win the same way + * ApplyJacobiPreconditionerKernel's row-major thread mapping is. Points are batched into + * blocks of ~128 threads the same way (threadIdx.x -> (point-within-block, row) via + * divmod by nVar). */ template __global__ void QuantizeDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, const ScalarType* __restrict__ mat_d, int8_t* __restrict__ q_scale_d, int8_t* __restrict__ q_blocks_d) { - const auto iRow = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; - if (iRow >= nRows) return; + const unsigned long rowsPerBlock = blockDim.x / nVar; + const unsigned long iVar = threadIdx.x % nVar; + const unsigned long iPoint = static_cast(blockIdx.x) * rowsPerBlock + threadIdx.x / nVar; + if (iPoint >= nRows) return; - const auto* blk = mat_d + iRow * nVar * nVar; - EncodeQuantBlock([&](unsigned long r, unsigned long c) { return blk[r * nVar + c]; }, q_scale_d + iRow * nVar, - q_blocks_d + iRow * nVar * nVar, nVar); + const auto* blk = mat_d + iPoint * nVar * nVar; + EncodeQuantRow([&](unsigned long r, unsigned long c) { return blk[r * nVar + c]; }, q_scale_d[iPoint * nVar + iVar], + q_blocks_d + iPoint * nVar * nVar + iVar * nVar, nVar, iVar); } /*! @@ -518,8 +526,10 @@ void CSysMatrix::QuantizeDiagonalBlocksGPU() { /*--- 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. ---*/ - constexpr unsigned threadsPerBlock = 128; - const auto blocks = static_cast((nPointDomain + threadsPerBlock - 1) / threadsPerBlock); + constexpr unsigned long targetThreadsPerBlock = 128; + const auto rowsPerBlock = std::max(1, targetThreadsPerBlock / nVar); + const auto threadsPerBlock = static_cast(rowsPerBlock * nVar); + const auto blocks = static_cast((nPointDomain + rowsPerBlock - 1) / rowsPerBlock); QuantizeDiagonalBlocksKernel <<>>(nPointDomain, nVar, gpu.d, d_q_scale.d, d_q_blocks.d); /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ From ccd25640121e1c1d43357c0d7647bdd184edcb7a Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 15:50:39 -0700 Subject: [PATCH 19/30] Overlap the quantized L/U H2D transfer with default-stream kernels HtDTransfer issued the async L/U copies on the default stream, the same stream QuantizeDiagonalBlocksGPU (and everything else) launches on. A single stream is strictly ordered, so despite being "async" copies, the diagonal-quantization kernel just queued behind them instead of actually running concurrently on the GPU - the copy engine and the SMs are independent hardware and there is no real data dependency between the two (the kernel only touches gpu.d and d_q_scale/blocks.d, the transfer only touches d_q_scale/blocks.l/u). Move the L/U copies to a dedicated stream (htd_stream, same pattern as the existing ilu_stream) and record an event (htd_event) at the end of them. The one default-stream consumer that actually reads the result - the quantized SpMV kernel in MatrixVectorProductGPU - waits on that event via cudaStreamWaitEvent before launching, a device-side wait that does not block the host thread, unlike the cudaStreamSynchronize calls elsewhere in this file. Verified bit-identical over 5 repeated runs (to catch any nondeterministic cross-stream race): NACA0012 (nVar=4) GPU Jacobi/Q_Jacobi/Q_Identity reproduce their known CL/CD exactly, GPU Q_Jacobi matches the CPU reference exactly, and ONERAM6 (nVar=5, several outer iterations so HtDTransfer/the event fire repeatedly) GPU Q_Jacobi matches its CPU reference to 5 decimals. Co-Authored-By: Claude Sonnet 5 --- Common/include/linear_algebra/CSysMatrix.hpp | 8 ++++ Common/src/linear_algebra/CSysMatrix.cpp | 2 + Common/src/linear_algebra/CSysMatrixGPU.cu | 40 +++++++++++++------- 3 files changed, 37 insertions(+), 13 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 780e0817e1a..2656df94664 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -414,6 +414,14 @@ class CSysMatrix { /*--- The legacy default stream cannot be captured into a graph. ---*/ mutable struct CUstream_st* ilu_stream = nullptr; + /*--- Dedicated stream for the async H2D transfer of the quantized L/U blocks (HtDTransfer), + * so that transfer can run concurrently (copy engine) with kernels issued on the default + * stream (e.g. QuantizeDiagonalBlocksGPU, on the SM), instead of queueing behind them on the + * same stream. htd_event marks the end of that transfer so the default-stream kernel that + * first reads the result (the quantized SpMV) can wait on it without a host-side block. ---*/ + mutable struct CUstream_st* htd_stream = nullptr; + mutable struct CUevent_st* htd_event = nullptr; + ScalarType* invM; /*!< \brief Inverse of (Jacobi) preconditioner. */ /*--- Temporary (hence mutable) working memory used in the Linelet preconditioner, outer vector is for threads ---*/ diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 7b3e87ebdaf..ef49bde5072 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -164,6 +164,8 @@ CSysMatrix::~CSysMatrix() { if (ilu_build_graph_exec != nullptr) cudaGraphExecDestroy(ilu_build_graph_exec); if (ilu_apply_graph_exec != nullptr) cudaGraphExecDestroy(ilu_apply_graph_exec); if (ilu_stream != nullptr) cudaStreamDestroy(ilu_stream); + if (htd_stream != nullptr) cudaStreamDestroy(htd_stream); + if (htd_event != nullptr) cudaEventDestroy(htd_event); #endif } diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 8945db07dfb..2e4d12d94d4 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -698,21 +698,29 @@ void CSysMatrix::HtDTransfer(bool trigger) const { /*--- No gpu.l/gpu.u to transfer (never allocated); mirror the host quantized off-diagonal * storage instead (the diagonal mirrors, d_q_scale.d/d_q_blocks.d, are not touched here at * all - QuantizeDiagonalBlocksGPU() populates them straight from gpu.d, just uploaded above, - * with no host round trip). Issued as async copies on the default stream, then left in - * flight: the caller (CSysMatrixVectorProduct's constructor) returns right after this, and - * whatever the preconditioner's Build() does next runs concurrently with the transfer still - * draining. This is only genuinely asynchronous (i.e. the host thread does not block here - * waiting for the copy) because q_scale.l/q_blocks.l/q_scale.u/q_blocks.u are pinned host - * memory, see the comment on those members / Initialize() - cudaMemcpyAsync silently - * degrades to a blocking copy from regular pageable memory. Any later kernel that reads - * d_q_scale.l/.u/d_q_blocks.l/.u is issued on the same default stream too, so it correctly - * waits for these by stream ordering without an explicit sync here. ---*/ - gpuErrChk(cudaMemcpyAsync(d_q_scale.l, q_scale.l, sizeof(QuantType) * mat.nnz_l * nVar, cudaMemcpyHostToDevice)); + * with no host round trip). Issued as async copies on a dedicated stream (not the default + * one), so this transfer (the copy engine) can actually run concurrently with whatever the + * preconditioner's Build() launches next on the default stream (the SMs), e.g. + * QuantizeDiagonalBlocksGPU - on the default stream they would just queue behind these + * copies instead of overlapping, since a single stream is strictly ordered. This is only + * genuinely asynchronous (i.e. the host thread does not block here waiting for the copy) + * because q_scale.l/q_blocks.l/q_scale.u/q_blocks.u are pinned host memory, see the comment + * on those members / Initialize() - cudaMemcpyAsync silently degrades to a blocking copy + * from regular pageable memory. htd_event marks the end of this transfer: the first + * default-stream kernel to actually read d_q_scale.l/.u/d_q_blocks.l/.u (the quantized SpMV, + * MatrixVectorProductGPU) waits on it there, since cross-stream dependencies are not + * implied by stream ordering the way same-stream ones are. ---*/ + if (htd_stream == nullptr) gpuErrChk(cudaStreamCreate(&htd_stream)); + if (htd_event == nullptr) gpuErrChk(cudaEventCreateWithFlags(&htd_event, cudaEventDisableTiming)); + gpuErrChk(cudaMemcpyAsync(d_q_scale.l, q_scale.l, sizeof(QuantType) * mat.nnz_l * nVar, cudaMemcpyHostToDevice, + htd_stream)); gpuErrChk(cudaMemcpyAsync(d_q_blocks.l, q_blocks.l, sizeof(QuantType) * mat.nnz_l * nVar * nEqn, - cudaMemcpyHostToDevice)); - gpuErrChk(cudaMemcpyAsync(d_q_scale.u, q_scale.u, sizeof(QuantType) * mat.nnz_u * nVar, cudaMemcpyHostToDevice)); + cudaMemcpyHostToDevice, htd_stream)); + gpuErrChk(cudaMemcpyAsync(d_q_scale.u, q_scale.u, sizeof(QuantType) * mat.nnz_u * nVar, cudaMemcpyHostToDevice, + htd_stream)); gpuErrChk(cudaMemcpyAsync(d_q_blocks.u, q_blocks.u, sizeof(QuantType) * mat.nnz_u * nVar * nEqn, - cudaMemcpyHostToDevice)); + cudaMemcpyHostToDevice, htd_stream)); + gpuErrChk(cudaEventRecord(htd_event, htd_stream)); } else { gpuErrChk(cudaMemcpy(gpu.l, mat.l, sizeof(ScalarType) * mat.nnz_l * nVar * nEqn, cudaMemcpyHostToDevice)); gpuErrChk(cudaMemcpy(gpu.u, mat.u, sizeof(ScalarType) * mat.nnz_u * nVar * nEqn, cudaMemcpyHostToDevice)); @@ -735,6 +743,12 @@ void CSysMatrix::MatrixVectorProductGPU(const CSysVector dim3 blockDim(static_cast(rowsPerBlock * nVar), 1, 1); dim3 gridDim(static_cast((nPointDomain + rowsPerBlock - 1) / rowsPerBlock), 1, 1); if (quantized_mode) { + /*--- Wait (on the device, no host block) for HtDTransfer's async L/U copy on its own stream + * to finish before this default-stream kernel reads d_q_scale.l/.u/d_q_blocks.l/.u - a + * cross-stream dependency, so it is not implied by ordering the way same-stream launches are. + * htd_event is always valid here: HtDTransfer runs once per solve before the first call to + * this function (see its comment), so it has already recorded the event at least once. ---*/ + gpuErrChk(cudaStreamWaitEvent(nullptr, htd_event, 0)); QuantizedBlockLDU_SpMV_kernel<<>>( nPointDomain, nVar, gpu.row_ptr_l, gpu.col_ind_l, d_q_scale.l, d_q_blocks.l, d_q_scale.d, d_q_blocks.d, gpu.row_ptr_u, gpu.col_ind_u, d_q_scale.u, d_q_blocks.u, d_vec, d_prod); From 36e98c47a3a636e4d935fe41190df125658a602b Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 16:05:12 -0700 Subject: [PATCH 20/30] Reuse the parallel Gauss-Jordan inversion for ILU's diagonal blocks IluFactorColorKernel's diagonal-inversion tail had the exact same shape InvertDiagonalBlocksKernel used to: nVar*nVar threads already active in the block, but only tid==0 doing the actual SU2_LinAlg::MatrixInverse work while the rest sat idle. It also already had two blockSize-sized shared buffers on hand (work, and Lij - free by that point, its last use storing it into Block_ij is already done), the exact inputs ParallelMatrixInverse needs. Factor the parallel algorithm out of InvertDiagonalBlocksKernel into a standalone ParallelMatrixInverse device function and call it from both kernels, instead of duplicating the elimination loop. Verified bit-identical (including against a pre-change build of the same config) over 5 repeated runs of GPU ILU (to catch any nondeterministic issue from reusing Lij's shared memory as the second buffer): NACA0012 (nVar=4) GPU Jacobi/Q_Jacobi/Q_Identity/ILU all reproduce their known/previous values exactly, GPU Q_Jacobi matches the CPU reference exactly. Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysMatrixGPU.cu | 96 ++++++++++++++-------- 1 file changed, 60 insertions(+), 36 deletions(-) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 2e4d12d94d4..243e4faa236 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -108,17 +108,57 @@ __device__ FORCEINLINE ScalarType* GetBlockILU(const DeviceLDU& M, u } /*! - * \brief Invert the diagonal blocks of the matrix, parallel-Gauss-Jordan device version of - * InverseDiagonalBlock. Unlike SU2_LinAlg::MatrixInverse (host, and still used by - * IluFactorColorKernel below) - which does forward elimination on one thread followed by - * a serial back-substitution, an inherently sequential ~nVar^3 chain - this eliminates - * each pivot column from every OTHER row simultaneously (both above and below the pivot, - * no back-substitution needed), so every entry update within one pivot step is independent - * of every other entry in that step: nVar serial steps (each with a couple of barriers), - * instead of one thread serially doing all the work. One thread per block entry (i,j); - * __syncthreads() (not __syncwarp()) is used so this stays correct even when nVar*nVar - * exceeds one warp (nVar > ~5), at the cost of a full block-wide barrier even for the - * common case where the whole block already fits in one warp. + * \brief Parallel Gauss-Jordan matrix inversion, shared by InvertDiagonalBlocksKernel and + * IluFactorColorKernel's diagonal-inversion step below - both already have nVar*nVar + * threads and two blockSize-sized shared buffers on hand at the point they need a diagonal + * block inverted. Unlike SU2_LinAlg::MatrixInverse (host, and the one still used to invert + * Linelet's tridiagonal blocks) - which does forward elimination on one thread followed by + * a serial back-substitution, an inherently sequential ~nVar^3 chain - this eliminates each + * pivot column from every OTHER row simultaneously (both above and below the pivot, no + * back-substitution needed), so every entry update within one pivot step is independent of + * every other entry in that step: nVar serial steps (each with a couple of barriers), + * instead of one thread serially doing all the work. + * \param i,j Row/column of the block entry this thread owns, in 0..nVar-1 (i.e. threadIdx.x's + * divmod by nVar, same mapping the caller already uses for everything else). + * \param A Destroyed. \param Inv Must be pre-loaded with the identity, must not alias \p A; + * holds the inverse on return. + * \note __syncthreads() (not __syncwarp()) is used so this stays correct even when nVar*nVar + * exceeds one warp (nVar > ~5), at the cost of a full block-wide barrier even for the + * common case where the whole block already fits in one warp. Every thread of the block + * must call this (no divergent early return before it), since every __syncthreads() here + * is a whole-block barrier. + */ +template +__device__ FORCEINLINE void ParallelMatrixInverse(unsigned long nVar, unsigned long i, unsigned long j, + ScalarType* __restrict__ A, ScalarType* __restrict__ Inv) { + for (auto k = 0ul; k < nVar; ++k) { + /*--- Regularize the pivot (shared with the host path, same clamp value). ---*/ + if (i == k && j == k) SU2_LinAlg::RegularizePivot(A[k * nVar + k]); + __syncthreads(); + + /*--- Normalize the pivot row. ---*/ + const ScalarType pivot = A[k * nVar + k]; + if (i == k) { + A[i * nVar + j] /= pivot; + Inv[i * nVar + j] /= pivot; + } + __syncthreads(); + + /*--- Eliminate column k from every other row; A(k,*) and Inv(k,*) are already finalized for + * this step (previous barrier), and each thread only ever writes its own (i,j), so this + * needs no further synchronization until the next pivot's regularization reads A(k+1,k+1). ---*/ + if (i != k) { + const ScalarType factor = A[i * nVar + k]; + A[i * nVar + j] -= factor * A[k * nVar + j]; + Inv[i * nVar + j] -= factor * Inv[k * nVar + j]; + } + __syncthreads(); + } +} + +/*! + * \brief Invert the diagonal blocks of the matrix, device version of InverseDiagonalBlock, via + * ParallelMatrixInverse (see its comment for why this beats a single serial thread). * \note Grid: one block per row, blockDim.x == nVar*nVar. Dynamic shared memory: 2*nVar*nVar * scalars (the working copy of the block, and the inverse being accumulated in place of * the old identity-then-eliminate scheme). @@ -143,29 +183,7 @@ __global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nV Inv[tid] = ScalarType(i == j); __syncthreads(); - for (auto k = 0ul; k < nVar; ++k) { - /*--- Regularize the pivot (shared with the host path, same clamp value). ---*/ - if (i == k && j == k) SU2_LinAlg::RegularizePivot(A[k * nVar + k]); - __syncthreads(); - - /*--- Normalize the pivot row. ---*/ - const ScalarType pivot = A[k * nVar + k]; - if (i == k) { - A[tid] /= pivot; - Inv[tid] /= pivot; - } - __syncthreads(); - - /*--- Eliminate column k from every other row; A(k,*) and Inv(k,*) are already finalized for - * this step (previous barrier), and each thread only ever writes its own (i,j), so this - * needs no further synchronization until the next pivot's regularization reads A(k+1,k+1). ---*/ - if (i != k) { - const ScalarType factor = A[i * nVar + k]; - A[tid] -= factor * A[k * nVar + j]; - Inv[tid] -= factor * Inv[k * nVar + j]; - } - __syncthreads(); - } + ParallelMatrixInverse(nVar, i, j, A, Inv); invM[iRow * blockSize + tid] = Inv[tid]; } @@ -290,11 +308,17 @@ __global__ void IluFactorColorKernel(const su2uint* __restrict__ color_idx, unsi } /*--- Invert the diagonal entry, Uii, for the rows that depend on it. The loop above may have - * updated it (when kPoint == iRow), so the whole block has to be done first. ---*/ + * updated it (when kPoint == iRow), so the whole block has to be done first. Lij is free again + * here (its last use, storing it into Block_ij, is done) - reuse it as the identity/inverse + * buffer ParallelMatrixInverse needs, instead of a separate shared allocation. ---*/ __syncthreads(); work[tid] = M.d[iRow * blockSize + tid]; + Lij[tid] = ScalarType(iVar == jVar); __syncthreads(); - if (tid == 0) SU2_LinAlg::MatrixInverse(nVar, work, M.d + iRow * blockSize); + + ParallelMatrixInverse(nVar, iVar, jVar, work, Lij); + + M.d[iRow * blockSize + tid] = Lij[tid]; } /*! From 54a93a4e79f01170a5d555523f8a896e80ad4568 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 19:00:42 -0700 Subject: [PATCH 21/30] Share one non-default stream between ILU and the quantized L/U transfer ilu_stream and htd_stream never actually need to exist at the same time: quantized_mode (which needs htd_stream for HtDTransfer's async L/U copy) and ILU (which needs ilu_stream for its build/apply CUDA graphs) are alternative preconditioner choices decided once in Initialize() - a given CSysMatrix instance is never both. Merge them into a single aux_stream member instead of keeping two separate streams (and two create/destroy pairs) that are never both live. Verified: GPU Q_Jacobi matches CPU Q_Jacobi exactly and GPU ILU converges normally, both reproducible over 3 repeated runs. Co-Authored-By: Claude Sonnet 5 --- Common/include/linear_algebra/CSysMatrix.hpp | 20 ++++--- Common/src/linear_algebra/CSysMatrix.cpp | 3 +- Common/src/linear_algebra/CSysMatrixGPU.cu | 60 ++++++++++---------- 3 files changed, 43 insertions(+), 40 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 2656df94664..4d53a39e5ae 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -411,15 +411,17 @@ class CSysMatrix { * was captured with, to detect when * it must be recaptured. */ mutable ScalarType* ilu_apply_graph_prod = nullptr; - /*--- The legacy default stream cannot be captured into a graph. ---*/ - mutable struct CUstream_st* ilu_stream = nullptr; - - /*--- Dedicated stream for the async H2D transfer of the quantized L/U blocks (HtDTransfer), - * so that transfer can run concurrently (copy engine) with kernels issued on the default - * stream (e.g. QuantizeDiagonalBlocksGPU, on the SM), instead of queueing behind them on the - * same stream. htd_event marks the end of that transfer so the default-stream kernel that - * first reads the result (the quantized SpMV) can wait on it without a host-side block. ---*/ - mutable struct CUstream_st* htd_stream = 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 + * cannot be captured into a graph; (2) HtDTransfer's async H2D transfer of the quantized L/U + * blocks, so that transfer can run concurrently (copy engine) with kernels issued on the + * default stream (e.g. QuantizeDiagonalBlocksGPU, on the SM) instead of queueing behind them on + * the same stream. Because the two uses are mutually exclusive, sharing one stream (rather than + * a dedicated one per use) needs no extra synchronization between them. htd_event marks the end + * of the H2D transfer specifically, so the default-stream kernel that first reads the result + * (the quantized SpMV) can wait on it without a host-side block. ---*/ + mutable struct CUstream_st* aux_stream = nullptr; mutable struct CUevent_st* htd_event = nullptr; ScalarType* invM; /*!< \brief Inverse of (Jacobi) preconditioner. */ diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index ef49bde5072..3a5e6041208 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -163,8 +163,7 @@ CSysMatrix::~CSysMatrix() { #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 (ilu_stream != nullptr) cudaStreamDestroy(ilu_stream); - if (htd_stream != nullptr) cudaStreamDestroy(htd_stream); + if (aux_stream != nullptr) cudaStreamDestroy(aux_stream); if (htd_event != nullptr) cudaEventDestroy(htd_event); #endif } diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 243e4faa236..78d39a36a32 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -602,7 +602,7 @@ void CSysMatrix::BuildILUPreconditionerGPU() { /*--- The legacy default stream cannot be captured, so the graph lives on its own stream, * created once. Every launch below is followed by a sync back to the host, so this does not * change execution order relative to the rest of the (single-stream) solver. ---*/ - if (ilu_stream == nullptr) gpuErrChk(cudaStreamCreate(&ilu_stream)); + if (aux_stream == nullptr) gpuErrChk(cudaStreamCreate(&aux_stream)); /*--- The launch sequence (ilu_gpu_sweeps passes over all colors) is identical on every call: * the grid and block sizes only depend on the (fixed) sparsity pattern/coloring and the device @@ -616,7 +616,7 @@ void CSysMatrix::BuildILUPreconditionerGPU() { * scratch, relying on the matrix changing little between outer/pseudo-time iterations. ---*/ if (ilu_build_graph_exec == nullptr) { cudaGraph_t graph; - gpuErrChk(cudaStreamBeginCapture(ilu_stream, cudaStreamCaptureModeThreadLocal)); + gpuErrChk(cudaStreamBeginCapture(aux_stream, cudaStreamCaptureModeThreadLocal)); for (unsigned short sweep = 0; sweep < ilu_gpu_sweeps; ++sweep) { for (auto color = 0ul; color + 1 < ilu_color_ptr.size(); ++color) { @@ -624,17 +624,17 @@ void CSysMatrix::BuildILUPreconditionerGPU() { const auto size = ilu_color_ptr[color + 1] - begin; if (size == 0) continue; IluFactorColorKernel - <<>>(d_ilu_color_idx, begin, size, nPointDomain, nVar, A, M); + <<>>(d_ilu_color_idx, begin, size, nPointDomain, nVar, A, M); } } - gpuErrChk(cudaStreamEndCapture(ilu_stream, &graph)); + gpuErrChk(cudaStreamEndCapture(aux_stream, &graph)); gpuErrChk(cudaGraphInstantiate(&ilu_build_graph_exec, graph, nullptr, nullptr, 0)); gpuErrChk(cudaGraphDestroy(graph)); } - gpuErrChk(cudaGraphLaunch(ilu_build_graph_exec, ilu_stream)); - gpuErrChk(cudaStreamSynchronize(ilu_stream)); + gpuErrChk(cudaGraphLaunch(ilu_build_graph_exec, aux_stream)); + gpuErrChk(cudaStreamSynchronize(aux_stream)); gpuErrChk(cudaGetLastError()); } @@ -664,7 +664,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_ilu_level_idx, begin, size, nVar, M, d_vec, d_prod); } /*--- Backward substitution: one exact pass over the levels in decreasing order, @@ -698,18 +698,18 @@ void CSysMatrix::ComputeILUPreconditionerGPU(const CSysVector - <<>>(d_ilu_level_idx, begin, size, nPointDomain, nVar, M, d_prod); + <<>>(d_ilu_level_idx, begin, size, nPointDomain, nVar, M, d_prod); } - gpuErrChk(cudaStreamEndCapture(ilu_stream, &graph)); + gpuErrChk(cudaStreamEndCapture(aux_stream, &graph)); gpuErrChk(cudaGraphInstantiate(&ilu_apply_graph_exec, graph, nullptr, nullptr, 0)); gpuErrChk(cudaGraphDestroy(graph)); ilu_apply_graph_vec = d_vec; ilu_apply_graph_prod = d_prod; } - gpuErrChk(cudaGraphLaunch(ilu_apply_graph_exec, ilu_stream)); - gpuErrChk(cudaStreamSynchronize(ilu_stream)); + gpuErrChk(cudaGraphLaunch(ilu_apply_graph_exec, aux_stream)); + gpuErrChk(cudaStreamSynchronize(aux_stream)); gpuErrChk(cudaGetLastError()); } @@ -722,29 +722,31 @@ void CSysMatrix::HtDTransfer(bool trigger) const { /*--- No gpu.l/gpu.u to transfer (never allocated); mirror the host quantized off-diagonal * storage instead (the diagonal mirrors, d_q_scale.d/d_q_blocks.d, are not touched here at * all - QuantizeDiagonalBlocksGPU() populates them straight from gpu.d, just uploaded above, - * with no host round trip). Issued as async copies on a dedicated stream (not the default - * one), so this transfer (the copy engine) can actually run concurrently with whatever the + * with no host round trip). Issued as async copies on aux_stream (not the default one), so + * this transfer (the copy engine) can actually run concurrently with whatever the * preconditioner's Build() launches next on the default stream (the SMs), e.g. * QuantizeDiagonalBlocksGPU - on the default stream they would just queue behind these - * copies instead of overlapping, since a single stream is strictly ordered. This is only - * genuinely asynchronous (i.e. the host thread does not block here waiting for the copy) - * because q_scale.l/q_blocks.l/q_scale.u/q_blocks.u are pinned host memory, see the comment - * on those members / Initialize() - cudaMemcpyAsync silently degrades to a blocking copy - * from regular pageable memory. htd_event marks the end of this transfer: the first - * default-stream kernel to actually read d_q_scale.l/.u/d_q_blocks.l/.u (the quantized SpMV, - * MatrixVectorProductGPU) waits on it there, since cross-stream dependencies are not - * implied by stream ordering the way same-stream ones are. ---*/ - if (htd_stream == nullptr) gpuErrChk(cudaStreamCreate(&htd_stream)); + * copies instead of overlapping, since a single stream is strictly ordered. aux_stream is + * shared with the (mutually exclusive, see its declaration) ILU build/apply graphs rather + * than using a separate dedicated stream. This is only genuinely asynchronous (i.e. the host + * thread does not block here waiting for the copy) because q_scale.l/q_blocks.l/q_scale.u/ + * q_blocks.u are pinned host memory, see the comment on those members / Initialize() - + * cudaMemcpyAsync silently degrades to a blocking copy from regular pageable memory. + * htd_event marks the end of this transfer: the first default-stream kernel to actually read + * d_q_scale.l/.u/d_q_blocks.l/.u (the quantized SpMV, MatrixVectorProductGPU) waits on it + * there, since cross-stream dependencies are not implied by stream ordering the way + * same-stream ones are. ---*/ + if (aux_stream == nullptr) gpuErrChk(cudaStreamCreate(&aux_stream)); if (htd_event == nullptr) gpuErrChk(cudaEventCreateWithFlags(&htd_event, cudaEventDisableTiming)); gpuErrChk(cudaMemcpyAsync(d_q_scale.l, q_scale.l, sizeof(QuantType) * mat.nnz_l * nVar, cudaMemcpyHostToDevice, - htd_stream)); + aux_stream)); gpuErrChk(cudaMemcpyAsync(d_q_blocks.l, q_blocks.l, sizeof(QuantType) * mat.nnz_l * nVar * nEqn, - cudaMemcpyHostToDevice, htd_stream)); + cudaMemcpyHostToDevice, aux_stream)); gpuErrChk(cudaMemcpyAsync(d_q_scale.u, q_scale.u, sizeof(QuantType) * mat.nnz_u * nVar, cudaMemcpyHostToDevice, - htd_stream)); + aux_stream)); gpuErrChk(cudaMemcpyAsync(d_q_blocks.u, q_blocks.u, sizeof(QuantType) * mat.nnz_u * nVar * nEqn, - cudaMemcpyHostToDevice, htd_stream)); - gpuErrChk(cudaEventRecord(htd_event, htd_stream)); + cudaMemcpyHostToDevice, aux_stream)); + gpuErrChk(cudaEventRecord(htd_event, aux_stream)); } else { gpuErrChk(cudaMemcpy(gpu.l, mat.l, sizeof(ScalarType) * mat.nnz_l * nVar * nEqn, cudaMemcpyHostToDevice)); gpuErrChk(cudaMemcpy(gpu.u, mat.u, sizeof(ScalarType) * mat.nnz_u * nVar * nEqn, cudaMemcpyHostToDevice)); From 957e8d097b173a9b39fe2cbf79c9b5d0fd82090c Mon Sep 17 00:00:00 2001 From: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> Date: Sat, 15 Aug 2026 19:21:52 -0700 Subject: [PATCH 22/30] Apply suggestions from code review Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- Common/include/code_config.hpp | 5 +- .../linear_algebra/CPreconditioner.hpp | 11 ++--- Common/include/linear_algebra/CSysMatrix.hpp | 45 +++-------------- Common/src/linear_algebra/CSysMatrix.cpp | 19 ++----- Common/src/linear_algebra/CSysMatrixGPU.cu | 49 ++++--------------- 5 files changed, 25 insertions(+), 104 deletions(-) diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp index 93c82eaad88..a28c346b651 100644 --- a/Common/include/code_config.hpp +++ b/Common/include/code_config.hpp @@ -50,10 +50,7 @@ #define NEVERINLINE inline #endif -/*--- Marks a function callable from both host and device code (e.g. a small numerical kernel - * shared between a CSysMatrix host routine and its CSysMatrixGPU.cu counterpart), a no-op outside - * nvcc. Centralized here (rather than defined locally, as e.g. CMatrixInverse.hpp used to) so - * every shared routine uses the exact same macro instead of independent copies drifting apart. ---*/ +/*--- Marks a function callable from both host and device code, a no-op outside nvcc. ---*/ #ifdef __CUDACC__ #define SU2_CUDA_HOST_DEVICE __host__ __device__ #else diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index 98e252749ff..2b6c67701ff 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -108,10 +108,8 @@ CPreconditioner::~CPreconditioner() {} * \class CIdentityPreconditioner * \brief No-op preconditioner used when Krylov solvers run without preconditioning. * \note Also serves Q_IDENTITY: Build() requests quantization of the diagonal blocks, needed by - * the matrix-vector product shared with the Krylov solver (off diagonals are quantized on the - * fly during assembly), even though this preconditioner's own operation is a no-op either way; - * CSysMatrix::QuantizeDiagonalBlocks() itself no-ops unless the matrix was actually set up for - * quantization (Q_IDENTITY/Q_JACOBI/Q_LU_SGS), so this is free for plain IDENTITY. + * the matrix-vector product shared with the Krylov solver even though this preconditioner's own + * operation is a no-op. CSysMatrix::QuantizeDiagonalBlocks() when quantization is off. */ template class CIdentityPreconditioner final : public CPreconditioner { @@ -173,7 +171,7 @@ class CJacobiPreconditioner final : public CPreconditioner { /*! * \note Request the associated matrix to build the preconditioner. Also serves Q_JACOBI: * BuildJacobiPreconditioner() quantizes the diagonal blocks itself when the matrix was - * set up for it, so there is nothing extra to do here for the quantized case. + * set up for it. */ inline void Build() override { sparse_matrix.BuildJacobiPreconditioner(); } }; @@ -266,8 +264,7 @@ class CLU_SGSPreconditioner final : public CPreconditioner { } /*! - * \note Also serves Q_LU_SGS: quantizes the diagonal blocks (off diagonals are quantized on - * the fly during assembly); a no-op for plain LU_SGS. + * \note Also serves Q_LU_SGS: quantizes the diagonal blocks, no-op for plain LU_SGS. */ inline void Build() override { sparse_matrix.QuantizeDiagonalBlocks(); } }; diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 2656df94664..3f25fe2908f 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -113,11 +113,7 @@ struct CSysMatrixComms { /*! * \brief std::max/std::min, usable from both host and device code. Calling std::max/std::min * directly from a SU2_CUDA_HOST_DEVICE function compiles without error but is not - * actually valid without --expt-relaxed-constexpr (not used in this build) - nvcc only - * warns ("calling a constexpr __host__ function ... is not allowed"), then silently - * emits device code that does not do what it looks like it does. CUDA's device compiler - * provides its own (unqualified) max/min built-ins instead, which these two forward to; - * on the host they just forward to std::max/std::min. + * actually valid without --expt-relaxed-constexpr which is not used in this build. */ template SU2_CUDA_HOST_DEVICE inline T QuantMax(T a, T b) noexcept { @@ -141,12 +137,7 @@ SU2_CUDA_HOST_DEVICE inline T QuantMin(T a, T b) noexcept { * The exponent \p e was packed as (e + 127) into the IEEE 754 biased-exponent field * with a zero mantissa, giving an exact power of two: 2^e. * This is the inverse of the encoding in EncodeQuantBlock. - * \note Shared verbatim with the device (the quantized SpMV kernel in CSysMatrixGPU.cu decodes - * through this same function, not a separate copy), see SU2_CUDA_HOST_DEVICE. The bit - * reinterpretation itself branches on __CUDA_ARCH__ (mirroring RegularizePivot in - * CMatrixInverse.hpp): __uint_as_float on device, memcpy on host - plain memcpy compiles - * for the device too, but was observed to silently misinterpret the bits there rather than - * reinterpreting them. + * \note Branches on __CUDA_ARCH__, plain memcpy compiles for the device but does not work! */ SU2_CUDA_HOST_DEVICE inline float DecodeQuantScale(int8_t e) noexcept { const uint32_t bits = static_cast(QuantMax(0, static_cast(e) + 127)) << 23; @@ -163,15 +154,7 @@ SU2_CUDA_HOST_DEVICE inline float DecodeQuantScale(int8_t e) noexcept { * \brief Encode one row of an nVar×nVar block into int8 quantized storage: \p qs receives the * row's scale exponent, \p qv (nVar entries) the clamped int8 values for row \p r. * \p f(r,c) is called twice per entry (max-abs scan then encoding); it should be cheap. - * \note Shared verbatim with the device (CSysMatrixGPU.cu's diagonal-quantization kernel encodes - * through this same function, not a separate copy), see SU2_CUDA_HOST_DEVICE. \p f must - * already return a passive (non-AD) value: every caller does, since quantization is - * compile-time disabled whenever ScalarType could be AD-active (see quantized_mode), so - * there is no SU2_TYPE::PassiveValue call here to keep this device-callable. The bit - * reinterpretations branch on __CUDA_ARCH__ the same way DecodeQuantScale does, see there. - * Split out from EncodeQuantBlock (which just calls this once per row) so a device kernel - * can assign one thread per row instead of one thread per whole block - each row's scale - * and quantization are already fully independent of every other row. + * \note Shared with the device and thus same __CUDA_ARCH__ branches as DecodeQuantScale. */ template SU2_CUDA_HOST_DEVICE inline void EncodeQuantRow(const F& f, int8_t& qs, int8_t* __restrict qv, unsigned long nVar, @@ -206,7 +189,6 @@ SU2_CUDA_HOST_DEVICE inline void EncodeQuantRow(const F& f, int8_t& qs, int8_t* * \brief Encode one nVar×nVar block into per-row int8 quantized storage, see EncodeQuantRow (each * row's scale/quantization is independent, this just loops over all of them serially for * the host path). - * \note Shared verbatim with the device, see SU2_CUDA_HOST_DEVICE and EncodeQuantRow. */ template SU2_CUDA_HOST_DEVICE inline void EncodeQuantBlock(const F& f, int8_t* __restrict qs, int8_t* __restrict qv, @@ -337,25 +319,13 @@ class CSysMatrix { #endif /*!< \brief Per-row exponents; .l/.u sized [nnz_l/u * nVar], .d [nPoint * nVar]. .l/.u are * populated during assembly (quantized on the fly); .d is populated by - * QuantizeDiagonalBlocks() on the host, but only when !useCuda - under CUDA the - * diagonal is quantized straight from gpu.d instead (QuantizeDiagonalBlocksGPU()), - * so q_scale.d/q_blocks.d are simply unused there. .l/.u are pinned (cudaMallocHost) - * rather than aligned_alloc when useCuda, see Initialize(), so HtDTransfer()'s async - * uploads of them are genuinely asynchronous instead of silently blocking - * (cudaMemcpyAsync only overlaps with the host from pinned memory); .d is never - * uploaded, so it is always plain aligned_alloc regardless of useCuda. */ + * QuantizeDiagonalBlocks(). .l/.u are pinned (cudaMallocHost) rather than + * aligned_alloc when useCuda, so HtDTransfer()'s async uploads them. */ LDU q_scale; /*!< \brief Quantized block entries; .l/.u sized [nnz_l/u * nVar * nEqn], .d [nPoint * nVar * nEqn]. */ LDU q_blocks; - /*!< \brief Device mirrors of the quantized storage, only allocated when quantized_mode && - * useCuda (currently only reachable for Q_JACOBI/Q_IDENTITY, Q_LU_SGS stays host-only). - * d_q_scale.l/.u and d_q_blocks.l/.u are device-side copies of q_scale.l/.u and q_blocks.l/.u, - * uploaded *asynchronously* by HtDTransfer(). d_q_scale.d/d_q_blocks.d are populated directly - * on the device from gpu.d by QuantizeDiagonalBlocksGPU() instead - gpu.d is uploaded - * unconditionally by HtDTransfer() anyway (Jacobi's own build needs the full precision - * diagonal regardless of quantization), so quantizing it there avoids a host quantize + - * upload round trip for the diagonal specifically. */ + /*!< \brief Device mirrors of the quantized storage, only allocated when quantized_mode. */ LDU d_q_scale; LDU d_q_blocks; @@ -650,8 +620,7 @@ class CSysMatrix { const CConfig* config) const; /*! - * \brief Quantize the diagonal blocks directly on the device, from gpu.d into - * d_q_scale.d/d_q_blocks.d (EncodeQuantBlock, shared with the host path). + * \brief Quantize the diagonal blocks directly on the device. * \note Requires the device matrix to be up to date, see HtDTransfer. */ void QuantizeDiagonalBlocksGPU(); diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index ef49bde5072..8e5baa8ae4c 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -223,10 +223,7 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi * the host, so only plain (or quantized) Jacobi can keep them exclusively on the device. ---*/ jacobi_on_device = useCuda && (prec == JACOBI || prec == Q_JACOBI); #ifndef CODI_REVERSE_TYPE - /*--- Q_LU_SGS stays host-only: its forward/backward sweeps are inherently sequential and - * already bracketed to the host via ApplyPreconditionerOnHost, so there is nothing to gain - * from a device SpMV there. Q_JACOBI and Q_IDENTITY have no such sequential preconditioner - * step, so their (shared) quantized matrix-vector product is allowed on the device too. ---*/ + /*--- Q_LU_SGS is still host-only. ---*/ const bool quantized_offdiag_needed = allow_quant && (prec == Q_JACOBI || prec == Q_IDENTITY || (prec == Q_LU_SGS && !useCuda)); #else @@ -263,12 +260,8 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi #ifndef CODI_REVERSE_TYPE quantized_mode = true; #endif - /*--- .l/.u are pinned (page-locked) when useCuda: HtDTransfer() uploads them with - * cudaMemcpyAsync, which is only genuinely asynchronous from pinned host memory (see - * GPUMemoryAllocation::pinned_alloc); from regular pageable memory it silently degrades to a - * blocking copy. The diagonal (.d) is never uploaded - under CUDA it is quantized straight - * from gpu.d on the device instead (QuantizeDiagonalBlocksGPU()) - so it stays plain - * aligned_alloc regardless of useCuda. ---*/ + /*--- .l/.u are pinned (page-locked) when useCuda because HtDTransfer() uploads them with + * cudaMemcpyAsync, which is only genuinely asynchronous from pinned host memory. ---*/ auto allocQ = [](QuantType*& ptr, unsigned long n) { ptr = MemoryAllocation::aligned_alloc(64, n * sizeof(QuantType)); }; @@ -783,11 +776,7 @@ void CSysMatrix::QuantizeDiagonalBlocks() { /*--- gpu.d is already on the device - HtDTransfer() uploads it unconditionally, since * Jacobi's own build needs the full precision diagonal regardless of quantization - so * quantize straight from it here instead of quantizing on the host and uploading the - * result: one less host/device round trip for the diagonal specifically (L/U still have - * to be quantized on the host, during assembly, since only the host ever touches the - * matrix as it is being built). This calls the exact same EncodeQuantBlock routine the - * host path below uses (CSysMatrix.hpp, SU2_CUDA_HOST_DEVICE), not a separate device copy. - * ---*/ + * result. ---*/ SU2_DEVICE_REGION(QuantizeDiagonalBlocksGPU();) return; } else { diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 243e4faa236..f8795b7a260 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -111,13 +111,7 @@ __device__ FORCEINLINE ScalarType* GetBlockILU(const DeviceLDU& M, u * \brief Parallel Gauss-Jordan matrix inversion, shared by InvertDiagonalBlocksKernel and * IluFactorColorKernel's diagonal-inversion step below - both already have nVar*nVar * threads and two blockSize-sized shared buffers on hand at the point they need a diagonal - * block inverted. Unlike SU2_LinAlg::MatrixInverse (host, and the one still used to invert - * Linelet's tridiagonal blocks) - which does forward elimination on one thread followed by - * a serial back-substitution, an inherently sequential ~nVar^3 chain - this eliminates each - * pivot column from every OTHER row simultaneously (both above and below the pivot, no - * back-substitution needed), so every entry update within one pivot step is independent of - * every other entry in that step: nVar serial steps (each with a couple of barriers), - * instead of one thread serially doing all the work. + * block inverted. SU2_LinAlg::MatrixInverse would have to run on a single thread. * \param i,j Row/column of the block entry this thread owns, in 0..nVar-1 (i.e. threadIdx.x's * divmod by nVar, same mapping the caller already uses for everything else). * \param A Destroyed. \param Inv Must be pre-loaded with the identity, must not alias \p A; @@ -191,15 +185,8 @@ __global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nV /*! * \brief Quantize the diagonal blocks straight from the device diagonal (gpu.d), device * counterpart of CSysMatrix::QuantizeBlock (CSysMatrix.cpp). Calls the exact same - * EncodeQuantRow encoding routine (CSysMatrix.hpp, SU2_CUDA_HOST_DEVICE) that the host - * path's EncodeQuantBlock loops over, rather than a separate device copy of the encoding - * logic. One thread per (point, row) - each block-row's scale and quantization are - * already independent of every other row (see EncodeQuantRow), so this is nVar times - * more parallel than one-thread-per-point, and consecutive threads land nVar elements - * apart in mat_d instead of nVar*nVar apart, a partial coalescing win the same way - * ApplyJacobiPreconditionerKernel's row-major thread mapping is. Points are batched into - * blocks of ~128 threads the same way (threadIdx.x -> (point-within-block, row) via - * divmod by nVar). + * EncodeQuantRow encoding function. One thread per (point, row) - each block-row's + * scale and quantization are independent of every other row. */ template __global__ void QuantizeDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, @@ -459,14 +446,8 @@ __global__ void BlockLDU_SpMV_kernel(unsigned long nRows, unsigned long nVar, } /*! - * \brief Quantized block-LDU SpMV kernel: y[iRow] = (L + D + U) * x per block-row, reading int8 - * row-scaled quantized blocks instead of full precision ones. Device version of - * QuantizedRowProduct/QuantizedMatVecAdd (CSysMatrix.inl). Rows are batched per block the - * same way as BlockLDU_SpMV_kernel (see its comment). Each row of quantized mantissas is - * also read 4 bytes at a time (one 32-bit load instead of four 8-bit ones): dp4a does not - * apply here since only the matrix is quantized, not x, so the accumulation itself is - * still one scalar FMA per element in the exact same order as the scalar loop (bit- - * identical result) - the win is fewer, wider load instructions, not fewer FMAs. + * \brief Device version of QuantizedRowProduct/QuantizedMatVecAdd (CSysMatrix.inl). Rows are + * batched per block the same way as BlockLDU_SpMV_kernel. */ template __global__ void QuantizedBlockLDU_SpMV_kernel( @@ -722,18 +703,9 @@ void CSysMatrix::HtDTransfer(bool trigger) const { /*--- No gpu.l/gpu.u to transfer (never allocated); mirror the host quantized off-diagonal * storage instead (the diagonal mirrors, d_q_scale.d/d_q_blocks.d, are not touched here at * all - QuantizeDiagonalBlocksGPU() populates them straight from gpu.d, just uploaded above, - * with no host round trip). Issued as async copies on a dedicated stream (not the default - * one), so this transfer (the copy engine) can actually run concurrently with whatever the - * preconditioner's Build() launches next on the default stream (the SMs), e.g. - * QuantizeDiagonalBlocksGPU - on the default stream they would just queue behind these - * copies instead of overlapping, since a single stream is strictly ordered. This is only - * genuinely asynchronous (i.e. the host thread does not block here waiting for the copy) - * because q_scale.l/q_blocks.l/q_scale.u/q_blocks.u are pinned host memory, see the comment - * on those members / Initialize() - cudaMemcpyAsync silently degrades to a blocking copy - * from regular pageable memory. htd_event marks the end of this transfer: the first - * default-stream kernel to actually read d_q_scale.l/.u/d_q_blocks.l/.u (the quantized SpMV, - * MatrixVectorProductGPU) waits on it there, since cross-stream dependencies are not - * implied by stream ordering the way same-stream ones are. ---*/ + * with no host round trip). Issued as async copies on a dedicated stream so this transfer can + * run concurrently with whatever the preconditioner's Build() launches next on the default + * stream. ---*/ if (htd_stream == nullptr) gpuErrChk(cudaStreamCreate(&htd_stream)); if (htd_event == nullptr) gpuErrChk(cudaEventCreateWithFlags(&htd_event, cudaEventDisableTiming)); gpuErrChk(cudaMemcpyAsync(d_q_scale.l, q_scale.l, sizeof(QuantType) * mat.nnz_l * nVar, cudaMemcpyHostToDevice, @@ -768,10 +740,7 @@ void CSysMatrix::MatrixVectorProductGPU(const CSysVector dim3 gridDim(static_cast((nPointDomain + rowsPerBlock - 1) / rowsPerBlock), 1, 1); if (quantized_mode) { /*--- Wait (on the device, no host block) for HtDTransfer's async L/U copy on its own stream - * to finish before this default-stream kernel reads d_q_scale.l/.u/d_q_blocks.l/.u - a - * cross-stream dependency, so it is not implied by ordering the way same-stream launches are. - * htd_event is always valid here: HtDTransfer runs once per solve before the first call to - * this function (see its comment), so it has already recorded the event at least once. ---*/ + * to finish before this default-stream kernel reads d_q_scale.l/.u/d_q_blocks.l/.u. ---*/ gpuErrChk(cudaStreamWaitEvent(nullptr, htd_event, 0)); QuantizedBlockLDU_SpMV_kernel<<>>( nPointDomain, nVar, gpu.row_ptr_l, gpu.col_ind_l, d_q_scale.l, d_q_blocks.l, d_q_scale.d, d_q_blocks.d, From 81fc7d6af16cd2bc3b26d0f8a783e8041d5a0a78 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 21:21:36 -0700 Subject: [PATCH 23/30] Fix forward-mode AD build: restore PassiveValue in EncodeQuantRow EncodeQuantRow cast F's return straight to float, on the assumption that quantization is compile-time disabled whenever ScalarType could be AD-active. That's only true for reverse-mode AD - quantized_mode (CSysMatrix.cpp's quantized_offdiag_needed) is gated on `#ifndef CODI_REVERSE_TYPE`, so it's still compiled in for forward-mode AD (CODI_FORWARD_TYPE), where ScalarType is codi::ActiveType<...> and a direct static_cast doesn't compile. Branch on __CUDA_ARCH__ instead: on device (only ever instantiated for plain ScalarType, never AD-active) cast directly; on host go through SU2_TYPE::PassiveValue first (host-only, not SU2_CUDA_HOST_DEVICE) to safely extract the passive value before the float cast. Verified: SU2_CFD_DIRECTDIFF (forward-mode AD) now builds, and the GPU regression set (Jacobi/Q_Jacobi/Q_Identity/ILU) is unaffected - bit-identical to before this fix. Co-Authored-By: Claude Sonnet 5 --- Common/include/linear_algebra/CSysMatrix.hpp | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index bf37ace7c16..ccd9cae7877 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -154,15 +154,25 @@ SU2_CUDA_HOST_DEVICE inline float DecodeQuantScale(int8_t e) noexcept { * \brief Encode one row of an nVar×nVar block into int8 quantized storage: \p qs receives the * row's scale exponent, \p qv (nVar entries) the clamped int8 values for row \p r. * \p f(r,c) is called twice per entry (max-abs scan then encoding); it should be cheap. - * \note Shared with the device and thus same __CUDA_ARCH__ branches as DecodeQuantScale. + * \note Shared with the device and thus same __CUDA_ARCH__ branches as DecodeQuantScale. \p f's + * return type is cast to float directly on device (only ever instantiated there for plain + * ScalarType, never AD-active); on host it goes through SU2_TYPE::PassiveValue first, since + * ScalarType can be AD-active there (quantized_mode is only compiled out for reverse-mode + * AD, not forward-mode, see quantized_offdiag_needed in CSysMatrix.cpp) and PassiveValue is + * host-only (not SU2_CUDA_HOST_DEVICE). */ template SU2_CUDA_HOST_DEVICE inline void EncodeQuantRow(const F& f, int8_t& qs, int8_t* __restrict qv, unsigned long nVar, unsigned long r) noexcept { +#ifdef __CUDA_ARCH__ + auto passive = [&](unsigned long row, unsigned long col) { return f(row, col); }; +#else + auto passive = [&](unsigned long row, unsigned long col) { return SU2_TYPE::PassiveValue(f(row, col)); }; +#endif constexpr uint32_t eps_bits = 0x34000000u; uint32_t max_abs_bits = eps_bits; for (auto c = 0ul; c < nVar; ++c) { - const float fv = static_cast(f(r, c)); + const float fv = static_cast(passive(r, c)); #ifdef __CUDA_ARCH__ const uint32_t fb = __float_as_uint(fv); #else @@ -181,7 +191,8 @@ SU2_CUDA_HOST_DEVICE inline void EncodeQuantRow(const F& f, int8_t& qs, int8_t* memcpy(&inv_rscale, &inv_bits, sizeof(inv_rscale)); #endif for (auto c = 0ul; c < nVar; ++c) { - qv[c] = static_cast(QuantMax(-128.f, QuantMin(127.f, roundf(static_cast(f(r, c)) * inv_rscale)))); + qv[c] = + static_cast(QuantMax(-128.f, QuantMin(127.f, roundf(static_cast(passive(r, c)) * inv_rscale)))); } } From da428d02b9e864435c5cf2a2053957cfa6e2e567 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 22:17:53 -0700 Subject: [PATCH 24/30] Replace multiDotGPU's cublasgemmBatched with a fused reduction kernel multiDotGPU computed the n*m matrix of dot products C(i,j) = via cublasgemmBatched, treating each (i,j) pair as its own independent (M=1, N=1, K=size) GEMM. Measured at ~12ms for a modest n*m, roughly an order of magnitude above the bandwidth-bound ideal - GEMM kernels parallelize across output (M,N) tiles and iterate the K (reduction) dimension within one thread block's launch, so an M=N=1 tile with K in the tens of millions gets little of the parallelism a dedicated reduction primitive would, and the batched form also re-read every input vector once per (i,j) pair instead of once total (n*m reads of length size instead of n+m). The setup also did five blocking cudaMemcpy round trips per call (three pointer-array uploads, a memset, one result download), serializing host and device for no reason. Replace it with a hand-written MultiDotKernel: threads grid-stride over the shared K=size dimension, each accumulating its own private n*m running sums (reading every V/W vector element exactly once, reused across all n*m pairs), then combine via warp-shuffle reduction, one shared-memory combine per block, and one atomicAdd per (i,j) per block - O(nm * numBlocks) atomics instead of O(nm * size). The per-thread accumulator needs a compile-time size (CUDA has no runtime-sized register/local arrays), capped at MULTIDOT_MAX_NM (1024, comfortably above FGCRODR's n=m+1 deflation-matrix usage even for a generously raised LINEAR_SOLVER_RESTART_DEFLATION) with a clear error if exceeded rather than silent truncation. Also switched the remaining setup (pointer-array uploads, output zeroing, result download) to a dedicated stream with async copies instead of five blocking calls. Verified: GPU FGCRODR (which exercises multiDot for both the classical Gram-Schmidt orthogonalization and the Ritz-value deflation matrix, up to n*m=56 at LINEAR_SOLVER_RESTART_DEFLATION=7) matches CPU FGCRODR exactly. The standard GPU regression set (Jacobi/Q_Jacobi/ILU) is unaffected; Q_Identity shows a small (~0.1%) trajectory difference from CPU by iteration 39, the same order of magnitude as the pre-existing GPU(float)/CPU(double) gap seen throughout this branch's work and consistent with this kernel's different (but equally valid) floating- point summation order, not a correctness bug - not helped by this being an unconverged, chaotically-sensitive nonlinear residual at that point. Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysVectorGPU.cu | 226 ++++++++++++++------- 1 file changed, 149 insertions(+), 77 deletions(-) diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index ba53b0cbc29..8d3eec0b34c 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -47,6 +47,83 @@ cublasHandle_t GetBlasHandle() { return solver_blas_handle; } +/*! + * \brief Cap on n*m (the number of dot products one MultiDotKernel launch computes). Bounds the + * per-thread accumulator array below (necessarily a compile-time size - CUDA has no + * runtime-sized register/local arrays) and the block's shared-memory reduction buffer. + * FGCRODR's Ritz-value path (CSysSolve.cpp) calls this with n = m+1 where m is + * LINEAR_SOLVER_RESTART_DEFLATION (user-configurable, default 4 but not uncommon to raise + * into the tens), so n*m grows quadratically with that setting - e.g. m=10 already needs + * 110. Sized generously above that; raise it further if a caller ever needs more (a clear + * SU2_MPI::Error fires instead of silently producing a wrong/truncated result). + */ +constexpr unsigned int MULTIDOT_MAX_NM = 1024; + +/*! + * \brief Compute the n*m matrix of dot products C(i,j) = in a single pass over the + * data, replacing a previous cublasgemmBatched implementation that treated each (i,j) + * pair as its own independent (M=1, N=1, K=size) GEMM. That was a poor fit twice over: + * (a) it re-read every input vector once per pair instead of once total (n*m reads of + * length size instead of n+m), and (b) GEMM kernels parallelize across output (M,N) + * tiles and iterate the reduction (K) dimension within one thread block's kernel + * invocation - for an M=N=1 tile with K in the tens of millions, that leaves most of the + * GPU idle, unlike a dedicated reduction primitive (cublasdot, or this kernel's + * grid-stride loop across many blocks) which parallelizes across all of K. + * \note Each thread accumulates its own private running n*m sums while striding over k (kept in + * thread-local storage, capped at MULTIDOT_MAX_NM - small enough that it should stay + * resident in registers or L1 for realistic n*m, cheap either way next to the K-length + * main loop's DRAM traffic), so the single pass over V/W stays memory-bound. Only after + * that loop do threads combine: a warp-shuffle reduction, then one shared-memory combine + * and one atomicAdd per (i,j) per block, so total atomics are O(nm * numBlocks), not + * O(nm * size). + */ +template +__global__ void MultiDotKernel(const ScalarType* const* __restrict__ V, unsigned int n, + const ScalarType* const* __restrict__ W, unsigned int m, unsigned long size, + ScalarType* __restrict__ C) { + const unsigned int nm = n * m; + + ScalarType local[MULTIDOT_MAX_NM]; + for (unsigned int t = 0; t < nm; ++t) local[t] = ScalarType(0); + + const auto stride = static_cast(blockDim.x) * gridDim.x; + for (auto k = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; k < size; k += stride) { + for (unsigned int i = 0; i < n; ++i) { + const ScalarType v = V[i][k]; + for (unsigned int j = 0; j < m; ++j) local[i * m + j] += v * W[j][k]; + } + } + + /*--- Warp-level reduction: after this, lane 0 of every warp holds that warp's true sum. ---*/ + for (unsigned int t = 0; t < nm; ++t) { + ScalarType val = local[t]; + for (int offset = 16; offset > 0; offset >>= 1) val += __shfl_down_sync(0xFFFFFFFFu, val, offset); + local[t] = val; + } + + /*--- Block-level combine: warp leaders stage their partials in shared memory, thread 0 sums + * them and issues the one atomicAdd per (i,j) this block contributes to the global result + * (pre-zeroed by the caller). ---*/ + extern __shared__ __align__(sizeof(double)) char smem[]; + auto* warpPartials = reinterpret_cast(smem); + const unsigned int lane = threadIdx.x % 32u; + const unsigned int warpId = threadIdx.x / 32u; + const unsigned int warpsPerBlock = blockDim.x / 32u; + + if (lane == 0) { + for (unsigned int t = 0; t < nm; ++t) warpPartials[warpId * nm + t] = local[t]; + } + __syncthreads(); + + if (threadIdx.x == 0) { + for (unsigned int t = 0; t < nm; ++t) { + ScalarType sum = 0; + for (unsigned int w = 0; w < warpsPerBlock; ++w) sum += warpPartials[w * nm + t]; + atomicAdd(&C[t], sum); + } + } +} + } // namespace namespace VecExpr { @@ -112,102 +189,97 @@ ScalarType CSysVector::dotGPU(const CSysVector& other) const { } /*! - * \brief multi vector product with cublasgemmBatched + * \brief Multi vector dot product, C(i,j) = , via MultiDotKernel (see its comment + * for why that beats the batched-GEMM approach this replaced). */ template su2matrix CSysVector::multiDotGPU(const std::vector>& V, const size_t i0, const size_t n, const std::vector>& W, const size_t m) { - /*--- The multiDot product between n V[size] and m W[size] vectors is performed as - * a General Matrix Multiplication between two tall-skinny matrices: - * C = \alpha * A^T * B + \beta * C - * being A = V[ size * n ] and B = W[ size * m ] the batched vectors ---*/ - cublasHandle_t handle = GetBlasHandle(); - cublasStatus_t status = CUBLAS_STATUS_SUCCESS; + const size_t nm = n * m; + if (nm > MULTIDOT_MAX_NM) { + SU2_MPI::Error("CSysVector::multiDotGPU: n*m exceeds MULTIDOT_MAX_NM, raise that constant in " + "CSysVectorGPU.cu.", + CURRENT_FUNCTION); + } const size_t size = V[0].nElmDomain; - const size_t batch = n * m; - /*--- Persistent device workspace, cached across calls and freed automatically - * when the program exits (static local destruction), instead of leaking. ---*/ + su2matrix local; + local.resize(n, m); + if (nm == 0) return local; + + /*--- Persistent device workspace, cached across calls and freed automatically when the + * program exits (static local destruction), instead of leaking. Only n+m pointers (not n*m, + * as the previous per-pair-GEMM implementation needed) since MultiDotKernel reads each vector + * once and reuses it across all pairs internally. ---*/ struct Workspace { - ScalarType* d_local = nullptr; - const ScalarType** d_A = nullptr; - const ScalarType** d_B = nullptr; - ScalarType** d_C = nullptr; - size_t capacity = 0; - - void EnsureCapacity(size_t batch) { - if (batch <= capacity) return; - cudaFree(d_local); - cudaFree(d_A); - cudaFree(d_B); - cudaFree(d_C); - gpuErrChk(cudaMalloc(&d_local, batch * sizeof(ScalarType))); - gpuErrChk(cudaMalloc(&d_A, batch * sizeof(ScalarType*))); - gpuErrChk(cudaMalloc(&d_B, batch * sizeof(ScalarType*))); - gpuErrChk(cudaMalloc(&d_C, batch * sizeof(ScalarType*))); - capacity = batch; + const ScalarType** d_V = nullptr; + const ScalarType** d_W = nullptr; + ScalarType* d_C = nullptr; + cudaStream_t stream = nullptr; + size_t vCapacity = 0, wCapacity = 0, cCapacity = 0; + + void EnsureCapacity(size_t n, size_t m, size_t nm) { + if (stream == nullptr) gpuErrChk(cudaStreamCreate(&stream)); + if (n > vCapacity) { + cudaFree(d_V); + gpuErrChk(cudaMalloc(&d_V, n * sizeof(ScalarType*))); + vCapacity = n; + } + if (m > wCapacity) { + cudaFree(d_W); + gpuErrChk(cudaMalloc(&d_W, m * sizeof(ScalarType*))); + wCapacity = m; + } + if (nm > cCapacity) { + cudaFree(d_C); + gpuErrChk(cudaMalloc(&d_C, nm * sizeof(ScalarType))); + cCapacity = nm; + } } ~Workspace() { - cudaFree(d_local); - cudaFree(d_A); - cudaFree(d_B); + cudaFree(d_V); + cudaFree(d_W); cudaFree(d_C); + if (stream != nullptr) cudaStreamDestroy(stream); } }; static Workspace ws; - - // allocate persistent result buffer local on host and device, is resized if needed - su2matrix local; - local.resize(n, m); - ws.EnsureCapacity(batch); - - // zero out the result buffer - gpuErrChk(cudaMemset(ws.d_local, 0, batch * sizeof(ScalarType))); - - // prepare the arrays A,B,C on host - static std::vector h_A, h_B; - static std::vector h_C; - h_A.resize(batch); h_B.resize(batch); h_C.resize(batch); - - for (size_t i = 0; i < n; ++i) { - for (size_t j =0; j < m; ++j) { - const size_t idx = i * m + j; - h_A[idx] = V[i0 + i].GetDevicePointer(); - h_B[idx] = W[j].GetDevicePointer(); - h_C[idx] = ws.d_local + idx; // C maps to d_local to store the coefficients in the 2D array - } + ws.EnsureCapacity(n, m, nm); + + static std::vector h_V, h_W; + h_V.resize(n); + h_W.resize(m); + for (size_t i = 0; i < n; ++i) h_V[i] = V[i0 + i].GetDevicePointer(); + for (size_t j = 0; j < m; ++j) h_W[j] = W[j].GetDevicePointer(); + + gpuErrChk(cudaMemcpyAsync(ws.d_V, h_V.data(), n * sizeof(ScalarType*), cudaMemcpyHostToDevice, ws.stream)); + gpuErrChk(cudaMemcpyAsync(ws.d_W, h_W.data(), m * sizeof(ScalarType*), cudaMemcpyHostToDevice, ws.stream)); + gpuErrChk(cudaMemsetAsync(ws.d_C, 0, nm * sizeof(ScalarType), ws.stream)); + + constexpr unsigned int threadsPerBlock = 256; + const auto blocks = static_cast(std::min((size + threadsPerBlock - 1) / threadsPerBlock, 1024)); + const auto sharedBytes = static_cast(threadsPerBlock / 32u) * nm * sizeof(ScalarType); + + /*--- sharedBytes can exceed the default 48KB static shared-memory limit for large nm (the + * FGCRODR deflation matrix in particular, see MULTIDOT_MAX_NM); opt in to the device's larger + * "dynamic" limit once, the first time it is actually needed, rather than always paying for + * the query. ---*/ + static size_t optedInSharedBytes = 0; + if (sharedBytes > optedInSharedBytes) { + gpuErrChk(cudaFuncSetAttribute(MultiDotKernel, cudaFuncAttributeMaxDynamicSharedMemorySize, + static_cast(sharedBytes))); + optedInSharedBytes = sharedBytes; } - // copy pointers to device - gpuErrChk(cudaMemcpy(ws.d_A, h_A.data(), batch * sizeof(ScalarType*), cudaMemcpyHostToDevice)); - gpuErrChk(cudaMemcpy(ws.d_B, h_B.data(), batch * sizeof(ScalarType*), cudaMemcpyHostToDevice)); - gpuErrChk(cudaMemcpy(ws.d_C, h_C.data(), batch * sizeof(ScalarType*), cudaMemcpyHostToDevice)); - - // define alpha = 1.0 and beta = 0.0 - const auto alpha = ScalarType(1.0); - const auto beta = ScalarType(0.0); - - if constexpr (std::is_same_v) { - status = cublasSgemmBatched(handle, CUBLAS_OP_T, CUBLAS_OP_N, 1, 1, size, &alpha, ws.d_A, static_cast(size), - ws.d_B, static_cast(size), &beta, ws.d_C, 1, static_cast(batch)); - } else if constexpr (std::is_same_v) { - status = cublasDgemmBatched(handle, CUBLAS_OP_T, CUBLAS_OP_N, 1, 1, size, &alpha, ws.d_A, static_cast(size), - ws.d_B, static_cast(size), &beta, ws.d_C, 1, static_cast(batch)); - } else { - SU2_MPI::Error("Unsupported ScalarType in CSysVector::multiDotGPU.", CURRENT_FUNCTION); - return local; - } - - if (status != CUBLAS_STATUS_SUCCESS) { - SU2_MPI::Error("cuBLAS cublasgemmBatched failed in CSysVector::multiDotGPU.", CURRENT_FUNCTION); - return local; - } + MultiDotKernel<<>>( + ws.d_V, static_cast(n), ws.d_W, static_cast(m), size, ws.d_C); - // copy result to host for MPI reduce - gpuErrChk(cudaMemcpy(local.data(), ws.d_local, batch * sizeof(ScalarType), cudaMemcpyDeviceToHost)); + gpuErrChk(cudaMemcpyAsync(local.data(), ws.d_C, nm * sizeof(ScalarType), cudaMemcpyDeviceToHost, ws.stream)); + gpuErrChk(cudaStreamSynchronize(ws.stream)); + gpuErrChk(cudaGetLastError()); return local; } From e52584800ea1f7ab9b659b547ed55490c576aa90 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 15 Aug 2026 22:45:25 -0700 Subject: [PATCH 25/30] Pass MultiDotKernel's V/W pointers as launch parameters, not a device array The previous version uploaded the n and m device pointers into device buffers (d_V/d_W) via cudaMemcpyAsync from plain (pageable) host std::vectors. That's the same trap called out for HtDTransfer's L/U copy: cudaMemcpyAsync from pageable host memory silently degrades to a blocking copy, so the "async" upload wasn't actually async without also pinning h_V/h_W. Pass the pointers as ordinary by-value kernel launch parameters instead (MultiDotPointers, a fixed-size array wrapper) - CUDA already marshals launch parameters as part of the launch itself, regardless of host memory pinning, so this sidesteps the whole question rather than fixing it with a pinned buffer. Drops the d_V/d_W device buffers and their setup entirely; only the output accumulator (d_C) still needs device storage, since atomicAdd needs a device-resident target. This needed a cap on n and m individually (MULTIDOT_MAX_VEC, 128) to bound the fixed-size parameter struct, on top of the existing cap on n*m (MULTIDOT_MAX_NM, bounding the per-thread accumulator array) - a different constraint, since n or m can independently reach the Krylov restart length depending on call site (ModGramSchmidt grows m with n=1; FGCRODR's Ritz-value path grows n with m bounded by the deflation count), so both n*m and n,m individually need their own headroom. Verified: GPU FGCRODR (LINEAR_SOLVER_RESTART_DEFLATION=7, exercising both ModGramSchmidt and the Ritz-value path) matches CPU FGCRODR exactly, reproducible over repeated runs. Standard GPU regression set (Jacobi/Q_Jacobi/Q_Identity/ILU) unaffected. Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysVectorGPU.cu | 87 +++++++++++++--------- 1 file changed, 52 insertions(+), 35 deletions(-) diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 8d3eec0b34c..20fdf0767fe 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -55,10 +55,39 @@ cublasHandle_t GetBlasHandle() { * LINEAR_SOLVER_RESTART_DEFLATION (user-configurable, default 4 but not uncommon to raise * into the tens), so n*m grows quadratically with that setting - e.g. m=10 already needs * 110. Sized generously above that; raise it further if a caller ever needs more (a clear - * SU2_MPI::Error fires instead of silently producing a wrong/truncated result). + * SU2_MPI::Error fires instead of silently producing a wrong/truncated result). This is a + * different constraint from MULTIDOT_MAX_VEC below (n and m individually can each reach + * that cap, but not both at once - n*m up to MULTIDOT_MAX_VEC^2 would need an + * impractically large per-thread array). */ constexpr unsigned int MULTIDOT_MAX_NM = 1024; +/*! + * \brief Cap on the individual vector counts n and m. Bounds the fixed-size pointer arrays + * passed into MultiDotKernel by value as ordinary launch parameters, rather than as a + * device-side array of pointers uploaded via a separate cudaMemcpy - CUDA already + * marshals launch parameters for you as part of the launch itself, regardless of host + * memory pinning, which sidesteps needing pinned host buffers here the way HtDTransfer's + * L/U copy does (see aux_stream/htd_event in CSysMatrix, CUDA silently downgrades + * cudaMemcpyAsync to a blocking copy from ordinary pageable host memory). n or m can + * independently reach the Krylov restart length in this codebase's two call sites + * (ModGramSchmidt in CSysSolve.cpp grows m up to the restart length with n=1; FGCRODR's + * Ritz-value path grows n up to the restart length with m bounded by + * LINEAR_SOLVER_RESTART_DEFLATION), so this has to cover restart length, not just the + * (usually much smaller) deflation count - 128 is generous headroom over realistic + * restart lengths (typically 10-50). + */ +constexpr unsigned int MULTIDOT_MAX_VEC = 128; + +/*! + * \brief Fixed-size array of device pointers, passed to MultiDotKernel by value (see + * MULTIDOT_MAX_VEC for why). + */ +template +struct MultiDotPointers { + const ScalarType* ptr[MULTIDOT_MAX_VEC]; +}; + /*! * \brief Compute the n*m matrix of dot products C(i,j) = in a single pass over the * data, replacing a previous cublasgemmBatched implementation that treated each (i,j) @@ -78,9 +107,8 @@ constexpr unsigned int MULTIDOT_MAX_NM = 1024; * O(nm * size). */ template -__global__ void MultiDotKernel(const ScalarType* const* __restrict__ V, unsigned int n, - const ScalarType* const* __restrict__ W, unsigned int m, unsigned long size, - ScalarType* __restrict__ C) { +__global__ void MultiDotKernel(MultiDotPointers V, unsigned int n, MultiDotPointers W, + unsigned int m, unsigned long size, ScalarType* __restrict__ C) { const unsigned int nm = n * m; ScalarType local[MULTIDOT_MAX_NM]; @@ -89,8 +117,8 @@ __global__ void MultiDotKernel(const ScalarType* const* __restrict__ V, unsigned const auto stride = static_cast(blockDim.x) * gridDim.x; for (auto k = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; k < size; k += stride) { for (unsigned int i = 0; i < n; ++i) { - const ScalarType v = V[i][k]; - for (unsigned int j = 0; j < m; ++j) local[i * m + j] += v * W[j][k]; + const ScalarType v = V.ptr[i][k]; + for (unsigned int j = 0; j < m; ++j) local[i * m + j] += v * W.ptr[j][k]; } } @@ -196,6 +224,11 @@ template su2matrix CSysVector::multiDotGPU(const std::vector>& V, const size_t i0, const size_t n, const std::vector>& W, const size_t m) { + if (n > MULTIDOT_MAX_VEC || m > MULTIDOT_MAX_VEC) { + SU2_MPI::Error("CSysVector::multiDotGPU: n or m exceeds MULTIDOT_MAX_VEC, raise that constant in " + "CSysVectorGPU.cu.", + CURRENT_FUNCTION); + } const size_t nm = n * m; if (nm > MULTIDOT_MAX_NM) { SU2_MPI::Error("CSysVector::multiDotGPU: n*m exceeds MULTIDOT_MAX_NM, raise that constant in " @@ -209,29 +242,19 @@ su2matrix CSysVector::multiDotGPU(const std::vector vCapacity) { - cudaFree(d_V); - gpuErrChk(cudaMalloc(&d_V, n * sizeof(ScalarType*))); - vCapacity = n; - } - if (m > wCapacity) { - cudaFree(d_W); - gpuErrChk(cudaMalloc(&d_W, m * sizeof(ScalarType*))); - wCapacity = m; - } if (nm > cCapacity) { cudaFree(d_C); gpuErrChk(cudaMalloc(&d_C, nm * sizeof(ScalarType))); @@ -240,23 +263,17 @@ su2matrix CSysVector::multiDotGPU(const std::vector h_V, h_W; - h_V.resize(n); - h_W.resize(m); - for (size_t i = 0; i < n; ++i) h_V[i] = V[i0 + i].GetDevicePointer(); - for (size_t j = 0; j < m; ++j) h_W[j] = W[j].GetDevicePointer(); + MultiDotPointers vPtrs{}, wPtrs{}; + for (size_t i = 0; i < n; ++i) vPtrs.ptr[i] = V[i0 + i].GetDevicePointer(); + for (size_t j = 0; j < m; ++j) wPtrs.ptr[j] = W[j].GetDevicePointer(); - gpuErrChk(cudaMemcpyAsync(ws.d_V, h_V.data(), n * sizeof(ScalarType*), cudaMemcpyHostToDevice, ws.stream)); - gpuErrChk(cudaMemcpyAsync(ws.d_W, h_W.data(), m * sizeof(ScalarType*), cudaMemcpyHostToDevice, ws.stream)); gpuErrChk(cudaMemsetAsync(ws.d_C, 0, nm * sizeof(ScalarType), ws.stream)); constexpr unsigned int threadsPerBlock = 256; @@ -275,7 +292,7 @@ su2matrix CSysVector::multiDotGPU(const std::vector<<>>( - ws.d_V, static_cast(n), ws.d_W, static_cast(m), size, ws.d_C); + vPtrs, static_cast(n), wPtrs, static_cast(m), size, ws.d_C); gpuErrChk(cudaMemcpyAsync(local.data(), ws.d_C, nm * sizeof(ScalarType), cudaMemcpyDeviceToHost, ws.stream)); gpuErrChk(cudaStreamSynchronize(ws.stream)); From f9de6d47b0c968ffa4469d7ed60227918415ef06 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 16 Aug 2026 09:03:33 -0700 Subject: [PATCH 26/30] Split multiDotGPU's pointer-array cap into an asymmetric large/small pair MULTIDOT_MAX_VEC (128) bounded n and m symmetrically, but the two call sites in this codebase produce genuinely asymmetric shapes: ModGramSchmidt (CSysSolve.cpp) grows m up to the Krylov restart length with n=1 (really a row vector, not a rectangle), while FGCRODR's Ritz-value path grows n up to the restart length with m bounded by the usually much smaller LINEAR_SOLVER_RESTART_DEFLATION - a genuine, independently-sized rectangle where neither side is trivially 1. Split into MULTIDOT_MAX_VEC_LARGE (256) and MULTIDOT_MAX_VEC_SMALL (64) instead of one symmetric cap: multiDotGPU now puts whichever of n, m is larger into MultiDotKernel's "A" (large-capacity) argument - swapping V/W (and their counts) for the call when m > n - and transposes the small result matrix back into the (n,m) shape the caller expects if it had to swap. This gets more headroom on the genuinely-large dimension for about the same total kernel-parameter budget as before (256+64 pointers vs the old 128+128), rather than wasting half the budget covering a "large" case on an axis that in practice only ever needs to hold a handful of deflation vectors. MULTIDOT_MAX_NM (1024, bounding the per-thread accumulator array) is unrelated to this and stays as-is: it's a genuine per-thread storage cost paid regardless of the runtime n*m, unlike the two vector-count caps, which only bound a launch parameter's marshaling size - growing it to cover MULTIDOT_MAX_VEC_LARGE * MULTIDOT_MAX_VEC_SMALL would be an impractically large per-thread footprint for the common (small nm) case. Verified: GPU FGCRODR (LINEAR_SOLVER_RESTART_DEFLATION=7, exercising both the swapped ModGramSchmidt shape and FGCRODR's own Ritz-value shape) matches CPU FGCRODR exactly, reproducible over 4 repeated runs. Standard GPU regression set (Jacobi/Q_Jacobi/Q_Identity/ILU) unaffected. Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysVectorGPU.cu | 175 +++++++++++++-------- 1 file changed, 106 insertions(+), 69 deletions(-) diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 20fdf0767fe..f64a797df01 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -56,69 +56,80 @@ cublasHandle_t GetBlasHandle() { * into the tens), so n*m grows quadratically with that setting - e.g. m=10 already needs * 110. Sized generously above that; raise it further if a caller ever needs more (a clear * SU2_MPI::Error fires instead of silently producing a wrong/truncated result). This is a - * different constraint from MULTIDOT_MAX_VEC below (n and m individually can each reach - * that cap, but not both at once - n*m up to MULTIDOT_MAX_VEC^2 would need an - * impractically large per-thread array). + * different, unrelated constraint from MULTIDOT_MAX_VEC_LARGE/SMALL below - it is a + * genuine per-thread storage cost (every thread pays for the compile-time-sized array + * below regardless of the runtime n*m), unlike those two, which only bound a launch + * parameter's marshaling size, so it must stay independently modest rather than grow to + * cover MULTIDOT_MAX_VEC_LARGE * MULTIDOT_MAX_VEC_SMALL (which would be impractically + * large per thread). */ constexpr unsigned int MULTIDOT_MAX_NM = 1024; /*! - * \brief Cap on the individual vector counts n and m. Bounds the fixed-size pointer arrays - * passed into MultiDotKernel by value as ordinary launch parameters, rather than as a - * device-side array of pointers uploaded via a separate cudaMemcpy - CUDA already - * marshals launch parameters for you as part of the launch itself, regardless of host - * memory pinning, which sidesteps needing pinned host buffers here the way HtDTransfer's - * L/U copy does (see aux_stream/htd_event in CSysMatrix, CUDA silently downgrades - * cudaMemcpyAsync to a blocking copy from ordinary pageable host memory). n or m can - * independently reach the Krylov restart length in this codebase's two call sites - * (ModGramSchmidt in CSysSolve.cpp grows m up to the restart length with n=1; FGCRODR's - * Ritz-value path grows n up to the restart length with m bounded by - * LINEAR_SOLVER_RESTART_DEFLATION), so this has to cover restart length, not just the - * (usually much smaller) deflation count - 128 is generous headroom over realistic - * restart lengths (typically 10-50). + * \brief Caps on the individual vector counts n and m, asymmetric because the two call sites in + * this codebase produce genuinely asymmetric shapes: ModGramSchmidt (CSysSolve.cpp) grows + * m up to the Krylov restart length with n=1 (a row vector, not really a rectangle), while + * FGCRODR's Ritz-value path grows n up to the restart length with m bounded by the + * (usually much smaller) LINEAR_SOLVER_RESTART_DEFLATION - a genuine, independently-sized + * rectangle. The host wrapper always puts whichever of n, m is larger into the "large" + * kernel argument (swapping V/W and transposing the result back if needed, see + * multiDotGPU), so one axis only needs to cover the deflation-count case while the other + * covers the restart-length case, instead of both needing the same generous bound (as a + * single symmetric MULTIDOT_MAX_VEC used to). Bounds the fixed-size pointer arrays passed + * into MultiDotKernel by value as ordinary launch parameters, rather than as a device-side + * array of pointers uploaded via a separate cudaMemcpy - CUDA already marshals launch + * parameters for you as part of the launch itself, regardless of host memory pinning, + * which sidesteps needing pinned host buffers here the way HtDTransfer's L/U copy does + * (see aux_stream/htd_event in CSysMatrix, CUDA silently downgrades cudaMemcpyAsync to a + * blocking copy from ordinary pageable host memory). */ -constexpr unsigned int MULTIDOT_MAX_VEC = 128; +constexpr unsigned int MULTIDOT_MAX_VEC_LARGE = 256; +constexpr unsigned int MULTIDOT_MAX_VEC_SMALL = 64; /*! * \brief Fixed-size array of device pointers, passed to MultiDotKernel by value (see - * MULTIDOT_MAX_VEC for why). + * MULTIDOT_MAX_VEC_LARGE/SMALL for why). */ -template +template struct MultiDotPointers { - const ScalarType* ptr[MULTIDOT_MAX_VEC]; + const ScalarType* ptr[MaxCount]; }; /*! - * \brief Compute the n*m matrix of dot products C(i,j) = in a single pass over the - * data, replacing a previous cublasgemmBatched implementation that treated each (i,j) - * pair as its own independent (M=1, N=1, K=size) GEMM. That was a poor fit twice over: - * (a) it re-read every input vector once per pair instead of once total (n*m reads of - * length size instead of n+m), and (b) GEMM kernels parallelize across output (M,N) - * tiles and iterate the reduction (K) dimension within one thread block's kernel - * invocation - for an M=N=1 tile with K in the tens of millions, that leaves most of the - * GPU idle, unlike a dedicated reduction primitive (cublasdot, or this kernel's - * grid-stride loop across many blocks) which parallelizes across all of K. - * \note Each thread accumulates its own private running n*m sums while striding over k (kept in - * thread-local storage, capped at MULTIDOT_MAX_NM - small enough that it should stay - * resident in registers or L1 for realistic n*m, cheap either way next to the K-length - * main loop's DRAM traffic), so the single pass over V/W stays memory-bound. Only after - * that loop do threads combine: a warp-shuffle reduction, then one shared-memory combine - * and one atomicAdd per (i,j) per block, so total atomics are O(nm * numBlocks), not - * O(nm * size). + * \brief Compute the aCount*bCount matrix of dot products D(a,b) = in a single pass + * over the data, replacing a previous cublasgemmBatched implementation that treated + * each pair as its own independent (M=1, N=1, K=size) GEMM. That was a poor fit twice + * over: (a) it re-read every input vector once per pair instead of once total (aCount* + * bCount reads of length size instead of aCount+bCount), and (b) GEMM kernels parallelize + * across output (M,N) tiles and iterate the reduction (K) dimension within one thread + * block's kernel invocation - for an M=N=1 tile with K in the tens of millions, that + * leaves most of the GPU idle, unlike a dedicated reduction primitive (cublasdot, or + * this kernel's grid-stride loop across many blocks) which parallelizes across all of K. + * \note Each thread accumulates its own private running aCount*bCount sums while striding over k + * (kept in thread-local storage, capped at MULTIDOT_MAX_NM - small enough that it should + * stay resident in registers or L1 for realistic shapes, cheap either way next to the + * K-length main loop's DRAM traffic), so the single pass over A/B stays memory-bound. Only + * after that loop do threads combine: a warp-shuffle reduction, then one shared-memory + * combine and one atomicAdd per (a,b) per block, so total atomics are + * O(aCount*bCount * numBlocks), not O(aCount*bCount * size). + * \note A is the "large" argument (up to MULTIDOT_MAX_VEC_LARGE), B the "small" one (up to + * MULTIDOT_MAX_VEC_SMALL) - the caller (multiDotGPU) is responsible for putting the larger + * of its two vector counts into A, and transposing the result back if it had to swap. */ template -__global__ void MultiDotKernel(MultiDotPointers V, unsigned int n, MultiDotPointers W, - unsigned int m, unsigned long size, ScalarType* __restrict__ C) { - const unsigned int nm = n * m; +__global__ void MultiDotKernel(MultiDotPointers A, unsigned int aCount, + MultiDotPointers B, unsigned int bCount, + unsigned long size, ScalarType* __restrict__ D) { + const unsigned int nm = aCount * bCount; ScalarType local[MULTIDOT_MAX_NM]; for (unsigned int t = 0; t < nm; ++t) local[t] = ScalarType(0); const auto stride = static_cast(blockDim.x) * gridDim.x; for (auto k = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; k < size; k += stride) { - for (unsigned int i = 0; i < n; ++i) { - const ScalarType v = V.ptr[i][k]; - for (unsigned int j = 0; j < m; ++j) local[i * m + j] += v * W.ptr[j][k]; + for (unsigned int a = 0; a < aCount; ++a) { + const ScalarType v = A.ptr[a][k]; + for (unsigned int b = 0; b < bCount; ++b) local[a * bCount + b] += v * B.ptr[b][k]; } } @@ -147,7 +158,7 @@ __global__ void MultiDotKernel(MultiDotPointers V, unsigned int n, M for (unsigned int t = 0; t < nm; ++t) { ScalarType sum = 0; for (unsigned int w = 0; w < warpsPerBlock; ++w) sum += warpPartials[w * nm + t]; - atomicAdd(&C[t], sum); + atomicAdd(&D[t], sum); } } } @@ -217,18 +228,18 @@ ScalarType CSysVector::dotGPU(const CSysVector& other) const { } /*! - * \brief Multi vector dot product, C(i,j) = , via MultiDotKernel (see its comment - * for why that beats the batched-GEMM approach this replaced). + * \brief Multi vector dot product, local(i,j) = , via MultiDotKernel (see its + * comment for why that beats the batched-GEMM approach this replaced). + * \note Whichever of n, m is larger is passed as MultiDotKernel's "A" (large-capacity) argument + * and the other as "B" (small-capacity) - see MULTIDOT_MAX_VEC_LARGE/SMALL for why - so if + * m > n, V and W (and their counts) are swapped for the call, and the resulting m*n matrix + * is transposed back into the n*m shape the caller expects (cheap: this matrix is at most + * MULTIDOT_MAX_NM elements, nowhere near the size of the reduction itself). */ template su2matrix CSysVector::multiDotGPU(const std::vector>& V, const size_t i0, const size_t n, const std::vector>& W, const size_t m) { - if (n > MULTIDOT_MAX_VEC || m > MULTIDOT_MAX_VEC) { - SU2_MPI::Error("CSysVector::multiDotGPU: n or m exceeds MULTIDOT_MAX_VEC, raise that constant in " - "CSysVectorGPU.cu.", - CURRENT_FUNCTION); - } const size_t nm = n * m; if (nm > MULTIDOT_MAX_NM) { SU2_MPI::Error("CSysVector::multiDotGPU: n*m exceeds MULTIDOT_MAX_NM, raise that constant in " @@ -236,45 +247,60 @@ su2matrix CSysVector::multiDotGPU(const std::vector local; local.resize(n, m); if (nm == 0) return local; + const bool swap = m > n; + const size_t aCount = swap ? m : n; + const size_t bCount = swap ? n : m; + if (aCount > MULTIDOT_MAX_VEC_LARGE || bCount > MULTIDOT_MAX_VEC_SMALL) { + SU2_MPI::Error("CSysVector::multiDotGPU: n or m exceeds MULTIDOT_MAX_VEC_LARGE/SMALL, raise " + "those constants in CSysVectorGPU.cu.", + CURRENT_FUNCTION); + } + + const size_t size = V[0].nElmDomain; + /*--- Persistent device workspace for the output only, cached across calls and freed * automatically when the program exits (static local destruction), instead of leaking. The * V/W pointers themselves need no device buffer at all: they go to MultiDotKernel as ordinary - * by-value launch parameters (see MultiDotPointers/MULTIDOT_MAX_VEC), which CUDA marshals - * itself as part of the launch - unlike a separate cudaMemcpy, that needs no pinned host - * buffer to actually be asynchronous. ---*/ + * by-value launch parameters (see MultiDotPointers/MULTIDOT_MAX_VEC_LARGE/SMALL), which CUDA + * marshals itself as part of the launch - unlike a separate cudaMemcpy, that needs no pinned + * host buffer to actually be asynchronous. ---*/ struct Workspace { - ScalarType* d_C = nullptr; + ScalarType* d_D = nullptr; cudaStream_t stream = nullptr; - size_t cCapacity = 0; + size_t capacity = 0; void EnsureCapacity(size_t nm) { if (stream == nullptr) gpuErrChk(cudaStreamCreate(&stream)); - if (nm > cCapacity) { - cudaFree(d_C); - gpuErrChk(cudaMalloc(&d_C, nm * sizeof(ScalarType))); - cCapacity = nm; + if (nm > capacity) { + cudaFree(d_D); + gpuErrChk(cudaMalloc(&d_D, nm * sizeof(ScalarType))); + capacity = nm; } } ~Workspace() { - cudaFree(d_C); + cudaFree(d_D); if (stream != nullptr) cudaStreamDestroy(stream); } }; static Workspace ws; ws.EnsureCapacity(nm); - MultiDotPointers vPtrs{}, wPtrs{}; - for (size_t i = 0; i < n; ++i) vPtrs.ptr[i] = V[i0 + i].GetDevicePointer(); - for (size_t j = 0; j < m; ++j) wPtrs.ptr[j] = W[j].GetDevicePointer(); + MultiDotPointers aPtrs{}; + MultiDotPointers bPtrs{}; + if (!swap) { + for (size_t i = 0; i < n; ++i) aPtrs.ptr[i] = V[i0 + i].GetDevicePointer(); + for (size_t j = 0; j < m; ++j) bPtrs.ptr[j] = W[j].GetDevicePointer(); + } else { + for (size_t i = 0; i < m; ++i) aPtrs.ptr[i] = W[i].GetDevicePointer(); + for (size_t j = 0; j < n; ++j) bPtrs.ptr[j] = V[i0 + j].GetDevicePointer(); + } - gpuErrChk(cudaMemsetAsync(ws.d_C, 0, nm * sizeof(ScalarType), ws.stream)); + gpuErrChk(cudaMemsetAsync(ws.d_D, 0, nm * sizeof(ScalarType), ws.stream)); constexpr unsigned int threadsPerBlock = 256; const auto blocks = static_cast(std::min((size + threadsPerBlock - 1) / threadsPerBlock, 1024)); @@ -292,10 +318,21 @@ su2matrix CSysVector::multiDotGPU(const std::vector<<>>( - vPtrs, static_cast(n), wPtrs, static_cast(m), size, ws.d_C); + aPtrs, static_cast(aCount), bPtrs, static_cast(bCount), size, ws.d_D); - gpuErrChk(cudaMemcpyAsync(local.data(), ws.d_C, nm * sizeof(ScalarType), cudaMemcpyDeviceToHost, ws.stream)); - gpuErrChk(cudaStreamSynchronize(ws.stream)); + if (!swap) { + gpuErrChk(cudaMemcpyAsync(local.data(), ws.d_D, nm * sizeof(ScalarType), cudaMemcpyDeviceToHost, ws.stream)); + gpuErrChk(cudaStreamSynchronize(ws.stream)); + } else { + /*--- D is aCount*bCount = m*n row-major (D(a,b) = ); local is n*m with + * local(i,j) = = D(j,i), i.e. local is D transposed. ---*/ + static std::vector D; + D.resize(nm); + gpuErrChk(cudaMemcpyAsync(D.data(), ws.d_D, nm * sizeof(ScalarType), cudaMemcpyDeviceToHost, ws.stream)); + gpuErrChk(cudaStreamSynchronize(ws.stream)); + for (size_t i = 0; i < n; ++i) + for (size_t j = 0; j < m; ++j) local(i, j) = D[j * n + i]; + } gpuErrChk(cudaGetLastError()); return local; From e6554d0044cc9b1b8d4caf885378433993fb6641 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 16 Aug 2026 09:10:45 -0700 Subject: [PATCH 27/30] Use the default stream in multiDotGPU instead of a dedicated one Not just unnecessary overhead - a real latent race. V/W are written by VecExpr::AssignDeviceExpression (e.g. LinearCombination building w[i+1] right before ModGramSchmidt calls multiDot on it), which launches its kernel on the default stream with no synchronization of its own. multiDotGPU's kernel ran on a separate dedicated stream, so there was no guaranteed ordering between that write and this read - CUDA only orders operations within the same stream, not across different ones. It happened not to show up in testing (probably launch overhead vs. kernel runtime timing on this hardware/driver), but nothing in the CUDA execution model guaranteed that. There is also nothing to gain from keeping a dedicated stream now that the V/W pointer transfer is by-value launch parameters rather than an async cudaMemcpy: this function always synchronizes before returning (it needs to hand the caller real values), and its result buffer isn't pinned, so there was no async benefit being given up by using the default stream. Verified: GPU FGCRODR (LINEAR_SOLVER_RESTART_DEFLATION=7) still matches CPU FGCRODR exactly, reproducible over 5 repeated runs. Standard GPU regression set (Jacobi/Q_Jacobi/Q_Identity/ILU) unaffected. Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysVectorGPU.cu | 30 ++++++++++++---------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index f64a797df01..7f669cb31ad 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -267,14 +267,21 @@ su2matrix CSysVector::multiDotGPU(const std::vector capacity) { cudaFree(d_D); gpuErrChk(cudaMalloc(&d_D, nm * sizeof(ScalarType))); @@ -282,10 +289,7 @@ su2matrix CSysVector::multiDotGPU(const std::vector CSysVector::multiDotGPU(const std::vector(std::min((size + threadsPerBlock - 1) / threadsPerBlock, 1024)); @@ -317,19 +321,19 @@ su2matrix CSysVector::multiDotGPU(const std::vector<<>>( + MultiDotKernel<<>>( aPtrs, static_cast(aCount), bPtrs, static_cast(bCount), size, ws.d_D); if (!swap) { - gpuErrChk(cudaMemcpyAsync(local.data(), ws.d_D, nm * sizeof(ScalarType), cudaMemcpyDeviceToHost, ws.stream)); - gpuErrChk(cudaStreamSynchronize(ws.stream)); + gpuErrChk(cudaMemcpyAsync(local.data(), ws.d_D, nm * sizeof(ScalarType), cudaMemcpyDeviceToHost)); + gpuErrChk(cudaStreamSynchronize(nullptr)); } else { /*--- D is aCount*bCount = m*n row-major (D(a,b) = ); local is n*m with * local(i,j) = = D(j,i), i.e. local is D transposed. ---*/ static std::vector D; D.resize(nm); - gpuErrChk(cudaMemcpyAsync(D.data(), ws.d_D, nm * sizeof(ScalarType), cudaMemcpyDeviceToHost, ws.stream)); - gpuErrChk(cudaStreamSynchronize(ws.stream)); + gpuErrChk(cudaMemcpyAsync(D.data(), ws.d_D, nm * sizeof(ScalarType), cudaMemcpyDeviceToHost)); + gpuErrChk(cudaStreamSynchronize(nullptr)); for (size_t i = 0; i < n; ++i) for (size_t j = 0; j < m; ++j) local(i, j) = D[j * n + i]; } From d0e94996f9972cfeb46b4c824d1639ba898c579a Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 16 Aug 2026 09:18:07 -0700 Subject: [PATCH 28/30] Replace dotGPU's cuBLAS call with a custom kernel, drop the cuBLAS dependency dotGPU was the last user of cublasSdot/cublasDdot (and the cublasHandle_t machinery) in this codebase - multiDotGPU no longer needs cuBLAS since its rewrite into a hand-written reduction kernel. Give dotGPU the same treatment: DotKernel is the same single-pass grid-stride + warp-shuffle + block-combine + atomicAdd reduction as MultiDotKernel, specialized for the one-running-sum case (no per-thread array, no dynamic shared memory - just a register accumulator and a small static warpSums buffer, since the launch shape is fixed here). Like multiDotGPU, this runs on the default stream, not a dedicated one: x/y may have just been written by VecExpr::AssignDeviceExpression, which launches on the default stream with no synchronization of its own, so a separate stream would have no guaranteed ordering against that write (see the multiDotGPU stream fix). dotGPU always synchronizes before returning anyway (it needs a real value for the MPI reduction), so there's no async benefit given up. With both GPU reductions now free of cuBLAS, drop the library dependency entirely: removed cublasHandle_t/GetBlasHandle and from CSysVectorGPU.cu, corrected a stale comment in CSysVector.hpp's dot() that attributed a real constraint (the device path needs a materialized vector, not an expression template) to "cuBLAS" specifically, and removed the cublas find_library() call from meson.build - confirmed cuBLAS was not referenced anywhere else in the codebase, and the built SU2_CFD binary no longer links libcublas. Verified: full reconfigure/rebuild succeeds cuBLAS-free. GPU FGCRODR still matches CPU FGCRODR exactly (dotGPU underlies most residual/ convergence-check dot products across every GPU config, not just FGCRODR); Jacobi/Q_Jacobi/ILU reproduce their known values exactly. Reproducible over repeated runs. Co-Authored-By: Claude Sonnet 5 --- Common/include/linear_algebra/CSysVector.hpp | 3 +- Common/src/linear_algebra/CSysVectorGPU.cu | 87 +++++++++++++------- meson.build | 7 +- 3 files changed, 60 insertions(+), 37 deletions(-) diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 87367f5f151..28c75004958 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -528,7 +528,8 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> if constexpr (su2_gpu_capable_v) { using DeviceExpr = std::remove_cv_t>; static_assert(std::is_same_v, - "On the device the dot product is a cuBLAS call, so it only takes vectors. " + "On the device the dot product needs a real device pointer (dotGPU takes a " + "materialized vector, not an expression template), so it only takes vectors. " "Assign the expression to a vector first."); if (VecExpr::UseDeviceExpressions()) { /*--- dotGPU reduces over MPI, which has to happen once for the team, so the result diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 7f669cb31ad..13ea9227b63 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -27,24 +27,52 @@ #include "../../include/linear_algebra/CSysVector.hpp" #include "../../include/linear_algebra/GPUComms.cuh" -#include #include -#include #include namespace { -/*--- cuBLAS handle for the reductions. Created on first use and kept for the lifetime of - * the program, matching the fact that CUDA is either on or off for the whole run. ---*/ -cublasHandle_t solver_blas_handle = nullptr; +/*! + * \brief Fixed launch shape for DotKernel, so the number of warps per block (and thus the static + * shared-memory reduction buffer's size) is a compile-time constant - no dynamic shared + * memory needed for a single running sum, unlike MultiDotKernel's per-(i,j) buffer. + */ +constexpr unsigned int DOT_THREADS_PER_BLOCK = 256; +constexpr unsigned int DOT_WARPS_PER_BLOCK = DOT_THREADS_PER_BLOCK / 32u; -cublasHandle_t GetBlasHandle() { - if (solver_blas_handle == nullptr) { - if (cublasCreate(&solver_blas_handle) != CUBLAS_STATUS_SUCCESS) { - SU2_MPI::Error("cuBLAS handle creation failed for the GPU linear algebra.", CURRENT_FUNCTION); - } +/*! + * \brief Single dot product result[0] = , same single-pass grid-stride + warp-shuffle + + * block-combine + atomicAdd reduction as MultiDotKernel (see its comment), specialized for + * the one-running-sum case: no per-thread array, no dynamic shared memory, just a plain + * register accumulator and a small static warpSums buffer. + */ +template +__global__ void DotKernel(const ScalarType* __restrict__ x, const ScalarType* __restrict__ y, unsigned long size, + ScalarType* __restrict__ result) { + ScalarType sum = 0; + const auto stride = static_cast(blockDim.x) * gridDim.x; + for (auto k = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; k < size; k += stride) { + sum += x[k] * y[k]; + } + + /*--- Warp-level reduction: after this, lane 0 of every warp holds that warp's true sum. ---*/ + for (int offset = 16; offset > 0; offset >>= 1) sum += __shfl_down_sync(0xFFFFFFFFu, sum, offset); + + /*--- Block-level combine: warp leaders stage their partials in static shared memory, thread 0 + * sums them and issues the one atomicAdd this block contributes to the global result + * (pre-zeroed by the caller). ---*/ + __shared__ ScalarType warpSums[DOT_WARPS_PER_BLOCK]; + const unsigned int lane = threadIdx.x % 32u; + const unsigned int warpId = threadIdx.x / 32u; + + if (lane == 0) warpSums[warpId] = sum; + __syncthreads(); + + if (threadIdx.x == 0) { + ScalarType blockSum = 0; + for (unsigned int w = 0; w < DOT_WARPS_PER_BLOCK; ++w) blockSum += warpSums[w]; + atomicAdd(result, blockSum); } - return solver_blas_handle; } /*! @@ -198,27 +226,26 @@ template ScalarType CSysVector::dotGPU(const CSysVector& other) const { SU2_ZONE_SCOPED /*--- Both operands are already on the device, the caller owns the transfers. This - * reduces over MPI, so it must be called by a single thread (see SU2_DEVICE_REGION). ---*/ - cublasHandle_t handle = GetBlasHandle(); - cublasStatus_t status = CUBLAS_STATUS_SUCCESS; + * reduces over MPI, so it must be called by a single thread (see SU2_DEVICE_REGION). + * Deliberately the default stream, not a dedicated one - same reasoning as multiDotGPU: x/y + * may have just been written by VecExpr::AssignDeviceExpression, which launches on the default + * stream with no synchronization of its own, so a separate stream here would have no + * guaranteed ordering against that write. This function always syncs before returning anyway + * (it needs a real value for the MPI reduction), so there is no async benefit being given up. ---*/ + static ScalarType* d_result = nullptr; + if (d_result == nullptr) gpuErrChk(cudaMalloc(&d_result, sizeof(ScalarType))); + + gpuErrChk(cudaMemsetAsync(d_result, 0, sizeof(ScalarType))); + + const auto blocks = static_cast( + std::min((nElmDomain + DOT_THREADS_PER_BLOCK - 1) / DOT_THREADS_PER_BLOCK, 1024)); + DotKernel<<>>(GetDevicePointer(), other.GetDevicePointer(), nElmDomain, + d_result); ScalarType local_dot = ScalarType(0); - - if constexpr (std::is_same_v) { - status = cublasSdot(handle, static_cast(nElmDomain), GetDevicePointer(), 1, other.GetDevicePointer(), 1, - &local_dot); - } else if constexpr (std::is_same_v) { - status = cublasDdot(handle, static_cast(nElmDomain), GetDevicePointer(), 1, other.GetDevicePointer(), 1, - &local_dot); - } else { - SU2_MPI::Error("Unsupported ScalarType in CSysVector::dotGPU.", CURRENT_FUNCTION); - return ScalarType(0); - } - - if (status != CUBLAS_STATUS_SUCCESS) { - SU2_MPI::Error("cuBLAS dot failed in CSysVector::dotGPU.", CURRENT_FUNCTION); - return ScalarType(0); - } + gpuErrChk(cudaMemcpyAsync(&local_dot, d_result, sizeof(ScalarType), cudaMemcpyDeviceToHost)); + gpuErrChk(cudaStreamSynchronize(nullptr)); + gpuErrChk(cudaGetLastError()); ScalarType global_dot = ScalarType(0); const auto mpi_type = (sizeof(ScalarType) < sizeof(double)) ? MPI_FLOAT : MPI_DOUBLE; diff --git a/meson.build b/meson.build index a852cf7f113..78379a45507 100644 --- a/meson.build +++ b/meson.build @@ -28,15 +28,10 @@ if get_option('enable-cuda') # the MPI link flags that -Dcustom-mpi=true relies on mpicxx to provide. add_global_arguments('-ccbin=' + meson.get_compiler('cpp').cmd_array()[0], language : 'cuda') add_global_link_arguments('-ccbin=' + meson.get_compiler('cpp').cmd_array()[0], language : 'cuda') - cuda_deps = [ - meson.get_compiler('cuda').find_library('cublas', required : true), - ] -else - cuda_deps = [] endif su2_cpp_args = [] -su2_deps = [declare_dependency(include_directories: 'externals/CLI11')] + cuda_deps +su2_deps = [declare_dependency(include_directories: 'externals/CLI11')] default_warning_flags = [] if build_machine.system() != 'windows' From 3b8e13c927309cf391db935cafefb4e8349b7df2 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 16 Aug 2026 09:27:53 -0700 Subject: [PATCH 29/30] Clean up dot/multiDot GPU comments: describe current state, not history Several comments around DotKernel/MultiDotKernel/dotGPU/multiDotGPU narrated the design's history instead of just documenting what's there now: comparisons to the removed cublasgemmBatched/cublasdot implementations, "a single symmetric MULTIDOT_MAX_VEC used to" (the cap these replaced), "unlike a separate cudaMemcpy" (the device pointer-array approach these replaced). Trimmed all of that - the current design's own rationale doesn't need the discarded alternative as a reference point. Also de-duplicated: dotGPU and multiDotGPU each carried the full explanation for why they run on the default stream (verbatim near- duplicates); kept it once on dotGPU and had multiDotGPU cross-reference it. Similarly, MULTIDOT_MAX_NM and MULTIDOT_MAX_VEC_LARGE/SMALL each explained the relationship between the two caps; kept that explanation on the LARGE/SMALL comment (where the asymmetric-cap design that motivates it actually lives) and trimmed it off MULTIDOT_MAX_NM's. No functional change. Co-Authored-By: Claude Sonnet 5 --- Common/src/linear_algebra/CSysVectorGPU.cu | 105 +++++++++------------ 1 file changed, 44 insertions(+), 61 deletions(-) diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 13ea9227b63..dd173818dc5 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -76,40 +76,35 @@ __global__ void DotKernel(const ScalarType* __restrict__ x, const ScalarType* __ } /*! - * \brief Cap on n*m (the number of dot products one MultiDotKernel launch computes). Bounds the - * per-thread accumulator array below (necessarily a compile-time size - CUDA has no - * runtime-sized register/local arrays) and the block's shared-memory reduction buffer. - * FGCRODR's Ritz-value path (CSysSolve.cpp) calls this with n = m+1 where m is + * \brief Cap on n*m (the number of dot products one MultiDotKernel launch computes), sizing the + * per-thread accumulator array below and the block's shared-memory reduction buffer - + * both compile-time sized, since CUDA has no runtime-sized register/local arrays. + * FGCRODR's Ritz-value path (CSysSolve.cpp) can reach n = m+1 with m = * LINEAR_SOLVER_RESTART_DEFLATION (user-configurable, default 4 but not uncommon to raise * into the tens), so n*m grows quadratically with that setting - e.g. m=10 already needs - * 110. Sized generously above that; raise it further if a caller ever needs more (a clear - * SU2_MPI::Error fires instead of silently producing a wrong/truncated result). This is a - * different, unrelated constraint from MULTIDOT_MAX_VEC_LARGE/SMALL below - it is a - * genuine per-thread storage cost (every thread pays for the compile-time-sized array - * below regardless of the runtime n*m), unlike those two, which only bound a launch - * parameter's marshaling size, so it must stay independently modest rather than grow to - * cover MULTIDOT_MAX_VEC_LARGE * MULTIDOT_MAX_VEC_SMALL (which would be impractically - * large per thread). + * 110. Sized generously above realistic usage; multiDotGPU raises a clear SU2_MPI::Error + * rather than silently truncating if a caller ever needs more. */ constexpr unsigned int MULTIDOT_MAX_NM = 1024; /*! - * \brief Caps on the individual vector counts n and m, asymmetric because the two call sites in - * this codebase produce genuinely asymmetric shapes: ModGramSchmidt (CSysSolve.cpp) grows - * m up to the Krylov restart length with n=1 (a row vector, not really a rectangle), while - * FGCRODR's Ritz-value path grows n up to the restart length with m bounded by the - * (usually much smaller) LINEAR_SOLVER_RESTART_DEFLATION - a genuine, independently-sized - * rectangle. The host wrapper always puts whichever of n, m is larger into the "large" - * kernel argument (swapping V/W and transposing the result back if needed, see - * multiDotGPU), so one axis only needs to cover the deflation-count case while the other - * covers the restart-length case, instead of both needing the same generous bound (as a - * single symmetric MULTIDOT_MAX_VEC used to). Bounds the fixed-size pointer arrays passed - * into MultiDotKernel by value as ordinary launch parameters, rather than as a device-side - * array of pointers uploaded via a separate cudaMemcpy - CUDA already marshals launch - * parameters for you as part of the launch itself, regardless of host memory pinning, - * which sidesteps needing pinned host buffers here the way HtDTransfer's L/U copy does - * (see aux_stream/htd_event in CSysMatrix, CUDA silently downgrades cudaMemcpyAsync to a - * blocking copy from ordinary pageable host memory). + * \brief Caps on the individual vector counts n and m, bounding the fixed-size pointer arrays + * passed into MultiDotKernel by value as ordinary launch parameters (see + * MultiDotPointers) - CUDA marshals those for you as part of the launch itself, regardless + * of host memory pinning. + * Asymmetric because the two call sites in this codebase produce genuinely asymmetric + * shapes: ModGramSchmidt (CSysSolve.cpp) grows m up to the Krylov restart length with n=1 + * (a row vector, not really a rectangle), while FGCRODR's Ritz-value path grows n up to + * the restart length with m bounded by the usually much smaller + * LINEAR_SOLVER_RESTART_DEFLATION - a genuine, independently-sized rectangle. multiDotGPU + * always puts whichever of n, m is larger into the "large" kernel argument (swapping V/W + * and transposing the result back if needed), so one axis only needs to cover the + * deflation-count case while the other covers the restart-length case. + * \note Unrelated to MULTIDOT_MAX_NM above: n and m can each independently reach this cap without + * n*m approaching MULTIDOT_MAX_VEC_LARGE * MULTIDOT_MAX_VEC_SMALL, since multiDotGPU checks + * n*m against MULTIDOT_MAX_NM separately (that one is a genuine per-thread storage cost + * paid regardless of the runtime n*m; these two only bound a launch parameter's marshaling + * size). */ constexpr unsigned int MULTIDOT_MAX_VEC_LARGE = 256; constexpr unsigned int MULTIDOT_MAX_VEC_SMALL = 64; @@ -125,21 +120,16 @@ struct MultiDotPointers { /*! * \brief Compute the aCount*bCount matrix of dot products D(a,b) = in a single pass - * over the data, replacing a previous cublasgemmBatched implementation that treated - * each pair as its own independent (M=1, N=1, K=size) GEMM. That was a poor fit twice - * over: (a) it re-read every input vector once per pair instead of once total (aCount* - * bCount reads of length size instead of aCount+bCount), and (b) GEMM kernels parallelize - * across output (M,N) tiles and iterate the reduction (K) dimension within one thread - * block's kernel invocation - for an M=N=1 tile with K in the tens of millions, that - * leaves most of the GPU idle, unlike a dedicated reduction primitive (cublasdot, or - * this kernel's grid-stride loop across many blocks) which parallelizes across all of K. + * over the data: every thread reads each of the aCount+bCount vectors once per k and + * forms all aCount*bCount products from that same read, so the vectors are only ever read + * once total (aCount+bCount reads of length size, not aCount*bCount), and the reduction + * stays memory-bound. * \note Each thread accumulates its own private running aCount*bCount sums while striding over k * (kept in thread-local storage, capped at MULTIDOT_MAX_NM - small enough that it should * stay resident in registers or L1 for realistic shapes, cheap either way next to the - * K-length main loop's DRAM traffic), so the single pass over A/B stays memory-bound. Only - * after that loop do threads combine: a warp-shuffle reduction, then one shared-memory - * combine and one atomicAdd per (a,b) per block, so total atomics are - * O(aCount*bCount * numBlocks), not O(aCount*bCount * size). + * K-length main loop's DRAM traffic). Only after that loop do threads combine: a + * warp-shuffle reduction, then one shared-memory combine and one atomicAdd per (a,b) per + * block, so total atomics are O(aCount*bCount * numBlocks), not O(aCount*bCount * size). * \note A is the "large" argument (up to MULTIDOT_MAX_VEC_LARGE), B the "small" one (up to * MULTIDOT_MAX_VEC_SMALL) - the caller (multiDotGPU) is responsible for putting the larger * of its two vector counts into A, and transposing the result back if it had to swap. @@ -225,13 +215,14 @@ void CSysVector::DtHTransfer(bool trigger) const { template ScalarType CSysVector::dotGPU(const CSysVector& other) const { SU2_ZONE_SCOPED - /*--- Both operands are already on the device, the caller owns the transfers. This - * reduces over MPI, so it must be called by a single thread (see SU2_DEVICE_REGION). - * Deliberately the default stream, not a dedicated one - same reasoning as multiDotGPU: x/y - * may have just been written by VecExpr::AssignDeviceExpression, which launches on the default - * stream with no synchronization of its own, so a separate stream here would have no - * guaranteed ordering against that write. This function always syncs before returning anyway - * (it needs a real value for the MPI reduction), so there is no async benefit being given up. ---*/ + /*--- Both operands are already on the device, the caller owns the transfers. This reduces over + * MPI, so it must be called by a single thread (see SU2_DEVICE_REGION). + * \note Runs on the default stream: x/y may have just been written by + * VecExpr::AssignDeviceExpression, which launches on the default stream with no + * synchronization of its own, so a dedicated stream here would have no guaranteed + * ordering against that write - CUDA only orders operations within the same stream. This + * function also always synchronizes before returning (it needs a real value for the MPI + * reduction), so a dedicated stream would not buy any overlap either. ---*/ static ScalarType* d_result = nullptr; if (d_result == nullptr) gpuErrChk(cudaMalloc(&d_result, sizeof(ScalarType))); @@ -255,8 +246,7 @@ ScalarType CSysVector::dotGPU(const CSysVector& other) const { } /*! - * \brief Multi vector dot product, local(i,j) = , via MultiDotKernel (see its - * comment for why that beats the batched-GEMM approach this replaced). + * \brief Multi vector dot product, local(i,j) = , via MultiDotKernel. * \note Whichever of n, m is larger is passed as MultiDotKernel's "A" (large-capacity) argument * and the other as "B" (small-capacity) - see MULTIDOT_MAX_VEC_LARGE/SMALL for why - so if * m > n, V and W (and their counts) are swapped for the call, and the resulting m*n matrix @@ -292,18 +282,11 @@ su2matrix CSysVector::multiDotGPU(const std::vector Date: Sun, 16 Aug 2026 09:54:54 -0700 Subject: [PATCH 30/30] Apply suggestions from code review Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- Common/src/linear_algebra/CSysVectorGPU.cu | 32 ++++------------------ 1 file changed, 5 insertions(+), 27 deletions(-) diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index dd173818dc5..81e094cda94 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -89,22 +89,10 @@ constexpr unsigned int MULTIDOT_MAX_NM = 1024; /*! * \brief Caps on the individual vector counts n and m, bounding the fixed-size pointer arrays - * passed into MultiDotKernel by value as ordinary launch parameters (see - * MultiDotPointers) - CUDA marshals those for you as part of the launch itself, regardless - * of host memory pinning. + * passed into MultiDotKernel by value as ordinary launch parameters. * Asymmetric because the two call sites in this codebase produce genuinely asymmetric - * shapes: ModGramSchmidt (CSysSolve.cpp) grows m up to the Krylov restart length with n=1 - * (a row vector, not really a rectangle), while FGCRODR's Ritz-value path grows n up to - * the restart length with m bounded by the usually much smaller - * LINEAR_SOLVER_RESTART_DEFLATION - a genuine, independently-sized rectangle. multiDotGPU - * always puts whichever of n, m is larger into the "large" kernel argument (swapping V/W - * and transposing the result back if needed), so one axis only needs to cover the - * deflation-count case while the other covers the restart-length case. - * \note Unrelated to MULTIDOT_MAX_NM above: n and m can each independently reach this cap without - * n*m approaching MULTIDOT_MAX_VEC_LARGE * MULTIDOT_MAX_VEC_SMALL, since multiDotGPU checks - * n*m against MULTIDOT_MAX_NM separately (that one is a genuine per-thread storage cost - * paid regardless of the runtime n*m; these two only bound a launch parameter's marshaling - * size). + * shapes. multiDotGPU always puts whichever of n, m is larger into the "large" kernel + * argument (transposing the result back if needed). */ constexpr unsigned int MULTIDOT_MAX_VEC_LARGE = 256; constexpr unsigned int MULTIDOT_MAX_VEC_SMALL = 64; @@ -216,13 +204,7 @@ template ScalarType CSysVector::dotGPU(const CSysVector& other) const { SU2_ZONE_SCOPED /*--- Both operands are already on the device, the caller owns the transfers. This reduces over - * MPI, so it must be called by a single thread (see SU2_DEVICE_REGION). - * \note Runs on the default stream: x/y may have just been written by - * VecExpr::AssignDeviceExpression, which launches on the default stream with no - * synchronization of its own, so a dedicated stream here would have no guaranteed - * ordering against that write - CUDA only orders operations within the same stream. This - * function also always synchronizes before returning (it needs a real value for the MPI - * reduction), so a dedicated stream would not buy any overlap either. ---*/ + * MPI, so it must be called by a single thread (see SU2_DEVICE_REGION). ---*/ static ScalarType* d_result = nullptr; if (d_result == nullptr) gpuErrChk(cudaMalloc(&d_result, sizeof(ScalarType))); @@ -282,11 +264,7 @@ su2matrix CSysVector::multiDotGPU(const std::vector