Harden host-side sparse index width and allocation safety

Common/include/basic_types/datatype_structure.hpp

@@ -29,11 +29,14 @@
 
 #include <iostream>
 #include <complex>
+#include <cstdint>
 #include <cstdio>
 
 #include "../code_config.hpp"
 #include "ad_structure.hpp"
 
+using su2_index_t = std::uint64_t;
+
 /*!
  * \namespace SU2_TYPE
  * \brief Namespace for defining the datatype wrapper routines, this acts as a base

Common/include/geometry/CGeometry.hpp

@@ -1882,7 +1882,7 @@ class CGeometry {
    * \param[in] type - Finite volume or finite element.
    * \return Reference to the map.
    */
-  const su2vector<unsigned long>& GetTransposeSparsePatternMap(ConnectivityType type);
+  const su2vector<su2_index_t>& GetTransposeSparsePatternMap(ConnectivityType type);
 
   /*!
    * \brief Get the edge coloring.

Common/include/linear_algebra/CPastixWrapper.hpp

@@ -36,8 +36,11 @@
 
 #include <pastix.h>
 #include <spm.h>
+#include <limits>
 #include <vector>
 
+#include "../basic_types/datatype_structure.hpp"
+
 using namespace std;
 
 class CConfig;
@@ -66,14 +69,19 @@ class CPastixWrapper {
   double dparm[DPARM_SIZE];       /*!< \brief Floating point parameters for PaStiX. */
 
   struct {
-    unsigned long nVar = 0;
-    unsigned long nPoint = 0;
-    unsigned long nPointDomain = 0;
-    const unsigned long* rowptr = nullptr;
-    const unsigned long* colidx = nullptr;
+    su2_index_t nVar = 0;
+    su2_index_t nPoint = 0;
+    su2_index_t nPointDomain = 0;
+    const su2_index_t* rowptr = nullptr;
+    const su2_index_t* colidx = nullptr;
     const ScalarType* values = nullptr;
 
-    unsigned long size_rhs() const { return nPointDomain * nVar; }
+    su2_index_t size_rhs() const {
+      if (nPointDomain != 0 && nVar > std::numeric_limits<su2_index_t>::max() / nPointDomain) {
+        SU2_MPI::Error("Overflow while computing PaStiX rhs size.", CURRENT_FUNCTION);
+      }
+      return nPointDomain * nVar;
+    }
   } matrix; /*!< \brief Pointers and sizes of the input matrix. */
 
   bool issetup{};        /*!< \brief Signals that the matrix data has been provided. */
