Skip to content

Harden host-side sparse index width and allocation safety #2822

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
LwhJesse wants to merge 1 commit into
base: develop
Choose a base branch
from
Draft

Harden host-side sparse index width and allocation safety #2822

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion Common/include/code_config.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
12 changes: 6 additions & 6 deletions Common/include/geometry/CGeometry.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<su2uint> finiteVolumeLToUTranspMap; /*!< \brief FVM L-entry -> U-entry of its transpose. */
su2vector<su2uint> finiteVolumeUToLTranspMap; /*!< \brief FVM U-entry -> L-entry of its transpose. */
su2vector<su2uint> finiteElementLToUTranspMap; /*!< \brief FEM L-entry -> U-entry of its transpose. */
su2vector<su2uint> finiteElementUToLTranspMap; /*!< \brief FEM U-entry -> L-entry of its transpose. */
su2vector<su2_index_t> finiteVolumeLToUTranspMap; /*!< \brief FVM L-entry -> U-entry of its transpose. */
su2vector<su2_index_t> finiteVolumeUToLTranspMap; /*!< \brief FVM U-entry -> L-entry of its transpose. */
su2vector<su2_index_t> finiteElementLToUTranspMap; /*!< \brief FEM L-entry -> U-entry of its transpose. */
su2vector<su2_index_t> finiteElementUToLTranspMap; /*!< \brief FEM U-entry -> L-entry of its transpose. */

/*--- Edge and element colorings. ---*/

Expand Down Expand Up @@ -1892,15 +1892,15 @@ class CGeometry {
* \param[in] type - Finite volume or finite element.
* \return Reference to the l_to_u map.
*/
const su2vector<su2uint>& GetLToUTransposeSparsePatternMap(ConnectivityType type);
const su2vector<su2_index_t>& GetLToUTransposeSparsePatternMap(ConnectivityType type);

/*!
* \brief Get the bijective map from U-entry indices to L-entry indices of their transposes.
* \note Requires symmetric pattern. Builds both LU transpose maps if not already built.
* \param[in] type - Finite volume or finite element.
* \return Reference to the u_to_l map.
*/
const su2vector<su2uint>& GetUToLTransposeSparsePatternMap(ConnectivityType type);
const su2vector<su2_index_t>& GetUToLTransposeSparsePatternMap(ConnectivityType type);

/*!
* \brief Get the edge coloring.
Expand Down
94 changes: 64 additions & 30 deletions Common/include/linear_algebra/CPastixWrapper.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,8 @@

#include <pastix.h>
#include <spm.h>
#include <limits>
#include <string>
#include <vector>

using namespace std;
Expand All @@ -64,25 +66,30 @@ class CPastixWrapper {
vector<pastix_int_t> perm; /*!< \brief Ordering computed by PaStiX. */
vector<su2mixedfloat> workvec; /*!< \brief RHS vector which then becomes the solution. */

vector<unsigned long> csr_row_ptr; /*!< \brief Owned CSR row pointers (built from LDU). */
vector<unsigned long> csr_col_ind; /*!< \brief Owned CSR column indices (built from LDU). */
vector<su2_index_t> csr_row_ptr; /*!< \brief Owned CSR row pointers (built from LDU). */
vector<su2_index_t> 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<su2_index_t>::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. */
Expand All @@ -94,8 +101,30 @@ 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. */

static su2_index_t CheckedMul(su2_index_t lhs, su2_index_t rhs, const char* what) {
if (lhs != 0 && rhs > std::numeric_limits<su2_index_t>::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<su2_index_t>::max() - lhs) {
SU2_MPI::Error(std::string("Overflow while computing ") + what + ".", CURRENT_FUNCTION);
}
return lhs + rhs;
}

template <class TargetType>
static TargetType CheckedCast(su2_index_t value, const char* what) {
if (value > static_cast<su2_index_t>(std::numeric_limits<TargetType>::max())) {
SU2_MPI::Error(std::string("Overflow while converting ") + what + ".", CURRENT_FUNCTION);
}
return static_cast<TargetType>(value);
}

/*!
* \brief Run the "clean" task, releases all memory, leaves object in unusable state.
Expand Down Expand Up @@ -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;
Expand All @@ -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<unsigned long>(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<size_t>(CheckedAdd(nPointDomain, 1, "PaStiX CSR row pointer size"), "PaStiX CSR row pointer size"));
csr_col_ind.reserve(CheckedCast<size_t>(nnz_domain, "PaStiX CSR column index size"));
for (su2_index_t i = 0; i < nPointDomain; ++i) {
csr_row_ptr[i] = static_cast<su2_index_t>(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<unsigned long>(csr_col_ind.size());
csr_row_ptr[nPointDomain] = static_cast<su2_index_t>(csr_col_ind.size());
issetup = true;
}

Expand Down Expand Up @@ -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<pastix_int_t>(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
33 changes: 17 additions & 16 deletions Common/include/linear_algebra/CSysMatrix.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <typename IndexType>
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<su2_index_t> mat; /*!< \brief Host matrix (values owned via aligned_alloc; pattern from geometry). */
LDU<su2uint> gpu; /*!< \brief Device matrix (all pointers to GPU memory). */
LDU<su2_index_t> ilu; /*!< \brief ILU factorization, host (values owned; pattern from geometry). */
LDU<su2uint> 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). ---*/
Expand Down Expand Up @@ -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. */

Expand Down
Loading
Loading