diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp index 41d3c747cf8..a28c346b651 100644 --- a/Common/include/code_config.hpp +++ b/Common/include/code_config.hpp @@ -50,6 +50,13 @@ #define NEVERINLINE inline #endif +/*--- 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 +#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/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index e23f5de381f..2b6c67701ff 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -107,13 +107,25 @@ 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 even though this preconditioner's own + * operation is a no-op. CSysMatrix::QuantizeDiagonalBlocks() when quantization is off. */ template class CIdentityPreconditioner final : public CPreconditioner { + private: + CSysMatrix& sparse_matrix; + public: + inline explicit CIdentityPreconditioner(CSysMatrix& matrix_ref) : sparse_matrix(matrix_ref) {} + + CIdentityPreconditioner() = delete; + inline void operator()(const CSysVector& u, CSysVector& v) const override { v = u; } inline bool IsIdentity() const override { return true; } + + inline void Build() override { sparse_matrix.QuantizeDiagonalBlocks(); } }; /*! @@ -157,7 +169,9 @@ 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. */ inline void Build() override { sparse_matrix.BuildJacobiPreconditioner(); } }; @@ -248,36 +262,10 @@ 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, no-op for plain LU_SGS. + */ inline void Build() override { sparse_matrix.QuantizeDiagonalBlocks(); } }; @@ -404,19 +392,19 @@ CPreconditioner* CPreconditioner::Create(ENUM_LINEAR_SOL switch (kind) { case IDENTITY: - prec = new CIdentityPreconditioner(); + case Q_IDENTITY: + prec = new CIdentityPreconditioner(jacobian); break; case JACOBI: + case Q_JACOBI: 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/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 38f3c902214..ccd9cae7877 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,48 +110,103 @@ 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 which is not used in this build. + */ +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 Branches on __CUDA_ARCH__, plain memcpy compiles for the device but does not work! */ -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. + * \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 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 -FORCEINLINE 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)); - uint32_t fb; - memcpy(&fb, &fv, sizeof(fb)); - max_abs_bits = std::max(max_abs_bits, fb & 0x7FFFFFFFu); - } - const int e = std::min(127, std::max(-128, static_cast(max_abs_bits >> 23) - 133)); - qs[r] = static_cast(e); - const uint32_t inv_bits = static_cast(127 - e) << 23; - float inv_rscale; - memcpy(&inv_rscale, &inv_bits, sizeof(inv_rscale)); - 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)))); - } +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(passive(r, c)); +#ifdef __CUDA_ARCH__ + const uint32_t fb = __float_as_uint(fv); +#else + 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 = 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[c] = + static_cast(QuantMax(-128.f, QuantMin(127.f, roundf(static_cast(passive(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). + */ +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 @@ -231,14 +287,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,30 +308,37 @@ 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. - * mat.l and mat.u are NOT allocated; off-diagonal blocks live in the - * q_* arrays below. */ + /*! \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 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(). */ - QuantType* q_blocks_d; /*!< \brief Same as q_blocks_l for the diagonal entries. */ + /*!< \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(). .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. */ + LDU d_q_scale; + LDU d_q_blocks; bool useCuda = false; /*!< \brief Whether CUDA is enabled. */ @@ -324,8 +392,18 @@ 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; + /*--- 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. */ @@ -542,7 +620,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; @@ -554,6 +632,12 @@ class CSysMatrix { void MatrixVectorProductGPU(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const; + /*! + * \brief Quantize the diagonal blocks directly on the device. + * \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. @@ -667,10 +751,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; @@ -797,9 +881,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; } @@ -867,9 +951,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]; @@ -910,9 +994,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; } @@ -973,9 +1057,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/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 1b57a6ffe94..28c75004958 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 @@ -530,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 @@ -709,4 +708,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/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/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/CConfig.cpp b/Common/src/CConfig.cpp index 3fe4c98bbe9..c2633f84ec4 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -7481,6 +7481,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: @@ -7490,6 +7491,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 c9668aa0e88..7bee9a09347 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,7 +118,7 @@ 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); @@ -130,15 +126,25 @@ CSysMatrix::~CSysMatrix() { freeHostLDU(mat); freeHostLDU(ilu); 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); + + /*--- 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.l); + GPUMemoryAllocation::pinned_free(m.u); + } else { + MemoryAllocation::aligned_free(m.l); + MemoryAllocation::aligned_free(m.u); + } + }; + freeQuantLDU(q_scale); + freeQuantLDU(q_blocks); 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,13 +155,16 @@ 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); #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 (aux_stream != nullptr) cudaStreamDestroy(aux_stream); + if (htd_event != nullptr) cudaEventDestroy(htd_event); #endif } @@ -207,16 +216,18 @@ 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); + /*--- 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 /*--- 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,21 +253,30 @@ 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 + /*--- .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)); }; - 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); + 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)); + } + }; + 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 { allocAndInit(mat.l, mat.nnz_l * nVar * nEqn); allocAndInit(mat.u, mat.nnz_u * nVar * nEqn); @@ -274,12 +294,28 @@ 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 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); + } else { + GPUAllocAndInit(gpu.l, mat.nnz_l * nVar * nEqn); + GPUAllocAndInit(gpu.u, mat.nnz_u * nVar * nEqn); + } } if (type == ConnectivityType::FiniteVolume) { @@ -350,9 +386,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; } @@ -731,13 +767,31 @@ template void CSysMatrix::QuantizeDiagonalBlocks() { SU2_ZONE_SCOPED - if (quantized_mode) { - /*--- Q_LU_SGS: 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) { + /*--- 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. ---*/ + 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 @@ -758,10 +812,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 } @@ -863,10 +917,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) { @@ -938,6 +992,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) { @@ -1561,9 +1620,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 1cd5de51d5c..470b7a22bcf 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" @@ -33,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 @@ -98,9 +108,54 @@ __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 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. 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; + * 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). */ template __global__ void InvertDiagonalBlocksKernel(unsigned long nRows, unsigned long nVar, @@ -110,15 +165,41 @@ __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); + ParallelMatrixInverse(nVar, i, j, A, Inv); + + invM[iRow * blockSize + tid] = Inv[tid]; +} + +/*! + * \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 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, + const ScalarType* __restrict__ mat_d, int8_t* __restrict__ q_scale_d, + int8_t* __restrict__ q_blocks_d) { + 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 + 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); } /*! @@ -214,11 +295,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]; } /*! @@ -315,8 +402,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 +419,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 */ @@ -353,6 +445,58 @@ __global__ void BlockLDU_SpMV_kernel(unsigned long nRows, unsigned long nVar, y[iRow * nVar + iVar] = sum; } +/*! + * \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( + 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 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 = DecodeQuantScale(qs[iVar]); + const int8_t* __restrict__ row = qv + iVar * nVar; + ScalarType partial = 0; + 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; + }; + + 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 @@ -368,8 +512,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. ---*/ @@ -377,6 +523,25 @@ 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 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. ---*/ + gpuErrChk(cudaStreamSynchronize(nullptr)); + gpuErrChk(cudaGetLastError()); +} + template void CSysMatrix::BuildJacobiPreconditionerGPU() { SU2_ZONE_SCOPED @@ -389,9 +554,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()); @@ -419,7 +583,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 @@ -433,7 +597,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) { @@ -441,17 +605,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()); } @@ -481,7 +645,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, @@ -515,18 +679,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()); } @@ -535,8 +699,28 @@ 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 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 so this transfer can + * run concurrently with whatever the preconditioner's Build() launches next on the default + * stream. ---*/ + 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, + aux_stream)); + gpuErrChk(cudaMemcpyAsync(d_q_blocks.l, q_blocks.l, sizeof(QuantType) * mat.nnz_l * nVar * nEqn, + cudaMemcpyHostToDevice, aux_stream)); + gpuErrChk(cudaMemcpyAsync(d_q_scale.u, q_scale.u, sizeof(QuantType) * mat.nnz_u * nVar, cudaMemcpyHostToDevice, + aux_stream)); + gpuErrChk(cudaMemcpyAsync(d_q_blocks.u, q_blocks.u, sizeof(QuantType) * mat.nnz_u * nVar * nEqn, + 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)); + } } template @@ -547,11 +731,25 @@ 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); - 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); + /*--- 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) { + /*--- 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. ---*/ + 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); + } 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 +761,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, \ diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index dcd88c452c5..3c9176d8716 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -1566,6 +1566,8 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con break; case JACOBI: case LINELET: + case Q_JACOBI: + /*--- BuildJacobiPreconditioner() quantizes the diagonal itself when needed. ---*/ if (RequiresTranspose) Jacobian.BuildJacobiPreconditioner(); break; case LU_SGS: diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index ba53b0cbc29..81e094cda94 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -27,24 +27,146 @@ #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; + +/*! + * \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); + } +} -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 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 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, bounding the fixed-size pointer arrays + * passed into MultiDotKernel by value as ordinary launch parameters. + * Asymmetric because the two call sites in this codebase produce genuinely asymmetric + * 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; + +/*! + * \brief Fixed-size array of device pointers, passed to MultiDotKernel by value (see + * MULTIDOT_MAX_VEC_LARGE/SMALL for why). + */ +template +struct MultiDotPointers { + const ScalarType* ptr[MaxCount]; +}; + +/*! + * \brief Compute the aCount*bCount matrix of dot products D(a,b) = in a single pass + * 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). 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 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 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]; + } + } + + /*--- 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(&D[t], sum); } } - return solver_blas_handle; } } // namespace @@ -81,28 +203,22 @@ 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). ---*/ - cublasHandle_t handle = GetBlasHandle(); - cublasStatus_t status = CUBLAS_STATUS_SUCCESS; + /*--- 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). ---*/ + static ScalarType* d_result = nullptr; + if (d_result == nullptr) gpuErrChk(cudaMalloc(&d_result, sizeof(ScalarType))); - ScalarType local_dot = ScalarType(0); + gpuErrChk(cudaMemsetAsync(d_result, 0, sizeof(ScalarType))); - 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); - } + 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); - if (status != CUBLAS_STATUS_SUCCESS) { - SU2_MPI::Error("cuBLAS dot failed in CSysVector::dotGPU.", CURRENT_FUNCTION); - return ScalarType(0); - } + ScalarType local_dot = 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; @@ -112,102 +228,104 @@ ScalarType CSysVector::dotGPU(const CSysVector& other) const { } /*! - * \brief multi vector product with cublasgemmBatched + * \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 + * 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) { - /*--- 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); + } + + su2matrix 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; - 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. ---*/ + /*--- 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_LARGE/SMALL). ---*/ struct Workspace { - ScalarType* d_local = nullptr; - const ScalarType** d_A = nullptr; - const ScalarType** d_B = nullptr; - ScalarType** d_C = nullptr; + ScalarType* d_D = 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; + void EnsureCapacity(size_t nm) { + if (nm > capacity) { + cudaFree(d_D); + gpuErrChk(cudaMalloc(&d_D, nm * sizeof(ScalarType))); + capacity = nm; + } } - ~Workspace() { - cudaFree(d_local); - cudaFree(d_A); - cudaFree(d_B); - cudaFree(d_C); - } + ~Workspace() { cudaFree(d_D); } }; static Workspace ws; + ws.EnsureCapacity(nm); - // 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 - } - } - - // 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)); + 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 { - SU2_MPI::Error("Unsupported ScalarType in CSysVector::multiDotGPU.", CURRENT_FUNCTION); - return local; + 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(); } - if (status != CUBLAS_STATUS_SUCCESS) { - SU2_MPI::Error("cuBLAS cublasgemmBatched failed in CSysVector::multiDotGPU.", CURRENT_FUNCTION); - return local; + gpuErrChk(cudaMemsetAsync(ws.d_D, 0, nm * sizeof(ScalarType))); + + 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 result to host for MPI reduce - gpuErrChk(cudaMemcpy(local.data(), ws.d_local, batch * sizeof(ScalarType), cudaMemcpyDeviceToHost)); + 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)); + 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)); + 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]; + } + gpuErrChk(cudaGetLastError()); return local; } 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; } } 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. 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'