@@ -85,8 +93,8 @@ class CPastixWrapper {
   const int mpi_size = SU2_MPI::GetSize();
   const int mpi_rank = SU2_MPI::GetRank();
 
-  vector<unsigned long> sort_rows;          /*!< \brief List of rows with halo points. */
-  vector<vector<unsigned long>> sort_order; /*!< \brief How each of those rows needs to be sorted. */
+  vector<su2_index_t> sort_rows;          /*!< \brief List of rows with halo points. */
+  vector<vector<su2_index_t>> sort_order; /*!< \brief How each of those rows needs to be sorted. */
 
   /*!
    * \brief Run the "clean" task, releases all memory, leaves object in unusable state.
@@ -133,8 +141,8 @@ class CPastixWrapper {
    * \param[in] colidx - Non zeros column indices.
    * \param[in] values - Matrix coefficients.
    */
-  void SetMatrix(unsigned long nVar, unsigned long nPoint, unsigned long nPointDomain, const unsigned long* rowptr,
-                 const unsigned long* colidx, const ScalarType* values) {
+  void SetMatrix(su2_index_t nVar, su2_index_t nPoint, su2_index_t nPointDomain, const su2_index_t* rowptr,
+                 const su2_index_t* colidx, const ScalarType* values) {
     if (issetup) return;
     matrix.nVar = nVar;
     matrix.nPoint = nPoint;
@@ -176,11 +184,17 @@ 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();
+    if (rhs_size > static_cast<su2_index_t>(std::numeric_limits<pastix_int_t>::max())) {
+      SU2_MPI::Error("Overflow while converting PaStiX rhs size.", CURRENT_FUNCTION);
+    }
+    const auto pastix_rhs_size = static_cast<pastix_int_t>(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

Common/include/linear_algebra/CSysMatrix.hpp

@@ -138,25 +138,25 @@ class CSysMatrix {
   unsigned long nVar;         /*!< \brief Number of variables (and rows of the blocks). */
   unsigned long nEqn;         /*!< \brief Number of equations (and columns of the blocks). */
 
-  ScalarType* matrix;           /*!< \brief Entries of the sparse matrix. */
-  unsigned long nnz;            /*!< \brief Number of possible nonzero entries in the matrix. */
-  const unsigned long* row_ptr; /*!< \brief Pointers to the first element in each row. */
-  const unsigned long* dia_ptr; /*!< \brief Pointers to the diagonal element in each row. */
-  const unsigned long* col_ind; /*!< \brief Column index for each of the elements in val(). */
-  const unsigned long* col_ptr; /*!< \brief The transpose of col_ind, pointer to blocks with the same column index. */
+  ScalarType* matrix;         /*!< \brief Entries of the sparse matrix. */
+  su2_index_t nnz;            /*!< \brief Number of possible nonzero entries in the matrix. */
+  const su2_index_t* row_ptr; /*!< \brief Pointers to the first element in each row. */
+  const su2_index_t* dia_ptr; /*!< \brief Pointers to the diagonal element in each row. */
+  const su2_index_t* col_ind; /*!< \brief Column index for each of the elements in val(). */
+  const su2_index_t* col_ptr; /*!< \brief The transpose of col_ind, pointer to blocks with the same column index. */
 
   ScalarType* d_matrix;           /*!< \brief Device Pointer to store the matrix values on the GPU. */
   const unsigned long* d_row_ptr; /*!< \brief Device Pointers to the first element in each row. */
   const unsigned long* d_col_ind; /*!< \brief Device Column index for each of the elements in val(). */
   bool useCuda = false;           /*!< \brief Boolean that indicates whether user has enabled CUDA or not.
                                      Mainly used to conditionally free GPU memory in the class destructor. */
 
-  ScalarType* ILU_matrix;           /*!< \brief Entries of the ILU sparse matrix. */
-  unsigned long nnz_ilu;            /*!< \brief Number of possible nonzero entries in the matrix (ILU). */
-  const unsigned long* row_ptr_ilu; /*!< \brief Pointers to the first element in each row (ILU). */
-  const unsigned long* dia_ptr_ilu; /*!< \brief Pointers to the diagonal element in each row (ILU). */
-  const unsigned long* col_ind_ilu; /*!< \brief Column index for each of the elements in val() (ILU). */
-  unsigned short ilu_fill_in;       /*!< \brief Fill in level for the ILU preconditioner. */
+  ScalarType* ILU_matrix;         /*!< \brief Entries of the ILU sparse matrix. */
+  su2_index_t nnz_ilu;            /*!< \brief Number of possible nonzero entries in the matrix (ILU). */
+  const su2_index_t* row_ptr_ilu; /*!< \brief Pointers to the first element in each row (ILU). */
+  const su2_index_t* dia_ptr_ilu; /*!< \brief Pointers to the diagonal element in each row (ILU). */
+  const su2_index_t* col_ind_ilu; /*!< \brief Column index for each of the elements in val() (ILU). */
+  unsigned short ilu_fill_in;     /*!< \brief Fill in level for the ILU preconditioner. */
 
   /*!< \brief Level structure for alternative shared memory parallelization of ILU. */
   CCompressedSparsePatternUL levels_ilu;
@@ -192,14 +192,14 @@ class CSysMatrix {
    * \brief Auxilary object to wrap the edge map pointer used in fast block updates, i.e. without linear searches.
    */
   struct {
-    const unsigned long* ptr = nullptr;
+    const su2_index_t* ptr = nullptr;
     unsigned long nEdge = 0;
 
     operator bool() { return nEdge != 0; }
 
-    inline unsigned long operator()(unsigned long edge, unsigned long node) const { return ptr[2 * edge + node]; }
-    inline unsigned long ij(unsigned long edge) const { return ptr[2 * edge]; }
-    inline unsigned long ji(unsigned long edge) const { return ptr[2 * edge + 1]; }
+    inline su2_index_t operator()(unsigned long edge, unsigned long node) const { return ptr[2 * edge + node]; }
+    inline su2_index_t ij(unsigned long edge) const { return ptr[2 * edge]; }
+    inline su2_index_t ji(unsigned long edge) const { return ptr[2 * edge + 1]; }
 
   } edge_ptr;
 
@@ -817,7 +817,8 @@ class CSysMatrix {
    */
   template <class OtherType>
   inline void SetVal2Diag(unsigned long block_i, OtherType val_matrix) {
-    unsigned long iVar, index = dia_ptr[block_i] * nVar * nVar;
+    unsigned long iVar;
+    const su2_index_t index = dia_ptr[block_i] * nVar * nVar;
 
     /*--- Clear entire block before setting its diagonal. ---*/
     SU2_OMP_SIMD

Common/include/linear_algebra/CSysVector.hpp

@@ -34,6 +34,7 @@
 #include "../parallelization/mpi_structure.hpp"
 #include "../parallelization/omp_structure.hpp"
 #include "../parallelization/vectorization.hpp"
+#include "../basic_types/datatype_structure.hpp"
 #include "vector_expressions.hpp"
 #include "../../include/CConfig.hpp"
 
@@ -72,9 +73,9 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, 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. */
 
@@ -95,7 +96,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, 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);
 
   /*!
@@ -128,7 +129,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, 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).
@@ -137,7 +138,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, 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);
   }
 
@@ -146,7 +147,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, 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).
@@ -155,7 +156,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, 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);
   }
 
@@ -186,7 +187,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, 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);
   }
 
@@ -198,7 +199,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, 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);
   }
 
@@ -218,7 +219,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, 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
   }
 
@@ -248,35 +249,35 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, 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.
@@ -302,7 +303,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    */
   CSysVector& operator=(const CSysVector& other) {
     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;
   }
@@ -311,19 +312,19 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \brief Compound assignement operations with scalars and expressions.
    * \param[in] val/expr - Scalar value or expression.
    */
-#define MAKE_COMPOUND(OP)                                                 \
-  CSysVector& operator OP(ScalarType val) {                               \
-    CSYSVEC_PARFOR                                                        \
-    for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val;                  \
-    END_CSYSVEC_PARFOR                                                    \
-    return *this;                                                         \
-  }                                                                       \
-  template <class T>                                                      \
-  CSysVector& operator OP(const VecExpr::CVecExpr<T, ScalarType>& expr) { \
-    CSYSVEC_PARFOR                                                        \
-    for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP expr.derived()[i];    \
-    END_CSYSVEC_PARFOR                                                    \
-    return *this;                                                         \
+#define MAKE_COMPOUND(OP)                                                   \
+  CSysVector& operator OP(ScalarType val) {                                 \
+    CSYSVEC_PARFOR                                                          \
+    for (su2_index_t i = 0; i < nElm; ++i) vec_val[i] OP val;               \
+    END_CSYSVEC_PARFOR                                                      \
+    return *this;                                                           \
+  }                                                                         \
+  template <class T>                                                        \
+  CSysVector& operator OP(const VecExpr::CVecExpr<T, ScalarType>& expr) {   \
+    CSYSVEC_PARFOR                                                          \
+    for (su2_index_t i = 0; i < nElm; ++i) vec_val[i] OP expr.derived()[i]; \
+    END_CSYSVEC_PARFOR                                                      \
+    return *this;                                                           \
   }
   MAKE_COMPOUND(=)
   MAKE_COMPOUND(+=)
@@ -351,7 +352,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, 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
@@ -409,7 +410,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \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;
+    for (su2_index_t iVar = 0; iVar < nVar; iVar++) vec_val[iPoint * nVar + iVar] = 0.0;
   }
 
   /*!
@@ -422,9 +423,9 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
   template <class VectorType, bool Overwrite = true>
   FORCEINLINE void SetBlock(unsigned long 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];
     }
   }