From 780448de36963a2884cb29158195b7e39db8b06c Mon Sep 17 00:00:00 2001 From: LwhJesse <256257451+LwhJesse@users.noreply.github.com> Date: Mon, 10 Aug 2026 17:34:13 +0800 Subject: [PATCH] Use 64-bit host sparse indices --- Common/include/code_config.hpp | 5 +- Common/include/geometry/CGeometry.hpp | 12 +- .../include/linear_algebra/CPastixWrapper.hpp | 94 +++++--- Common/include/linear_algebra/CSysMatrix.hpp | 33 +-- Common/include/linear_algebra/CSysVector.hpp | 85 +++---- Common/include/toolboxes/graph_toolbox.hpp | 4 +- Common/src/geometry/CGeometry.cpp | 6 +- Common/src/linear_algebra/CPastixWrapper.cpp | 87 +++++--- Common/src/linear_algebra/CSysMatrix.cpp | 210 ++++++++++++------ Common/src/linear_algebra/CSysMatrixGPU.cu | 68 ++++-- Common/src/linear_algebra/CSysVector.cpp | 47 +++- Common/src/linear_algebra/CSysVectorGPU.cu | 39 +++- SU2_CFD/include/solvers/CFEASolver.hpp | 4 +- .../include/solvers/CFVMFlowSolverBase.hpp | 4 +- .../include/solvers/CFVMFlowSolverBase.inl | 2 +- SU2_CFD/include/solvers/CScalarSolver.hpp | 4 +- SU2_CFD/include/solvers/CScalarSolver.inl | 2 +- SU2_CFD/src/solvers/CFEASolver.cpp | 2 +- 18 files changed, 467 insertions(+), 241 deletions(-) diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp index 41d3c747cf8..f2cf6bc05f6 100644 --- a/Common/include/code_config.hpp +++ b/Common/include/code_config.hpp @@ -111,7 +111,10 @@ FORCEINLINE Out su2staticcast_p(In ptr) { #undef USE_SINGLE_PRECISION #endif -/*--- Default integer types. Currently used for rank-local sparse patterns. ---*/ +/*--- Default integer types. Host sparse metadata is explicit 64-bit so rank-local + * sparse patterns and allocation-size arithmetic do not depend on the platform + * width of unsigned long. Device sparse metadata keeps the compact su2uint type. ---*/ +using su2_index_t = uint64_t; using su2uint = uint32_t; using su2int = int32_t; diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 97d3f845537..5702ba5d62c 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -205,10 +205,10 @@ class CGeometry { LDUSparsePattern finiteElementPatternFill0; /*!< \brief FEM sparsity with 0-fill (structural pattern). */ LDUSparsePattern finiteElementPatternFillN; /*!< \brief FEM sparsity with N-fill (e.g. for ILU-N). */ - su2vector finiteVolumeLToUTranspMap; /*!< \brief FVM L-entry -> U-entry of its transpose. */ - su2vector finiteVolumeUToLTranspMap; /*!< \brief FVM U-entry -> L-entry of its transpose. */ - su2vector finiteElementLToUTranspMap; /*!< \brief FEM L-entry -> U-entry of its transpose. */ - su2vector finiteElementUToLTranspMap; /*!< \brief FEM U-entry -> L-entry of its transpose. */ + su2vector finiteVolumeLToUTranspMap; /*!< \brief FVM L-entry -> U-entry of its transpose. */ + su2vector finiteVolumeUToLTranspMap; /*!< \brief FVM U-entry -> L-entry of its transpose. */ + su2vector finiteElementLToUTranspMap; /*!< \brief FEM L-entry -> U-entry of its transpose. */ + su2vector finiteElementUToLTranspMap; /*!< \brief FEM U-entry -> L-entry of its transpose. */ /*--- Edge and element colorings. ---*/ @@ -1892,7 +1892,7 @@ class CGeometry { * \param[in] type - Finite volume or finite element. * \return Reference to the l_to_u map. */ - const su2vector& GetLToUTransposeSparsePatternMap(ConnectivityType type); + const su2vector& GetLToUTransposeSparsePatternMap(ConnectivityType type); /*! * \brief Get the bijective map from U-entry indices to L-entry indices of their transposes. @@ -1900,7 +1900,7 @@ class CGeometry { * \param[in] type - Finite volume or finite element. * \return Reference to the u_to_l map. */ - const su2vector& GetUToLTransposeSparsePatternMap(ConnectivityType type); + const su2vector& GetUToLTransposeSparsePatternMap(ConnectivityType type); /*! * \brief Get the edge coloring. diff --git a/Common/include/linear_algebra/CPastixWrapper.hpp b/Common/include/linear_algebra/CPastixWrapper.hpp index 44fb9d503db..a8646f6e6e2 100644 --- a/Common/include/linear_algebra/CPastixWrapper.hpp +++ b/Common/include/linear_algebra/CPastixWrapper.hpp @@ -38,6 +38,8 @@ #include #include +#include +#include #include using namespace std; @@ -64,25 +66,30 @@ class CPastixWrapper { vector perm; /*!< \brief Ordering computed by PaStiX. */ vector workvec; /*!< \brief RHS vector which then becomes the solution. */ - vector csr_row_ptr; /*!< \brief Owned CSR row pointers (built from LDU). */ - vector csr_col_ind; /*!< \brief Owned CSR column indices (built from LDU). */ + vector csr_row_ptr; /*!< \brief Owned CSR row pointers (built from LDU). */ + vector csr_col_ind; /*!< \brief Owned CSR column indices (built from LDU). */ pastix_int_t iparm[IPARM_SIZE]; /*!< \brief Integer parameters for PaStiX. */ double dparm[DPARM_SIZE]; /*!< \brief Floating point parameters for PaStiX. */ struct { - unsigned long nVar = 0; - unsigned long nPoint = 0; - unsigned long nPointDomain = 0; - unsigned long blkSz = 0; /*!< \brief Block size (nVar * nVar) for value assembly. */ - - const su2uint* row_ptr_l = nullptr; /*!< \brief LDU lower row pointers (geometry-owned). */ - const su2uint* row_ptr_u = nullptr; /*!< \brief LDU upper row pointers (geometry-owned). */ - const ScalarType* d = nullptr; /*!< \brief Diagonal blocks (matrix-owned). */ - const ScalarType* l = nullptr; /*!< \brief Lower blocks (matrix-owned). */ - const ScalarType* u = nullptr; /*!< \brief Upper blocks (matrix-owned). */ - - unsigned long size_rhs() const { return nPointDomain * nVar; } + su2_index_t nVar = 0; + su2_index_t nPoint = 0; + su2_index_t nPointDomain = 0; + su2_index_t blkSz = 0; /*!< \brief Block size (nVar * nVar) for value assembly. */ + + const su2_index_t* row_ptr_l = nullptr; /*!< \brief LDU lower row pointers (geometry-owned). */ + const su2_index_t* row_ptr_u = nullptr; /*!< \brief LDU upper row pointers (geometry-owned). */ + const ScalarType* d = nullptr; /*!< \brief Diagonal blocks (matrix-owned). */ + const ScalarType* l = nullptr; /*!< \brief Lower blocks (matrix-owned). */ + const ScalarType* u = nullptr; /*!< \brief Upper blocks (matrix-owned). */ + + su2_index_t size_rhs() const { + if (nPointDomain != 0 && nVar > std::numeric_limits::max() / nPointDomain) { + SU2_MPI::Error("Overflow while computing PaStiX rhs size.", CURRENT_FUNCTION); + } + return nPointDomain * nVar; + } } matrix; /*!< \brief Dimensions and LDU pointers captured from the owning CSysMatrix. */ bool issetup{}; /*!< \brief Signals that the structure has been provided. */ @@ -94,8 +101,30 @@ class CPastixWrapper { const int mpi_size = SU2_MPI::GetSize(); const int mpi_rank = SU2_MPI::GetRank(); - vector sort_rows; /*!< \brief List of rows with halo points. */ - vector> sort_order; /*!< \brief How each of those rows needs to be sorted. */ + vector sort_rows; /*!< \brief List of rows with halo points. */ + vector> sort_order; /*!< \brief How each of those rows needs to be sorted. */ + + static su2_index_t CheckedMul(su2_index_t lhs, su2_index_t rhs, const char* what) { + if (lhs != 0 && rhs > std::numeric_limits::max() / lhs) { + SU2_MPI::Error(std::string("Overflow while computing ") + what + ".", CURRENT_FUNCTION); + } + return lhs * rhs; + } + + static su2_index_t CheckedAdd(su2_index_t lhs, su2_index_t rhs, const char* what) { + if (rhs > std::numeric_limits::max() - lhs) { + SU2_MPI::Error(std::string("Overflow while computing ") + what + ".", CURRENT_FUNCTION); + } + return lhs + rhs; + } + + template + static TargetType CheckedCast(su2_index_t value, const char* what) { + if (value > static_cast(std::numeric_limits::max())) { + SU2_MPI::Error(std::string("Overflow while converting ") + what + ".", CURRENT_FUNCTION); + } + return static_cast(value); + } /*! * \brief Run the "clean" task, releases all memory, leaves object in unusable state. @@ -152,9 +181,9 @@ class CPastixWrapper { * \param[in] col_ind_l/u - LDU lower/upper column indices (geometry-owned). * \param[in] d/l/u - LDU value blocks (matrix-owned, must outlive wrapper). */ - void SetLDU(unsigned long nVar, unsigned long nPoint, unsigned long nPointDomain, const su2uint* row_ptr_l, - const su2uint* col_ind_l, const su2uint* row_ptr_u, const su2uint* col_ind_u, const ScalarType* d, - const ScalarType* l, const ScalarType* u) { + void SetLDU(su2_index_t nVar, su2_index_t nPoint, su2_index_t nPointDomain, const su2_index_t* row_ptr_l, + const su2_index_t* col_ind_l, const su2_index_t* row_ptr_u, const su2_index_t* col_ind_u, + const ScalarType* d, const ScalarType* l, const ScalarType* u) { if (issetup) return; matrix.nVar = nVar; matrix.nPoint = nPoint; @@ -164,18 +193,20 @@ class CPastixWrapper { matrix.d = d; matrix.l = l; matrix.u = u; - matrix.blkSz = nVar * nVar; - - const unsigned long nnz_domain = row_ptr_l[nPointDomain] + nPointDomain + row_ptr_u[nPointDomain]; - csr_row_ptr.resize(nPointDomain + 1); - csr_col_ind.reserve(nnz_domain); - for (auto i = 0ul; i < nPointDomain; ++i) { - csr_row_ptr[i] = static_cast(csr_col_ind.size()); + matrix.blkSz = CheckedMul(nVar, nVar, "PaStiX block size"); + + const auto nnz_domain = CheckedAdd(CheckedAdd(row_ptr_l[nPointDomain], nPointDomain, "PaStiX CSR index size"), + row_ptr_u[nPointDomain], "PaStiX CSR index size"); + csr_row_ptr.resize( + CheckedCast(CheckedAdd(nPointDomain, 1, "PaStiX CSR row pointer size"), "PaStiX CSR row pointer size")); + csr_col_ind.reserve(CheckedCast(nnz_domain, "PaStiX CSR column index size")); + for (su2_index_t i = 0; i < nPointDomain; ++i) { + csr_row_ptr[i] = static_cast(csr_col_ind.size()); for (auto k = row_ptr_l[i]; k < row_ptr_l[i + 1]; ++k) csr_col_ind.push_back(col_ind_l[k]); csr_col_ind.push_back(i); for (auto k = row_ptr_u[i]; k < row_ptr_u[i + 1]; ++k) csr_col_ind.push_back(col_ind_u[k]); } - csr_row_ptr[nPointDomain] = static_cast(csr_col_ind.size()); + csr_row_ptr[nPointDomain] = static_cast(csr_col_ind.size()); issetup = true; } @@ -210,11 +241,14 @@ class CPastixWrapper { } iparm[IPARM_VERBOSE] = PastixVerboseNot; - for (auto i = 0ul; i < matrix.size_rhs(); ++i) workvec[i] = rhs[i]; - if (pastix_task_solve(state, matrix.size_rhs(), 1, workvec.data(), matrix.size_rhs()) != PASTIX_SUCCESS) { + const auto rhs_size = matrix.size_rhs(); + const auto pastix_rhs_size = CheckedCast(rhs_size, "PaStiX rhs size"); + + for (su2_index_t i = 0; i < rhs_size; ++i) workvec[i] = rhs[i]; + if (pastix_task_solve(state, pastix_rhs_size, 1, workvec.data(), pastix_rhs_size) != PASTIX_SUCCESS) { SU2_MPI::Error("Error solving linear system.", CURRENT_FUNCTION); } - for (auto i = 0ul; i < matrix.size_rhs(); ++i) sol[i] = workvec[i]; + for (su2_index_t i = 0; i < rhs_size; ++i) sol[i] = workvec[i]; } }; #endif diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 38f3c902214..547797ae9b6 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -235,22 +235,23 @@ class CSysMatrix { * 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. */ + 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. */ - 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. */ - const su2uint* col_ind_u = nullptr; /*!< \brief Column indices for U. */ - unsigned long nnz_l = 0; /*!< \brief Number of L nonzeros. */ - unsigned long nnz_u = 0; /*!< \brief Number of U nonzeros. */ + ScalarType* d = nullptr; /*!< \brief Diagonal block values. */ + ScalarType* l = nullptr; /*!< \brief Strictly-lower block values. */ + ScalarType* u = nullptr; /*!< \brief Strictly-upper block values. */ + const IndexType* row_ptr_l = nullptr; /*!< \brief Row pointers for L (geometry-owned or GPU copy). */ + const IndexType* col_ind_l = nullptr; /*!< \brief Column indices for L. */ + const IndexType* row_ptr_u = nullptr; /*!< \brief Row pointers for U. */ + const IndexType* col_ind_u = nullptr; /*!< \brief Column indices for U. */ + su2_index_t nnz_l = 0; /*!< \brief Number of L nonzeros. */ + su2_index_t 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). ---*/ @@ -278,15 +279,15 @@ class CSysMatrix { * Linelet preconditioner, which builds the Jacobi one but reads invM on the host. */ bool jacobi_on_device = false; - const su2uint* l_to_u_transp; /*!< \brief L-entry index -> U-entry index of its transpose. */ - const su2uint* u_to_l_transp; /*!< \brief U-entry index -> L-entry index of its transpose. */ + const su2_index_t* l_to_u_transp; /*!< \brief L-entry index -> U-entry index of its transpose. */ + const su2_index_t* u_to_l_transp; /*!< \brief U-entry index -> L-entry index of its transpose. */ /*! * \brief Lookup table from edges to the L-index in the LDU split. * U-index == edge index by construction (edges are ordered 1:1 with the U pattern). * Therefore, edge_ptr_l == u_to_l_transp, but we keep a separate member for clarity. */ - const su2uint* edge_ptr_l; + const su2_index_t* edge_ptr_l; unsigned short ilu_fill_in; /*!< \brief Fill level for the ILU preconditioner. */ diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 2ffbb9a7967..a81da28cd95 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -28,7 +28,9 @@ #pragma once +#include #include +#include #include #include "../parallelization/mpi_structure.hpp" @@ -155,12 +157,12 @@ struct store_type> { }; template -void AssignDeviceExpression(Scalar* data, unsigned long size, const CVecExpr& expr); +void AssignDeviceExpression(Scalar* data, su2_index_t size, const CVecExpr& expr); #ifdef __CUDACC__ template -__global__ void DeviceAssignKernel(Scalar* data, unsigned long size, T expr) { - const unsigned long i = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; +__global__ void DeviceAssignKernel(Scalar* data, su2_index_t size, T expr) { + const su2_index_t i = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; if (i >= size) return; if constexpr (Op == DeviceAssignOp::Assign) { @@ -177,10 +179,14 @@ __global__ void DeviceAssignKernel(Scalar* data, unsigned long size, T expr) { } template -inline void AssignDeviceExpression(Scalar* data, unsigned long size, const CVecExpr& expr) { +inline void AssignDeviceExpression(Scalar* data, su2_index_t size, const CVecExpr& expr) { if (size == 0) return; constexpr unsigned block_size = 256; - const auto grid_size = static_cast((size + block_size - 1) / block_size); + const auto blocks = roundUpDiv(size, block_size); + if (blocks > static_cast(std::numeric_limits::max())) { + SU2_MPI::Error("CSysVector device assignment exceeds the CUDA grid dimension range.", CURRENT_FUNCTION); + } + const auto grid_size = static_cast(blocks); DeviceAssignKernel<<>>(data, size, expr.derived()); gpuErrChk(cudaPeekAtLastError()); } @@ -201,9 +207,9 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> /// NOTE: Update swap() if you add member variables. unsigned long omp_chunk_size = OMP_MAX_SIZE; /*!< \brief Static chunk size used in loops. */ ScalarType* vec_val = nullptr; /*!< \brief Storage, 64 byte aligned (do not use normal new/delete). */ - unsigned long nElm = 0; /*!< \brief Total number of elements (or number elements on this processor). */ - unsigned long nElmDomain = 0; /*!< \brief Total number of elements without Ghost cells. */ - unsigned long nVar = 1; /*!< \brief Number of elements in a block. */ + su2_index_t nElm = 0; /*!< \brief Total number of elements (or number elements on this processor). */ + su2_index_t nElmDomain = 0; /*!< \brief Total number of elements without Ghost cells. */ + su2_index_t nVar = 1; /*!< \brief Number of elements in a block. */ ScalarType* d_vec_val = nullptr; /*!< \brief Device Pointer to store the vector values on the GPU. */ @@ -224,7 +230,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] valIsArray - If true val is treated as array. * \param[in] errorIfParallel - Throw error if within parallel region (all ctors except the default one do this). */ - void Initialize(unsigned long numBlk, unsigned long numBlkDomain, unsigned long numVar, const ScalarType* val, + void Initialize(su2_index_t numBlk, su2_index_t numBlkDomain, su2_index_t numVar, const ScalarType* val, bool valIsArray, bool errorIfParallel = true); /*! @@ -288,7 +294,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] size - Number of elements locally. * \param[in] val - Default value for elements. */ - explicit CSysVector(unsigned long size, ScalarType val = 0.0) { Initialize(size, size, 1, &val, false); } + explicit CSysVector(su2_index_t size, ScalarType val = 0.0) { Initialize(size, size, 1, &val, false); } /*! * \brief Construct from size and value (block version). @@ -297,7 +303,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] numVar - Number of variables in each block. * \param[in] val - Default value for elements. */ - CSysVector(unsigned long numBlk, unsigned long numBlkDomain, unsigned long numVar, ScalarType val = 0.0) { + CSysVector(su2_index_t numBlk, su2_index_t numBlkDomain, su2_index_t numVar, ScalarType val = 0.0) { Initialize(numBlk, numBlkDomain, numVar, &val, false); } @@ -306,7 +312,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] size - Number of elements locally. * \param[in] u_array - Vector stored as array being copied. */ - CSysVector(unsigned long size, const ScalarType* u_array) { Initialize(size, size, 1, u_array, true); } + CSysVector(su2_index_t size, const ScalarType* u_array) { Initialize(size, size, 1, u_array, true); } /*! * \brief Constructor from array (block version). @@ -315,7 +321,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] numVar - number of variables in each block * \param[in] u_array - vector stored as array being copied */ - CSysVector(unsigned long numBlk, unsigned long numBlkDomain, unsigned long numVar, const ScalarType* u_array) { + CSysVector(su2_index_t numBlk, su2_index_t numBlkDomain, su2_index_t numVar, const ScalarType* u_array) { Initialize(numBlk, numBlkDomain, numVar, u_array, true); } @@ -346,7 +352,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] numVar - number of variables in each block * \param[in] val - default value for elements */ - void Initialize(unsigned long numBlk, unsigned long numBlkDomain, unsigned long numVar, ScalarType val = 0.0) { + void Initialize(su2_index_t numBlk, su2_index_t numBlkDomain, su2_index_t numVar, ScalarType val = 0.0) { Initialize(numBlk, numBlkDomain, numVar, &val, false, false); } @@ -358,7 +364,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] numVar - number of variables in each block * \param[in] ptr - pointer to data with which to initialize the vector */ - void Initialize(unsigned long numBlk, unsigned long numBlkDomain, unsigned long numVar, const ScalarType* ptr) { + void Initialize(su2_index_t numBlk, su2_index_t numBlkDomain, su2_index_t numVar, const ScalarType* ptr) { Initialize(numBlk, numBlkDomain, numVar, ptr, true, false); } @@ -378,7 +384,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> Initialize(other.GetNBlk(), other.GetNBlkDomain(), other.GetNVar(), nullptr, true, false);) CSYSVEC_PARFOR - for (auto i = 0ul; i < nElm; i++) vec_val[i] = SU2_TYPE::GetValue(other[i]); + for (su2_index_t i = 0; i < nElm; i++) vec_val[i] = SU2_TYPE::GetValue(other[i]); END_CSYSVEC_PARFOR } @@ -423,35 +429,35 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> /*! * \brief return the number of local elements in the CSysVector */ - inline unsigned long GetLocSize() const { return nElm; } + inline su2_index_t GetLocSize() const { return nElm; } /*! * \brief return the number of local elements in the CSysVector without ghost cells */ - inline unsigned long GetNElmDomain() const { return nElmDomain; } + inline su2_index_t GetNElmDomain() const { return nElmDomain; } /*! * \brief return the number of variables at each block (typically number per node) */ - inline unsigned long GetNVar() const { return nVar; } + inline su2_index_t GetNVar() const { return nVar; } /*! * \brief return the number of blocks (typically number of nodes locally) */ - inline unsigned long GetNBlk() const { return nElm / nVar; } + inline su2_index_t GetNBlk() const { return nElm / nVar; } /*! * \brief return the number of blocks (typically number of nodes locally) */ - inline unsigned long GetNBlkDomain() const { return nElmDomain / nVar; } + inline su2_index_t GetNBlkDomain() const { return nElmDomain / nVar; } /*! * \brief Access operator with assignment permitted. * \param[in] i - Local index to access. * \return Value at position i. */ - inline ScalarType& operator[](unsigned long i) { return vec_val[i]; } - inline const ScalarType& operator[](unsigned long i) const { return vec_val[i]; } + inline ScalarType& operator[](su2_index_t i) { return vec_val[i]; } + inline const ScalarType& operator[](su2_index_t i) const { return vec_val[i]; } /*! * \brief Iterators for range for loops. @@ -465,8 +471,8 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] iVar - Index of variable. * \return Value at position (i,j). */ - inline ScalarType& operator()(unsigned long iPoint, unsigned long iVar) { return vec_val[iPoint * nVar + iVar]; } - inline const ScalarType& operator()(unsigned long iPoint, unsigned long iVar) const { + inline ScalarType& operator()(su2_index_t iPoint, su2_index_t iVar) { return vec_val[iPoint * nVar + iVar]; } + inline const ScalarType& operator()(su2_index_t iPoint, su2_index_t iVar) const { return vec_val[iPoint * nVar + iVar]; } @@ -482,7 +488,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> } #endif CSYSVEC_PARFOR - for (auto i = 0ul; i < nElm; ++i) vec_val[i] = other.vec_val[i]; + for (su2_index_t i = 0; i < nElm; ++i) vec_val[i] = other.vec_val[i]; END_CSYSVEC_PARFOR return *this; } @@ -497,7 +503,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> if (VecExpr::UseDeviceExpressions()) return AssignDevice(val); \ } \ CSYSVEC_PARFOR \ - for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val; \ + for (su2_index_t i = 0; i < nElm; ++i) vec_val[i] OP val; \ END_CSYSVEC_PARFOR \ return *this; \ } \ @@ -507,7 +513,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> if (VecExpr::UseDeviceExpressions()) return AssignDevice(expr); \ } \ CSYSVEC_PARFOR \ - for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP expr.derived()[i]; \ + for (su2_index_t i = 0; i < nElm; ++i) vec_val[i] OP expr.derived()[i]; \ END_CSYSVEC_PARFOR \ return *this; \ } @@ -552,7 +558,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> ScalarType sum = 0.0; CSYSVEC_PARFOR - for (auto i = 0ul; i < nElmDomain; ++i) { + for (su2_index_t i = 0; i < nElmDomain; ++i) { sum += vec_val[i] * expr.derived()[i]; } END_CSYSVEC_PARFOR @@ -602,15 +608,15 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] iPoint - Index of block. * \return Pointer to start of block. */ - inline ScalarType* GetBlock(unsigned long iPoint) { return &vec_val[iPoint * nVar]; } - inline const ScalarType* GetBlock(unsigned long iPoint) const { return &vec_val[iPoint * nVar]; } + inline ScalarType* GetBlock(su2_index_t iPoint) { return &vec_val[iPoint * nVar]; } + inline const ScalarType* GetBlock(su2_index_t iPoint) const { return &vec_val[iPoint * nVar]; } /*! * \brief Set the values to zero for one block. * \param[in] iPoint - Index of the block being set to zero. */ - inline void SetBlock_Zero(unsigned long iPoint) { - for (auto iVar = 0ul; iVar < nVar; iVar++) vec_val[iPoint * nVar + iVar] = 0.0; + inline void SetBlock_Zero(su2_index_t iPoint) { + for (su2_index_t iVar = 0; iVar < nVar; iVar++) vec_val[iPoint * nVar + iVar] = 0.0; } /*! @@ -621,11 +627,11 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] alpha - Scale factor (axpy-type operation). */ template - FORCEINLINE void SetBlock(unsigned long iPoint, const VectorType& block, ScalarType alpha = 1) { + FORCEINLINE void SetBlock(su2_index_t iPoint, const VectorType& block, ScalarType alpha = 1) { if (Overwrite) { - for (auto i = 0ul; i < nVar; ++i) vec_val[iPoint * nVar + i] = alpha * block[i]; + for (su2_index_t i = 0; i < nVar; ++i) vec_val[iPoint * nVar + i] = alpha * block[i]; } else { - for (auto i = 0ul; i < nVar; ++i) vec_val[iPoint * nVar + i] += alpha * block[i]; + for (su2_index_t i = 0; i < nVar; ++i) vec_val[iPoint * nVar + i] += alpha * block[i]; } } @@ -633,7 +639,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \brief Add "block" to the vector, see SetBlock. */ template - FORCEINLINE void AddBlock(unsigned long iPoint, const VectorType& block, ScalarType alpha = 1) { + FORCEINLINE void AddBlock(su2_index_t iPoint, const VectorType& block, ScalarType alpha = 1) { SetBlock(iPoint, block, alpha); } @@ -641,7 +647,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \brief Subtract "block" from the vector, see AddBlock. */ template - FORCEINLINE void SubtractBlock(unsigned long iPoint, const VectorType& block) { + FORCEINLINE void SubtractBlock(su2_index_t iPoint, const VectorType& block) { AddBlock(iPoint, block, -1); } @@ -649,8 +655,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \brief Add to iPoint, subtract from jPoint. */ template - FORCEINLINE void UpdateBlocks(unsigned long iPoint, unsigned long jPoint, const VectorType& block, - ScalarType alpha = 1) { + FORCEINLINE void UpdateBlocks(su2_index_t iPoint, su2_index_t jPoint, const VectorType& block, ScalarType alpha = 1) { AddBlock(iPoint, block, alpha); AddBlock(jPoint, block, -alpha); } diff --git a/Common/include/toolboxes/graph_toolbox.hpp b/Common/include/toolboxes/graph_toolbox.hpp index f1c36969ada..8634f77d058 100644 --- a/Common/include/toolboxes/graph_toolbox.hpp +++ b/Common/include/toolboxes/graph_toolbox.hpp @@ -334,9 +334,9 @@ class CCompressedSparsePattern { template using CEdgeToNonZeroMap = C2DContainer; -using CCompressedSparsePatternUL = CCompressedSparsePattern; +using CCompressedSparsePatternUL = CCompressedSparsePattern; using CCompressedSparsePatternL = CCompressedSparsePattern; -using CEdgeToNonZeroMapUL = CEdgeToNonZeroMap; +using CEdgeToNonZeroMapUL = CEdgeToNonZeroMap; /*! * \brief Build a sparse pattern from geometry information, of type FVM or FEM, diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 45153456086..39507a3603d 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -4128,7 +4128,7 @@ const CGeometry::LDUSparsePattern& CGeometry::GetSparsePattern(ConnectivityType auto& grp = fillLvl == 0 ? (fvm ? finiteVolumePatternFill0 : finiteElementPatternFill0) : (fvm ? finiteVolumePatternFillN : finiteElementPatternFillN); if (grp.empty()) { - grp.csr = buildCSRPattern(*this, type, static_cast(fillLvl)); + grp.csr = buildCSRPattern(*this, type, static_cast(fillLvl)); grp.csr.buildDiagPtr(); grp.l = buildLowerPattern(grp.csr); grp.u = buildUpperPattern(grp.csr); @@ -4136,7 +4136,7 @@ const CGeometry::LDUSparsePattern& CGeometry::GetSparsePattern(ConnectivityType return grp; } -const su2vector& CGeometry::GetLToUTransposeSparsePatternMap(ConnectivityType type) { +const su2vector& CGeometry::GetLToUTransposeSparsePatternMap(ConnectivityType type) { bool fvm = (type == ConnectivityType::FiniteVolume); auto& l_to_u = fvm ? finiteVolumeLToUTranspMap : finiteElementLToUTranspMap; if (l_to_u.empty()) { @@ -4147,7 +4147,7 @@ const su2vector& CGeometry::GetLToUTransposeSparsePatternMap(Connectivi return l_to_u; } -const su2vector& CGeometry::GetUToLTransposeSparsePatternMap(ConnectivityType type) { +const su2vector& CGeometry::GetUToLTransposeSparsePatternMap(ConnectivityType type) { bool fvm = (type == ConnectivityType::FiniteVolume); auto& u_to_l = fvm ? finiteVolumeUToLTranspMap : finiteElementUToLTranspMap; if (u_to_l.empty()) { diff --git a/Common/src/linear_algebra/CPastixWrapper.cpp b/Common/src/linear_algebra/CPastixWrapper.cpp index f213939eaf6..a265137acb5 100644 --- a/Common/src/linear_algebra/CPastixWrapper.cpp +++ b/Common/src/linear_algebra/CPastixWrapper.cpp @@ -34,26 +34,29 @@ #include "../../include/geometry/CGeometry.hpp" #include "../../include/linear_algebra/CPastixWrapper.hpp" +#include #include template void CPastixWrapper::Initialize(CGeometry* geometry, const CConfig* config) { if (isinitialized) return; // only need to do this once - const unsigned long nVar = matrix.nVar, nPoint = matrix.nPoint, nPointDomain = matrix.nPointDomain; - const unsigned long *row_ptr = csr_row_ptr.data(), *col_ind = csr_col_ind.data(); - const unsigned long nNonZero = row_ptr[nPointDomain]; + const auto nVar = matrix.nVar, nPoint = matrix.nPoint, nPointDomain = matrix.nPointDomain; + const auto *row_ptr = csr_row_ptr.data(), *col_ind = csr_col_ind.data(); + const auto nNonZero = row_ptr[nPointDomain]; /*--- Allocate ---*/ - nCols = static_cast(nPointDomain); - colptr.resize(nPointDomain + 1); + nCols = CheckedCast(nPointDomain, "PaStiX local column count"); + colptr.resize(CheckedCast(nPointDomain + 1, "PaStiX colptr size")); rowidx.clear(); - rowidx.reserve(nNonZero); - values.resize(nNonZero * nVar * nVar); - loc2glb.resize(nPointDomain); - perm.resize(nPointDomain); - workvec.resize(nPointDomain * nVar); + rowidx.reserve(CheckedCast(nNonZero, "PaStiX rowidx size")); + values.resize(CheckedCast(CheckedMul(nNonZero, matrix.blkSz, "PaStiX value storage size"), + "PaStiX value storage size")); + loc2glb.resize(CheckedCast(nPointDomain, "PaStiX loc2glb size")); + perm.resize(CheckedCast(nPointDomain, "PaStiX permutation size")); + workvec.resize( + CheckedCast(CheckedMul(nPointDomain, nVar, "PaStiX work vector size"), "PaStiX work vector size")); /*--- Set default parameter values ---*/ @@ -90,19 +93,22 @@ void CPastixWrapper::Initialize(CGeometry* geometry, const CConfig* /*--- 1 - Determine position in the linear partitioning ---*/ - unsigned long offset = 0; + su2_index_t offset = 0; #ifdef HAVE_MPI - vector domain_sizes(mpi_size); - MPI_Allgather(&nPointDomain, 1, MPI_UNSIGNED_LONG, domain_sizes.data(), 1, MPI_UNSIGNED_LONG, SU2_MPI::GetComm()); - for (int i = 0; i < mpi_rank; ++i) offset += domain_sizes[i]; + vector domain_sizes(mpi_size); + MPI_Allgather(&nPointDomain, 1, MPI_UINT64_T, domain_sizes.data(), 1, MPI_UINT64_T, SU2_MPI::GetComm()); + for (int i = 0; i < mpi_rank; ++i) offset = CheckedAdd(offset, domain_sizes[i], "PaStiX global offset"); #endif - iota(loc2glb.begin(), loc2glb.end(), offset + 1); + for (su2_index_t i = 0; i < nPointDomain; ++i) { + loc2glb[static_cast(i)] = CheckedCast(CheckedAdd(offset, i + 1, "PaStiX global mapped index"), + "PaStiX global mapped index"); + } /*--- 2 - Communicate global indices of halo points to then renumber column indices from local to global when unpacking halos. ---*/ - vector map(nPoint - nPointDomain, 0); + vector map(CheckedCast(nPoint - nPointDomain, "PaStiX halo map size"), 0); #ifdef HAVE_MPI for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { @@ -116,15 +122,16 @@ void CPastixWrapper::Initialize(CGeometry* geometry, const CConfig* unsigned long nVertexR = geometry->nVertex[MarkerR]; /*--- Allocate Send/Receive buffers ---*/ - vector Buffer_Recv(nVertexR), Buffer_Send(nVertexS); + vector Buffer_Recv(nVertexR), Buffer_Send(nVertexS); /*--- Prepare data to send ---*/ for (unsigned long iVertex = 0; iVertex < nVertexS; iVertex++) - Buffer_Send[iVertex] = geometry->vertex[MarkerS][iVertex]->GetNode() + offset; + Buffer_Send[iVertex] = + CheckedAdd(geometry->vertex[MarkerS][iVertex]->GetNode(), offset, "PaStiX halo send index"); /*--- Send and Receive data ---*/ - MPI_Sendrecv(Buffer_Send.data(), nVertexS, MPI_UNSIGNED_LONG, sender, 0, Buffer_Recv.data(), nVertexR, - MPI_UNSIGNED_LONG, recver, 0, SU2_MPI::GetComm(), MPI_STATUS_IGNORE); + MPI_Sendrecv(Buffer_Send.data(), nVertexS, MPI_UINT64_T, sender, 0, Buffer_Recv.data(), nVertexR, MPI_UINT64_T, + recver, 0, SU2_MPI::GetComm(), MPI_STATUS_IGNORE); /*--- Store received data---*/ for (unsigned long iVertex = 0; iVertex < nVertexR; iVertex++) @@ -135,45 +142,51 @@ void CPastixWrapper::Initialize(CGeometry* geometry, const CConfig* /*--- 3 - Copy, map the sparsity, and put it in Fortran numbering ---*/ - for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) { - colptr[iPoint] = static_cast(row_ptr[iPoint] + 1); + for (su2_index_t iPoint = 0; iPoint < nPointDomain; ++iPoint) { + colptr[static_cast(iPoint)] = CheckedCast( + CheckedAdd(row_ptr[iPoint], 1, "PaStiX row pointer entry"), "PaStiX row pointer entry"); - const unsigned long begin = row_ptr[iPoint], end = row_ptr[iPoint + 1]; + const su2_index_t begin = row_ptr[iPoint], end = row_ptr[iPoint + 1]; /*--- If last point of row is halo ---*/ const bool sort_required = (col_ind[end - 1] >= nPointDomain); if (sort_required) { - const unsigned long nnz_row = end - begin; + const su2_index_t nnz_row = end - begin; sort_rows.push_back(iPoint); - sort_order.emplace_back(nnz_row); + sort_order.emplace_back(CheckedCast(nnz_row, "PaStiX row sort order size")); /*--- Sort mapped indices ("first") and keep track of source ("second") for when we later need to swap blocks for these rows. ---*/ - vector > aux(nnz_row); + vector > aux(CheckedCast(nnz_row, "PaStiX row sort auxiliary size")); for (auto j = begin; j < end; ++j) { if (col_ind[j] < nPointDomain) { - aux[j - begin].first = static_cast(offset + col_ind[j] + 1); + aux[j - begin].first = CheckedCast(CheckedAdd(offset, col_ind[j] + 1, "PaStiX column index"), + "PaStiX column index"); } else { - aux[j - begin].first = static_cast(map[col_ind[j] - nPointDomain] + 1); + aux[j - begin].first = CheckedCast( + CheckedAdd(map[col_ind[j] - nPointDomain], 1, "PaStiX halo column index"), "PaStiX halo column index"); } aux[j - begin].second = j; } sort(aux.begin(), aux.end()); - for (auto j = 0ul; j < nnz_row; ++j) { + for (su2_index_t j = 0; j < nnz_row; ++j) { rowidx.push_back(aux[j].first); sort_order.back()[j] = aux[j].second; } } else { /*--- These are all internal, no need to go through map. ---*/ - for (auto j = begin; j < end; ++j) rowidx.push_back(static_cast(offset + col_ind[j] + 1)); + for (auto j = begin; j < end; ++j) + rowidx.push_back(CheckedCast(CheckedAdd(offset, col_ind[j] + 1, "PaStiX column index"), + "PaStiX column index")); } } - colptr[nPointDomain] = static_cast(nNonZero + 1); + colptr[static_cast(nPointDomain)] = + CheckedCast(CheckedAdd(nNonZero, 1, "PaStiX final colptr entry"), "PaStiX final colptr entry"); if (rowidx.size() != nNonZero) SU2_MPI::Error("Error during preparation of PaStiX data", CURRENT_FUNCTION); @@ -187,8 +200,8 @@ void CPastixWrapper::Initialize(CGeometry* geometry, const CConfig* spm.baseval = 1; spm.n = nCols; - spm.nnz = nNonZero; - spm.dof = nVar; + spm.nnz = CheckedCast(nNonZero, "PaStiX nonzero count"); + spm.dof = CheckedCast(nVar, "PaStiX block degree of freedom"); spm.colptr = colptr.data(); spm.rowptr = rowidx.data(); @@ -214,7 +227,7 @@ void CPastixWrapper::AssembleValues() { const auto nDomain = matrix.nPointDomain; const auto blkSz = matrix.blkSz; const auto *d = matrix.d, *l = matrix.l, *u = matrix.u; - for (auto iPoint = 0ul; iPoint < nDomain; ++iPoint) { + for (su2_index_t iPoint = 0; iPoint < nDomain; ++iPoint) { auto* dst = values.data() + csr_row_ptr[iPoint] * blkSz; for (auto k = matrix.row_ptr_l[iPoint]; k < matrix.row_ptr_l[iPoint + 1]; ++k, dst += blkSz) for (auto b = 0ul; b < blkSz; ++b) dst[b] = SU2_TYPE::GetValue(l[k * blkSz + b]); @@ -276,14 +289,14 @@ void CPastixWrapper::Factorize(CGeometry* geometry, const CConfig* c AssembleValues wrote them in LDU order; copy to tmp then write back sorted. ---*/ vector tmp; - for (auto i = 0ul; i < sort_rows.size(); ++i) { + for (size_t i = 0; i < sort_rows.size(); ++i) { const auto iRow = sort_rows[i]; /*--- colptr is 1-based Fortran numbering: row start = colptr[iRow] - 1. ---*/ - const auto begin = static_cast(colptr[iRow] - 1); + const auto begin = static_cast(colptr[iRow] - 1); const auto nnz_row = sort_order[i].size(); tmp.assign(values.begin() + begin * blkSz, values.begin() + (begin + nnz_row) * blkSz); - for (auto j = 0ul; j < nnz_row; ++j) { + for (size_t j = 0; j < nnz_row; ++j) { const auto src_pos = sort_order[i][j] - begin; for (auto k = 0ul; k < blkSz; ++k) values[(begin + j) * blkSz + k] = tmp[src_pos * blkSz + k]; } diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index c9668aa0e88..f652c1c9413 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -33,10 +33,50 @@ #include #include +#include #include #include namespace { +su2_index_t CheckedMul(su2_index_t lhs, su2_index_t rhs, const char* what) { + if (lhs != 0 && rhs > std::numeric_limits::max() / lhs) { + SU2_MPI::Error(std::string("Overflow while computing ") + what + ".", CURRENT_FUNCTION); + } + return lhs * rhs; +} + +su2_index_t CheckedAdd(su2_index_t lhs, su2_index_t rhs, const char* what) { + if (rhs > std::numeric_limits::max() - lhs) { + SU2_MPI::Error(std::string("Overflow while computing ") + what + ".", CURRENT_FUNCTION); + } + return lhs + rhs; +} + +template +size_t CheckedBytes(su2_index_t count, const char* what) { + const auto bytes = CheckedMul(count, sizeof(T), what); + if (bytes > std::numeric_limits::max()) { + SU2_MPI::Error(std::string("Overflow while computing ") + what + " byte size.", CURRENT_FUNCTION); + } + return static_cast(bytes); +} + +template +TargetType CheckedCast(su2_index_t value, const char* what) { + if (value > static_cast(std::numeric_limits::max())) { + SU2_MPI::Error(std::string("Overflow while converting ") + what + ".", CURRENT_FUNCTION); + } + return static_cast(value); +} + +std::vector ConvertSparseIndicesForDevice(const su2_index_t* src, su2_index_t count, const char* name) { + std::vector dst(CheckedCast(count, name)); + for (su2_index_t i = 0; i < count; ++i) { + dst[static_cast(i)] = CheckedCast(src[i], name); + } + return dst; +} + /*--- Helper function to regularize small pivots ---*/ template FORCEINLINE void RegularizePivot(ScalarType& pivot, unsigned long row, unsigned long col, const char* context) { @@ -122,7 +162,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); @@ -138,7 +178,7 @@ CSysMatrix::~CSysMatrix() { 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); @@ -226,10 +266,17 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi nPointDomain = npointdomain; /*--- Allocate host data. ---*/ - auto allocAndInit = [](ScalarType*& ptr, unsigned long num) { - ptr = MemoryAllocation::aligned_alloc(64, num * sizeof(ScalarType)); + auto allocAndInit = [](ScalarType*& ptr, su2_index_t num, const char* what) { + ptr = MemoryAllocation::aligned_alloc(64, CheckedBytes(num, what)); + }; + auto allocQ = [](QuantType*& ptr, su2_index_t num, const char* what) { + ptr = MemoryAllocation::aligned_alloc(64, CheckedBytes(num, what)); }; + const auto block_entries = CheckedMul(nVar, nEqn, "CSysMatrix block size"); + const auto diag_entries = CheckedMul(nPoint, block_entries, "CSysMatrix diagonal storage size"); + const auto domain_diag_entries = CheckedMul(nPointDomain, block_entries, "CSysMatrix domain diagonal storage size"); + /*--- L/D/U index structures and value arrays. ---*/ { const auto& pat = geometry->GetSparsePattern(type, 0); @@ -240,7 +287,7 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi mat.col_ind_u = pat.u.innerIdx(); mat.nnz_u = pat.u.getNumNonZeros(); } - allocAndInit(mat.d, nPoint * nVar * nEqn); + allocAndInit(mat.d, diag_entries, "CSysMatrix diagonal host allocation"); if (q_lus_needed) { /*--- Q_LU_SGS: no full-precision L/U; off-diagonal blocks live in quantized storage. @@ -248,38 +295,47 @@ 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)); - }; - 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, CheckedMul(mat.nnz_l, nVar, "CSysMatrix quantized L scale size"), + "CSysMatrix quantized L scale host allocation"); + allocQ(q_blocks_l, CheckedMul(mat.nnz_l, block_entries, "CSysMatrix quantized L block size"), + "CSysMatrix quantized L block host allocation"); + allocQ(q_scale_u, CheckedMul(mat.nnz_u, nVar, "CSysMatrix quantized U scale size"), + "CSysMatrix quantized U scale host allocation"); + allocQ(q_blocks_u, CheckedMul(mat.nnz_u, block_entries, "CSysMatrix quantized U block size"), + "CSysMatrix quantized U block host allocation"); + allocQ(q_scale_d, CheckedMul(nPoint, nVar, "CSysMatrix quantized diagonal scale size"), + "CSysMatrix quantized diagonal scale host allocation"); + allocQ(q_blocks_d, diag_entries, "CSysMatrix quantized diagonal block host allocation"); } else { - allocAndInit(mat.l, mat.nnz_l * nVar * nEqn); - allocAndInit(mat.u, mat.nnz_u * nVar * nEqn); + allocAndInit(mat.l, CheckedMul(mat.nnz_l, block_entries, "CSysMatrix L storage size"), + "CSysMatrix L host allocation"); + allocAndInit(mat.u, CheckedMul(mat.nnz_u, block_entries, "CSysMatrix U storage size"), + "CSysMatrix U host allocation"); } - auto GPUAllocAndInit = [](ScalarType*& ptr, unsigned long num) { - ptr = GPUMemoryAllocation::gpu_alloc(num * sizeof(ScalarType)); + auto GPUAllocAndInit = [](ScalarType*& ptr, su2_index_t num, const char* what) { + ptr = GPUMemoryAllocation::gpu_alloc(CheckedBytes(num, what)); }; - auto GPUAllocAndCopy = [](const su2uint*& ptr, const su2uint* src_ptr, unsigned long num) { - ptr = GPUMemoryAllocation::gpu_alloc_cpy(src_ptr, num * sizeof(su2uint)); + auto GPUAllocAndCopy = [](const su2uint*& ptr, const su2_index_t* src_ptr, su2_index_t num, const char* what) { + const auto host = ConvertSparseIndicesForDevice(src_ptr, num, what); + ptr = GPUMemoryAllocation::gpu_alloc_cpy(host.data(), CheckedBytes(num, what)); }; if (useCuda) { if (nVar != nEqn) { 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); + GPUAllocAndInit(gpu.d, diag_entries, "CSysMatrix diagonal device allocation"); + GPUAllocAndInit(gpu.l, CheckedMul(mat.nnz_l, block_entries, "CSysMatrix L device storage size"), + "CSysMatrix L device allocation"); + GPUAllocAndInit(gpu.u, CheckedMul(mat.nnz_u, block_entries, "CSysMatrix U device storage size"), + "CSysMatrix U device allocation"); + GPUAllocAndCopy(gpu.row_ptr_l, mat.row_ptr_l, CheckedAdd(nPointDomain, 1, "CSysMatrix L row_ptr"), + "CSysMatrix L row_ptr"); + GPUAllocAndCopy(gpu.col_ind_l, mat.col_ind_l, mat.nnz_l, "CSysMatrix L col_ind"); + GPUAllocAndCopy(gpu.row_ptr_u, mat.row_ptr_u, CheckedAdd(nPointDomain, 1, "CSysMatrix U row_ptr"), + "CSysMatrix U row_ptr"); + GPUAllocAndCopy(gpu.col_ind_u, mat.col_ind_u, mat.nnz_u, "CSysMatrix U col_ind"); } if (type == ConnectivityType::FiniteVolume) { @@ -317,18 +373,18 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi * instead. This does not change the elimination order/pattern (nothing here affects L/U * membership), only how the build is scheduled on the device. ---*/ if (useCuda) { - std::vector adjPtr(nPointDomain + 1, 0); - std::vector adjIdx; + std::vector adjPtr(nPointDomain + 1, 0); + std::vector adjIdx; adjIdx.reserve(ilu.nnz_l + ilu.nnz_u); for (auto i = 0ul; i < nPointDomain; ++i) { - adjPtr[i] = static_cast(adjIdx.size()); + adjPtr[i] = static_cast(adjIdx.size()); for (auto k = ilu.row_ptr_l[i]; k < ilu.row_ptr_l[i + 1]; ++k) adjIdx.push_back(ilu.col_ind_l[k]); for (auto k = ilu.row_ptr_u[i]; k < ilu.row_ptr_u[i + 1]; ++k) { const auto j = ilu.col_ind_u[k]; - if (j < nPointDomain) adjIdx.push_back(static_cast(j)); + if (j < nPointDomain) adjIdx.push_back(j); } } - adjPtr[nPointDomain] = static_cast(adjIdx.size()); + adjPtr[nPointDomain] = static_cast(adjIdx.size()); color_ilu = colorSparsePattern(CCompressedSparsePatternUL(adjPtr, adjIdx), 1, true, false); /*--- Report, across ranks, how many colors/levels the GPU ILU ends up scheduled over and how @@ -362,12 +418,14 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi /*--- Preconditioners. ---*/ if (ilu_needed) { - allocAndInit(ilu.l, ilu.nnz_l * nVar * nEqn); - allocAndInit(ilu.d, nPointDomain * nVar * nEqn); - allocAndInit(ilu.u, ilu.nnz_u * nVar * nEqn); + allocAndInit(ilu.l, CheckedMul(ilu.nnz_l, block_entries, "CSysMatrix ILU L storage size"), + "CSysMatrix ILU L host allocation"); + allocAndInit(ilu.d, domain_diag_entries, "CSysMatrix ILU diagonal host allocation"); + allocAndInit(ilu.u, CheckedMul(ilu.nnz_u, block_entries, "CSysMatrix ILU U storage size"), + "CSysMatrix ILU U host allocation"); } - if (diag_needed) allocAndInit(invM, nPointDomain * nVar * nEqn); + if (diag_needed) allocAndInit(invM, domain_diag_entries, "CSysMatrix inverse diagonal host allocation"); if (jacobi_on_device) { if (nVar != nEqn) { @@ -377,7 +435,8 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi SU2_MPI::Error("CUDA Jacobi preconditioner uses one thread per block entry, nVar is too large.", CURRENT_FUNCTION); } - d_invM = GPUMemoryAllocation::gpu_alloc(nPointDomain * nVar * nEqn * sizeof(ScalarType)); + d_invM = GPUMemoryAllocation::gpu_alloc( + CheckedBytes(domain_diag_entries, "CSysMatrix inverse diagonal device allocation")); } if (useCuda && ilu_needed) { @@ -391,13 +450,17 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi * are uploaded (once, here) because they do not change. ---*/ gpu_ilu.nnz_l = ilu.nnz_l; gpu_ilu.nnz_u = ilu.nnz_u; - GPUAllocAndInit(gpu_ilu.d, nPointDomain * nVar * nEqn); - GPUAllocAndInit(gpu_ilu.l, ilu.nnz_l * nVar * nEqn); - GPUAllocAndInit(gpu_ilu.u, ilu.nnz_u * nVar * nEqn); - GPUAllocAndCopy(gpu_ilu.row_ptr_l, ilu.row_ptr_l, nPointDomain + 1); - GPUAllocAndCopy(gpu_ilu.col_ind_l, ilu.col_ind_l, ilu.nnz_l); - GPUAllocAndCopy(gpu_ilu.row_ptr_u, ilu.row_ptr_u, nPointDomain + 1); - GPUAllocAndCopy(gpu_ilu.col_ind_u, ilu.col_ind_u, ilu.nnz_u); + GPUAllocAndInit(gpu_ilu.d, domain_diag_entries, "CSysMatrix ILU diagonal device allocation"); + GPUAllocAndInit(gpu_ilu.l, CheckedMul(ilu.nnz_l, block_entries, "CSysMatrix ILU L device storage size"), + "CSysMatrix ILU L device allocation"); + GPUAllocAndInit(gpu_ilu.u, CheckedMul(ilu.nnz_u, block_entries, "CSysMatrix ILU U device storage size"), + "CSysMatrix ILU U device allocation"); + GPUAllocAndCopy(gpu_ilu.row_ptr_l, ilu.row_ptr_l, CheckedAdd(nPointDomain, 1, "CSysMatrix ILU L row_ptr"), + "CSysMatrix ILU L row_ptr"); + GPUAllocAndCopy(gpu_ilu.col_ind_l, ilu.col_ind_l, ilu.nnz_l, "CSysMatrix ILU L col_ind"); + GPUAllocAndCopy(gpu_ilu.row_ptr_u, ilu.row_ptr_u, CheckedAdd(nPointDomain, 1, "CSysMatrix ILU U row_ptr"), + "CSysMatrix ILU U row_ptr"); + GPUAllocAndCopy(gpu_ilu.col_ind_u, ilu.col_ind_u, ilu.nnz_u, "CSysMatrix ILU U col_ind"); /*--- Flatten the coloring, the index type differs from the one of the pattern. It drives * the factorization on the device. ---*/ @@ -407,11 +470,12 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi ilu_color_ptr.push_back(0); for (auto color = 0ul; color < color_ilu.getOuterSize(); ++color) { for (auto k = 0ul; k < color_ilu.getNumNonZeros(color); ++k) { - color_idx.push_back(static_cast(color_ilu.getInnerIdx(color, k))); + color_idx.push_back(CheckedCast(color_ilu.getInnerIdx(color, k), "CSysMatrix ILU color index")); } - ilu_color_ptr.push_back(static_cast(color_idx.size())); + ilu_color_ptr.push_back(CheckedCast(color_idx.size(), "CSysMatrix ILU color pointer")); } - d_ilu_color_idx = GPUMemoryAllocation::gpu_alloc_cpy(color_idx.data(), color_idx.size() * sizeof(su2uint)); + d_ilu_color_idx = GPUMemoryAllocation::gpu_alloc_cpy( + color_idx.data(), CheckedBytes(color_idx.size(), "CSysMatrix ILU color device copy")); /*--- Flatten levels_ilu the same way. It drives both triangular solves on the device. ---*/ std::vector level_idx; @@ -420,11 +484,12 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi ilu_level_ptr.push_back(0); for (auto level = 0ul; level < levels_ilu.getOuterSize(); ++level) { for (auto k = 0ul; k < levels_ilu.getNumNonZeros(level); ++k) { - level_idx.push_back(static_cast(levels_ilu.getInnerIdx(level, k))); + level_idx.push_back(CheckedCast(levels_ilu.getInnerIdx(level, k), "CSysMatrix ILU level index")); } - ilu_level_ptr.push_back(static_cast(level_idx.size())); + ilu_level_ptr.push_back(CheckedCast(level_idx.size(), "CSysMatrix ILU level pointer")); } - d_ilu_level_idx = GPUMemoryAllocation::gpu_alloc_cpy(level_idx.data(), level_idx.size() * sizeof(su2uint)); + d_ilu_level_idx = GPUMemoryAllocation::gpu_alloc_cpy( + level_idx.data(), CheckedBytes(level_idx.size(), "CSysMatrix ILU level device copy")); } /*--- Thread parallel initialization. ---*/ @@ -433,7 +498,7 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi /*--- Set suitable chunk sizes for light static for loops, and heavy dynamic ones, such that threads are approximately evenly loaded. ---*/ - omp_light_size = computeStaticChunkSize(nPoint * nVar * nEqn, num_threads, OMP_MAX_SIZE_L); + omp_light_size = computeStaticChunkSize(diag_entries, num_threads, OMP_MAX_SIZE_L); omp_heavy_size = computeStaticChunkSize(nPointDomain, num_threads, OMP_MAX_SIZE_H); omp_num_parts = config->GetLinear_Solver_Prec_Threads(); @@ -746,22 +811,23 @@ void CSysMatrix::SetValZero() { const auto nThreads = static_cast(omp_get_num_threads()); const auto iThread = static_cast(omp_get_thread_num()); - auto zeroChunk = [&](auto* arr, unsigned long n) { + auto zeroChunk = [&](auto* arr, su2_index_t n) { if (n == 0) return; const auto chunk = roundUpDiv(n, nThreads); const auto begin = min(chunk * iThread, n); const auto mySize = min(chunk, n - begin) * sizeof(std::remove_pointer_t); if (mySize) memset(&arr[begin], 0, mySize); }; - zeroChunk(mat.d, nPoint * nVar * nEqn); + const auto block_entries = CheckedMul(nVar, nEqn, "CSysMatrix zero-fill block size"); + zeroChunk(mat.d, CheckedMul(nPoint, block_entries, "CSysMatrix zero-fill diagonal size")); if (!quantized_mode) { - zeroChunk(mat.l, mat.nnz_l * nVar * nEqn); - zeroChunk(mat.u, mat.nnz_u * nVar * nEqn); + zeroChunk(mat.l, CheckedMul(mat.nnz_l, block_entries, "CSysMatrix zero-fill L size")); + zeroChunk(mat.u, CheckedMul(mat.nnz_u, block_entries, "CSysMatrix zero-fill U size")); } 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, CheckedMul(mat.nnz_l, nVar, "CSysMatrix zero-fill quantized L scale size")); + zeroChunk(q_scale_u, CheckedMul(mat.nnz_u, nVar, "CSysMatrix zero-fill quantized U scale size")); + zeroChunk(q_blocks_l, CheckedMul(mat.nnz_l, block_entries, "CSysMatrix zero-fill quantized L block size")); + zeroChunk(q_blocks_u, CheckedMul(mat.nnz_u, block_entries, "CSysMatrix zero-fill quantized U block size")); } SU2_OMP_BARRIER } @@ -769,8 +835,10 @@ void CSysMatrix::SetValZero() { template void CSysMatrix::SetValDiagonalZero() { SU2_ZONE_SCOPED + const auto nEntries = CheckedMul(CheckedMul(nPointDomain, nVar, "CSysMatrix diagonal zero-fill size"), nEqn, + "CSysMatrix diagonal zero-fill size"); SU2_OMP_FOR_STAT(omp_heavy_size) - for (auto iVar = 0ul; iVar < nPointDomain * nVar * nEqn; ++iVar) mat.d[iVar] = 0; + for (su2_index_t iVar = 0; iVar < nEntries; ++iVar) mat.d[iVar] = 0; END_SU2_OMP_FOR } @@ -1015,9 +1083,9 @@ void CSysMatrix::BuildILUPreconditioner() { if (ilu_fill_in == 0) { /*--- ILU0: Same sparse pattern, copy L and U blocks directly. ---*/ - auto copy = [&](const su2uint* row_ptr, const ScalarType* mat, ScalarType* ilu) { - const unsigned long begin = row_ptr[iPoint] * blockSize; - const unsigned long end = row_ptr[iPoint + 1] * blockSize; + auto copy = [&](const su2_index_t* row_ptr, const ScalarType* mat, ScalarType* ilu) { + const su2_index_t begin = row_ptr[iPoint] * blockSize; + const su2_index_t end = row_ptr[iPoint + 1] * blockSize; SU2_OMP_SIMD for (auto k = begin; k < end; ++k) ilu[k] = mat[k]; }; @@ -1026,8 +1094,8 @@ void CSysMatrix::BuildILUPreconditioner() { return; } /*--- ILUn: Merge-scan L and U via shared lambda. ---*/ - auto scatterPart = [&](const su2uint* mat_row_ptr, const su2uint* mat_col_ind, const ScalarType* mat_vals, - const su2uint* ilu_row_ptr, const su2uint* ilu_col_ind, ScalarType* ilu_vals) { + auto scatterPart = [&](const su2_index_t* mat_row_ptr, const su2_index_t* mat_col_ind, const ScalarType* mat_vals, + const su2_index_t* ilu_row_ptr, const su2_index_t* ilu_col_ind, ScalarType* ilu_vals) { auto km = mat_row_ptr[iPoint], km_end = mat_row_ptr[iPoint + 1]; for (auto k = ilu_row_ptr[iPoint]; k < ilu_row_ptr[iPoint + 1]; ++k) { const auto jPoint = ilu_col_ind[k]; @@ -1593,7 +1661,7 @@ void CSysMatrix::TransposeInPlace() { if (edge_ptr_l) { /*--- FV path: each edge maps to one U and one L block. ---*/ SU2_OMP_FOR_DYN(omp_heavy_size * 2) - for (auto iEdge = 0ul; iEdge < mat.nnz_l; ++iEdge) { + for (su2_index_t iEdge = 0; iEdge < mat.nnz_l; ++iEdge) { auto* bij_u = &mat.u[iEdge * nVar * nVar]; auto* bji_l = &mat.l[edge_ptr_l[iEdge] * nVar * nVar]; swapAndTransp(nVar, bij_u, bji_l); @@ -1650,14 +1718,18 @@ void CSysMatrix::MatrixMatrixAddition(ScalarType alpha, const CSysMa (nEqn == B.nEqn) && (nPoint == B.nPoint) && (mat.nnz_l == B.mat.nnz_l) && (mat.nnz_u == B.mat.nnz_u); if (!ok) SU2_MPI::Error("Matrices do not have compatible sparsity.", CURRENT_FUNCTION); + const auto block_entries = CheckedMul(nVar, nEqn, "CSysMatrix matrix sum block size"); + const auto diag_entries = CheckedMul(nPoint, block_entries, "CSysMatrix matrix sum diagonal size"); + const auto l_entries = CheckedMul(mat.nnz_l, block_entries, "CSysMatrix matrix sum L size"); + const auto u_entries = CheckedMul(mat.nnz_u, block_entries, "CSysMatrix matrix sum U size"); SU2_OMP_FOR_STAT(omp_light_size) - for (auto i = 0ul; i < nPoint * nVar * nEqn; ++i) mat.d[i] += alpha * B.mat.d[i]; + for (su2_index_t i = 0; i < diag_entries; ++i) mat.d[i] += alpha * B.mat.d[i]; END_SU2_OMP_FOR SU2_OMP_FOR_STAT(omp_light_size) - for (auto i = 0ul; i < mat.nnz_l * nVar * nEqn; ++i) mat.l[i] += alpha * B.mat.l[i]; + for (su2_index_t i = 0; i < l_entries; ++i) mat.l[i] += alpha * B.mat.l[i]; END_SU2_OMP_FOR SU2_OMP_FOR_STAT(omp_light_size) - for (auto i = 0ul; i < mat.nnz_u * nVar * nEqn; ++i) mat.u[i] += alpha * B.mat.u[i]; + for (su2_index_t i = 0; i < u_entries; ++i) mat.u[i] += alpha * B.mat.u[i]; END_SU2_OMP_FOR } diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index 1cd5de51d5c..fbab85a6ddb 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -26,6 +26,8 @@ */ #include +#include +#include #include "../../include/linear_algebra/CMatrixInverse.hpp" #include "../../include/linear_algebra/CSysMatrix.inl" @@ -353,6 +355,29 @@ __global__ void BlockLDU_SpMV_kernel(unsigned long nRows, unsigned long nVar, y[iRow * nVar + iVar] = sum; } +su2_index_t CheckedMul(su2_index_t lhs, su2_index_t rhs, const char* what) { + if (lhs != 0 && rhs > std::numeric_limits::max() / lhs) { + SU2_MPI::Error(std::string("Overflow while computing ") + what + ".", CURRENT_FUNCTION); + } + return lhs * rhs; +} + +template +size_t CheckedBytes(su2_index_t count, const char* what) { + const auto bytes = CheckedMul(count, sizeof(T), what); + if (bytes > std::numeric_limits::max()) { + SU2_MPI::Error(std::string("Overflow while computing ") + what + " byte size.", CURRENT_FUNCTION); + } + return static_cast(bytes); +} + +unsigned CheckedCudaDim(su2_index_t value, const char* what) { + if (value > static_cast(std::numeric_limits::max())) { + SU2_MPI::Error(std::string(what) + " exceeds the CUDA grid/block dimension range.", CURRENT_FUNCTION); + } + return static_cast(value); +} + } // namespace template @@ -369,7 +394,7 @@ void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector((nPointDomain + threadsPerBlock - 1) / threadsPerBlock); + const auto blocks = CheckedCudaDim(roundUpDiv(nPointDomain, threadsPerBlock), "CUDA Jacobi preconditioner grid"); ApplyJacobiPreconditionerKernel<<>>(d_invM, vec.GetDevicePointer(), prod.GetDevicePointer(), nPointDomain, nVar); /*--- Sync so the zone above actually times the kernel, not just the (async) launch call. ---*/ @@ -388,10 +413,12 @@ 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); + const auto blockSize = CheckedCudaDim(CheckedMul(nVar, nVar, "CUDA Jacobi preconditioner block size"), + "CUDA Jacobi preconditioner block"); InvertDiagonalBlocksKernel - <<(nPointDomain), blockSize, blockSize * sizeof(ScalarType)>>>(nPointDomain, nVar, gpu.d, - d_invM); + <<(blockSize, "CUDA Jacobi preconditioner shared memory")>>>(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()); @@ -413,8 +440,9 @@ void CSysMatrix::BuildILUPreconditionerGPU() { const DeviceLDU M{gpu_ilu.d, gpu_ilu.l, gpu_ilu.u, gpu_ilu.row_ptr_l, gpu_ilu.col_ind_l, gpu_ilu.row_ptr_u, gpu_ilu.col_ind_u}; - const auto blockSize = static_cast(nVar * nVar); - const auto shared = 2 * blockSize * sizeof(ScalarType); + const auto blockSize = + CheckedCudaDim(CheckedMul(nVar, nVar, "CUDA ILU factorization block size"), "CUDA ILU factorization block"); + const auto shared = CheckedBytes(2 * blockSize, "CUDA ILU factorization shared memory"); /*--- 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 @@ -477,9 +505,11 @@ void CSysMatrix::ComputeILUPreconditionerGPU(const CSysVector(nVar * nVar); - const auto sharedForward = threads * sizeof(ScalarType); - const auto sharedBackward = (threads + nVar) * sizeof(ScalarType); + const auto threads = + CheckedCudaDim(CheckedMul(nVar, nVar, "CUDA ILU solve block size"), "CUDA ILU solve block"); + const auto sharedForward = CheckedBytes(threads, "CUDA ILU forward shared memory"); + const auto sharedBackward = + CheckedBytes(threads + nVar, "CUDA ILU backward shared memory"); if (ilu_stream == nullptr) gpuErrChk(cudaStreamCreate(&ilu_stream)); @@ -534,9 +564,19 @@ template 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)); + const auto block_entries = CheckedMul(nVar, nEqn, "CSysMatrix GPU transfer block size"); + gpuErrChk(cudaMemcpy(gpu.d, mat.d, + CheckedBytes(CheckedMul(nPoint, block_entries, "CSysMatrix GPU diagonal copy"), + "CSysMatrix GPU diagonal copy"), + cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpy(gpu.l, mat.l, + CheckedBytes(CheckedMul(mat.nnz_l, block_entries, "CSysMatrix GPU L copy"), + "CSysMatrix GPU L copy"), + cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpy(gpu.u, mat.u, + CheckedBytes(CheckedMul(mat.nnz_u, block_entries, "CSysMatrix GPU U copy"), + "CSysMatrix GPU U copy"), + cudaMemcpyHostToDevice)); } template @@ -547,8 +587,8 @@ 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); + dim3 blockDim(CheckedCudaDim(nVar, "CUDA matrix-vector block"), 1, 1); + dim3 gridDim(CheckedCudaDim(nPointDomain, "CUDA matrix-vector grid"), 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); diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp index 54e47157d55..22d2afae9e8 100644 --- a/Common/src/linear_algebra/CSysVector.cpp +++ b/Common/src/linear_algebra/CSysVector.cpp @@ -28,8 +28,31 @@ #include "../../include/linear_algebra/CSysVector.hpp" #include "../../include/toolboxes/allocation_toolbox.hpp" +#include +#include + +namespace { + +su2_index_t CheckedMul(su2_index_t lhs, su2_index_t rhs, const char* what) { + if (lhs != 0 && rhs > std::numeric_limits::max() / lhs) { + SU2_MPI::Error(std::string("Overflow while computing ") + what + ".", CURRENT_FUNCTION); + } + return lhs * rhs; +} + +template +size_t CheckedBytes(su2_index_t count, const char* what) { + const auto bytes = CheckedMul(count, sizeof(T), what); + if (bytes > std::numeric_limits::max()) { + SU2_MPI::Error(std::string("Overflow while computing ") + what + " byte size.", CURRENT_FUNCTION); + } + return static_cast(bytes); +} + +} // namespace + template -void CSysVector::Initialize(unsigned long numBlk, unsigned long numBlkDomain, unsigned long numVar, +void CSysVector::Initialize(su2_index_t numBlk, su2_index_t numBlkDomain, su2_index_t numVar, const ScalarType* val, bool valIsArray, bool errorIfParallel) { if (errorIfParallel && omp_in_parallel()) { assert(false); @@ -39,23 +62,29 @@ void CSysVector::Initialize(unsigned long numBlk, unsigned long numB if (omp_get_thread_num()) SU2_MPI::Error("Only the master thread is allowed to initialize the vector.", CURRENT_FUNCTION); - if (nElm != numBlk * numVar) { + const auto new_nElm = CheckedMul(numBlk, numVar, "CSysVector storage size"); + const auto new_nElmDomain = CheckedMul(numBlkDomain, numVar, "CSysVector domain storage size"); + + if (nElm != new_nElm) { MemoryAllocation::aligned_free(vec_val); vec_val = nullptr; } - nElm = numBlk * numVar; - nElmDomain = numBlkDomain * numVar; + nElm = new_nElm; + nElmDomain = new_nElmDomain; nVar = numVar; omp_chunk_size = computeStaticChunkSize(nElm, omp_get_max_threads(), OMP_MAX_SIZE); - if (vec_val == nullptr) vec_val = MemoryAllocation::aligned_alloc(64, nElm * sizeof(ScalarType)); + if (vec_val == nullptr) + vec_val = MemoryAllocation::aligned_alloc( + 64, CheckedBytes(nElm, "CSysVector host allocation")); /*--- Device storage mirrors the host allocation; free first so that re-initializing a * vector does not leak it. ---*/ GPUMemoryAllocation::gpu_free(d_vec_val); - d_vec_val = GPUMemoryAllocation::gpu_alloc(nElm * sizeof(ScalarType)); + d_vec_val = + GPUMemoryAllocation::gpu_alloc(CheckedBytes(nElm, "CSysVector device allocation")); #ifdef HAVE_OMP dot_scratch.reset(new ScalarType[omp_get_max_threads()]); @@ -63,9 +92,9 @@ void CSysVector::Initialize(unsigned long numBlk, unsigned long numB if (val != nullptr) { if (!valIsArray) { - for (auto i = 0ul; i < nElm; i++) vec_val[i] = *val; + for (su2_index_t i = 0; i < nElm; i++) vec_val[i] = *val; } else { - for (auto i = 0ul; i < nElm; i++) vec_val[i] = val[i]; + for (su2_index_t i = 0; i < nElm; i++) vec_val[i] = val[i]; } } } @@ -155,7 +184,7 @@ const su2matrix& CSysVector::multiDot(const std::vector< template CSysVector::~CSysVector() { if constexpr (!std::is_trivial_v) { - for (auto i = 0ul; i < nElm; i++) vec_val[i].~ScalarType(); + for (su2_index_t i = 0; i < nElm; i++) vec_val[i].~ScalarType(); } MemoryAllocation::aligned_free(vec_val); diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 01b146baf76..0d3d3e07313 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -29,6 +29,8 @@ #include "../../include/linear_algebra/GPUComms.cuh" #include #include +#include +#include #include #include @@ -47,6 +49,29 @@ cublasHandle_t GetBlasHandle() { return solver_blas_handle; } +su2_index_t CheckedMul(su2_index_t lhs, su2_index_t rhs, const char* what) { + if (lhs != 0 && rhs > std::numeric_limits::max() / lhs) { + SU2_MPI::Error(std::string("Overflow while computing ") + what + ".", CURRENT_FUNCTION); + } + return lhs * rhs; +} + +template +size_t CheckedBytes(su2_index_t count, const char* what) { + const auto bytes = CheckedMul(count, sizeof(T), what); + if (bytes > std::numeric_limits::max()) { + SU2_MPI::Error(std::string("Overflow while computing ") + what + " byte size.", CURRENT_FUNCTION); + } + return static_cast(bytes); +} + +int CheckedCublasSize(su2_index_t size, const char* what) { + if (size > static_cast(std::numeric_limits::max())) { + SU2_MPI::Error(std::string(what) + " exceeds the cuBLAS int vector-length range.", CURRENT_FUNCTION); + } + return static_cast(size); +} + } // namespace namespace VecExpr { @@ -68,14 +93,16 @@ template void CSysVector::HtDTransfer(bool trigger) const { SU2_ZONE_SCOPED if (trigger) - gpuErrChk(cudaMemcpy((void*)(d_vec_val), (void*)&vec_val[0], (sizeof(ScalarType) * nElm), cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpy((void*)(d_vec_val), (void*)&vec_val[0], + CheckedBytes(nElm, "CSysVector host-to-device copy"), cudaMemcpyHostToDevice)); } template void CSysVector::DtHTransfer(bool trigger) const { SU2_ZONE_SCOPED if (trigger) - gpuErrChk(cudaMemcpy((void*)(&vec_val[0]), (void*)d_vec_val, (sizeof(ScalarType) * nElm), cudaMemcpyDeviceToHost)); + gpuErrChk(cudaMemcpy((void*)(&vec_val[0]), (void*)d_vec_val, + CheckedBytes(nElm, "CSysVector device-to-host copy"), cudaMemcpyDeviceToHost)); } template @@ -89,10 +116,12 @@ ScalarType CSysVector::GPUDot(const CSysVector& other) const { ScalarType local_dot = ScalarType(0); if constexpr (std::is_same_v) { - status = cublasSdot(handle, static_cast(nElmDomain), GetDevicePointer(), 1, other.GetDevicePointer(), 1, + status = cublasSdot(handle, CheckedCublasSize(nElmDomain, "CSysVector::GPUDot size"), 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, + status = cublasDdot(handle, CheckedCublasSize(nElmDomain, "CSysVector::GPUDot size"), GetDevicePointer(), 1, + other.GetDevicePointer(), 1, &local_dot); } else { SU2_MPI::Error("Unsupported ScalarType in CSysVector::GPUDot.", CURRENT_FUNCTION); @@ -190,7 +219,7 @@ using DeviceBcgsDir = VecExpr::add_, Vec, S>; #define INSTANTIATE_DEVICE_ASSIGN(SCALAR, OP, EXPR) \ template void VecExpr::AssignDeviceExpression>( \ - SCALAR*, unsigned long, const VecExpr::CVecExpr, SCALAR>&) + SCALAR*, su2_index_t, const VecExpr::CVecExpr, SCALAR>&) #define INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, EXPR) \ INSTANTIATE_DEVICE_ASSIGN(SCALAR, Assign, EXPR); \ diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index c2b460f8346..887888db2ff 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -86,11 +86,11 @@ class CFEASolver : public CFEASolverBase { CProperty** element_properties = nullptr; /*!< \brief Vector which stores the properties of each element */ #ifdef HAVE_OMP - vector > ElemColoring; /*!< \brief Element colors. */ + vector > ElemColoring; /*!< \brief Element colors. */ bool LockStrategy = false; /*!< \brief Whether to use an OpenMP lock to guard updates of the Jacobian. */ vector UpdateLocks; /*!< \brief Locks that may be used to protect accesses to CSysMatrix/Vector in element loops. */ #else - array,1> ElemColoring; /*--- Behaves like a normal integer type. ---*/ + array,1> ElemColoring; /*--- Behaves like a normal integer type. ---*/ static constexpr bool LockStrategy = false; /*--- Lock strategy is never needed for MPI-only. ---*/ DummyVectorOfLocks UpdateLocks; #endif diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index b2dd6441318..32d8c89ea91 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -197,10 +197,10 @@ class CFVMFlowSolverBase : public CSolver { /*--- Shallow copy of grid coloring for OpenMP parallelization. ---*/ #ifdef HAVE_OMP - vector > EdgeColoring; /*!< \brief Edge colors. */ + vector > EdgeColoring; /*!< \brief Edge colors. */ bool ReducerStrategy = false; /*!< \brief If the reducer strategy is in use. */ #else - array, 1> EdgeColoring; + array, 1> EdgeColoring; /*--- Never use the reducer strategy if compiling for MPI-only. ---*/ static constexpr bool ReducerStrategy = false; #endif diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 30dc89f0a43..b6535d0d1b5 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -372,7 +372,7 @@ void CFVMFlowSolverBase::HybridParallelInitialization(const CConfig& confi omp_chunk_size = computeStaticChunkSize(nPoint, omp_get_max_threads(), OMP_MAX_SIZE); #else - EdgeColoring[0] = DummyGridColor<>(geometry.GetnEdge()); + EdgeColoring[0] = DummyGridColor(geometry.GetnEdge()); #endif } diff --git a/SU2_CFD/include/solvers/CScalarSolver.hpp b/SU2_CFD/include/solvers/CScalarSolver.hpp index ce391414b0b..56a7cd86b39 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.hpp +++ b/SU2_CFD/include/solvers/CScalarSolver.hpp @@ -67,10 +67,10 @@ class CScalarSolver : public CSolver { /*--- Shallow copy of grid coloring for OpenMP parallelization. ---*/ #ifdef HAVE_OMP - vector > EdgeColoring; /*!< \brief Edge colors. */ + vector > EdgeColoring; /*!< \brief Edge colors. */ bool ReducerStrategy = false; /*!< \brief If the reducer strategy is in use. */ #else - array, 1> EdgeColoring; + array, 1> EdgeColoring; /*--- Never use the reducer strategy if compiling for MPI-only. ---*/ static constexpr bool ReducerStrategy = false; #endif diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index d6e05191b22..c2dcb34f728 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -74,7 +74,7 @@ CScalarSolver::CScalarSolver(CGeometry* geometry, CConfig* config, nPoint = geometry->GetnPoint(); omp_chunk_size = computeStaticChunkSize(nPoint, omp_get_max_threads(), OMP_MAX_SIZE); #else - EdgeColoring[0] = DummyGridColor<>(geometry->GetnEdge()); + EdgeColoring[0] = DummyGridColor(geometry->GetnEdge()); #endif /*--- Initialize lower and upper limits for solution clipping. Solvers might overwrite these values. ---*/ diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 27ad9c8fdff..a4ac44c7f83 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -282,7 +282,7 @@ void CFEASolver::HybridParallelInitialization(CGeometry* geometry) { omp_chunk_size = computeStaticChunkSize(nPointDomain, omp_get_max_threads(), OMP_MAX_SIZE); #else - ElemColoring[0] = DummyGridColor<>(nElement); + ElemColoring[0] = DummyGridColor(nElement); #endif }