From 53bacf193f4240d62492378c27319f517779bd4d Mon Sep 17 00:00:00 2001 From: LwhJesse <256257451+LwhJesse@users.noreply.github.com> Date: Mon, 18 May 2026 23:13:56 +0800 Subject: [PATCH 01/10] Cache CUDA SpMV cuSPARSE resources --- Common/include/linear_algebra/CSysMatrix.hpp | 11 +++ Common/src/linear_algebra/CSysMatrix.cpp | 6 ++ Common/src/linear_algebra/CSysMatrixGPU.cu | 79 +++++++++++++++----- 3 files changed, 76 insertions(+), 20 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index ecb26959a06..1b8a7565ce5 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -108,6 +108,16 @@ struct CSysMatrixComms { MPI_QUANTITIES commType = MPI_QUANTITIES::SOLUTION_MATRIX); }; +/*! + * \brief Opaque storage for CUDA/cuSPARSE resources used by the GPU SpMV path. + */ +struct CudaSpMVResources; + +/*! + * \brief Release cached CUDA/cuSPARSE resources used by the GPU SpMV path. + */ +void ReleaseCudaSpMVResources(CudaSpMVResources*& resources); + /*! * \class CSysMatrix * \ingroup SpLinSys @@ -150,6 +160,7 @@ class CSysMatrix { 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. */ + mutable CudaSpMVResources* spmv_resources = nullptr; /*!< \brief Cached cuSPARSE resources for GPU SpMV. */ ScalarType* ILU_matrix; /*!< \brief Entries of the ILU sparse matrix. */ unsigned long nnz_ilu; /*!< \brief Number of possible nonzero entries in the matrix (ILU). */ diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 00ef51c2023..30c5d52e014 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -67,6 +67,8 @@ template CSysMatrix::~CSysMatrix() { SU2_ZONE_SCOPED + ReleaseCudaSpMVResources(spmv_resources); + delete[] omp_partitions; MemoryAllocation::aligned_free(ILU_matrix); MemoryAllocation::aligned_free(matrix); @@ -86,6 +88,10 @@ CSysMatrix::~CSysMatrix() { #endif } +#ifndef HAVE_CUDA +void ReleaseCudaSpMVResources(CudaSpMVResources*& resources) { resources = nullptr; } +#endif + template void CSysMatrix::Initialize(unsigned long npoint, unsigned long npointdomain, unsigned short nvar, unsigned short neqn, bool EdgeConnect, CGeometry* geometry, diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index a3f1a77ca40..af31b9d8ec1 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -55,6 +55,38 @@ inline cusparseIndexType_t GetCusparseIndexType() { } } +struct CudaSpMVResources { + cusparseHandle_t handle = nullptr; + cusparseConstSpMatDescr_t mat = nullptr; + size_t buffer_size = 0; + void* buffer = nullptr; +}; + +void ReleaseCudaSpMVResources(CudaSpMVResources*& resources) { + if (resources == nullptr) { + return; + } + + if (resources->buffer != nullptr) { + gpuErrChk(cudaFree(resources->buffer)); + resources->buffer = nullptr; + resources->buffer_size = 0; + } + + if (resources->mat != nullptr) { + cusparseErrChk(cusparseDestroySpMat(resources->mat)); + resources->mat = nullptr; + } + + if (resources->handle != nullptr) { + cusparseErrChk(cusparseDestroy(resources->handle)); + resources->handle = nullptr; + } + + delete resources; + resources = nullptr; +} + template constexpr cudaDataType GetCudaDataType() { if constexpr (std::is_same::value) { @@ -100,40 +132,47 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector const ScalarType alpha = 1.0; const ScalarType beta = 0.0; - cusparseHandle_t handle = nullptr; - cusparseConstSpMatDescr_t matA = nullptr; - cusparseDnVecDescr_t vecX = nullptr; - cusparseDnVecDescr_t vecY = nullptr; + if (spmv_resources == nullptr) { + spmv_resources = new CudaSpMVResources; + + cusparseErrChk(cusparseCreate(&spmv_resources->handle)); - cusparseErrChk(cusparseCreate(&handle)); + cusparseErrChk(cusparseCreateConstBsr(&spmv_resources->mat, brows, bcols, bnnz, blockSize, blockSize, d_row_ptr, + d_col_ind, d_matrix, indexType, indexType, CUSPARSE_INDEX_BASE_ZERO, + valueType, CUSPARSE_ORDER_ROW)); + } - cusparseErrChk(cusparseCreateConstBsr(&matA, brows, bcols, bnnz, blockSize, blockSize, d_row_ptr, d_col_ind, d_matrix, - indexType, indexType, CUSPARSE_INDEX_BASE_ZERO, valueType, CUSPARSE_ORDER_ROW)); + cusparseDnVecDescr_t vecX = nullptr; + cusparseDnVecDescr_t vecY = nullptr; cusparseErrChk(cusparseCreateDnVec(&vecX, xSize, d_vec, valueType)); cusparseErrChk(cusparseCreateDnVec(&vecY, ySize, d_prod, valueType)); - size_t bufferSize = 0; - void* dBuffer = nullptr; + size_t required_buffer_size = 0; - cusparseErrChk(cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, - valueType, CUSPARSE_SPMV_BSR_ALG1, &bufferSize)); + cusparseErrChk(cusparseSpMV_bufferSize(spmv_resources->handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, + spmv_resources->mat, vecX, &beta, vecY, valueType, CUSPARSE_SPMV_BSR_ALG1, + &required_buffer_size)); - if (bufferSize > 0) { - gpuErrChk(cudaMalloc(&dBuffer, bufferSize)); - } + if (required_buffer_size > spmv_resources->buffer_size) { + if (spmv_resources->buffer != nullptr) { + gpuErrChk(cudaFree(spmv_resources->buffer)); + } - cusparseErrChk(cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta, vecY, valueType, - CUSPARSE_SPMV_BSR_ALG1, dBuffer)); + if (required_buffer_size > 0) { + gpuErrChk(cudaMalloc(&spmv_resources->buffer, required_buffer_size)); + } else { + spmv_resources->buffer = nullptr; + } - if (dBuffer != nullptr) { - gpuErrChk(cudaFree(dBuffer)); + spmv_resources->buffer_size = required_buffer_size; } + cusparseErrChk(cusparseSpMV(spmv_resources->handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, spmv_resources->mat, + vecX, &beta, vecY, valueType, CUSPARSE_SPMV_BSR_ALG1, spmv_resources->buffer)); + cusparseErrChk(cusparseDestroyDnVec(vecY)); cusparseErrChk(cusparseDestroyDnVec(vecX)); - cusparseErrChk(cusparseDestroySpMat(matA)); - cusparseErrChk(cusparseDestroy(handle)); prod.DtHTransfer(); } From 08fde80e1e8caf04f966e578e6bdb28a74ec1ec5 Mon Sep 17 00:00:00 2001 From: LwhJesse <256257451+LwhJesse@users.noreply.github.com> Date: Mon, 1 Jun 2026 22:18:21 +0800 Subject: [PATCH 02/10] Add CUDA FGMRES and Jacobi scaffolding --- .../linear_algebra/CPreconditioner.hpp | 15 ++++ Common/include/linear_algebra/CSysMatrix.hpp | 16 ++++ Common/include/linear_algebra/CSysSolve.hpp | 10 +++ Common/include/linear_algebra/CSysVector.hpp | 39 ++++++++++ Common/include/linear_algebra/GPUComms.cuh | 7 ++ Common/include/option_structure.hpp | 2 + Common/src/linear_algebra/CSysMatrix.cpp | 32 ++++++++ Common/src/linear_algebra/CSysMatrixGPU.cu | 17 +++-- .../linear_algebra/CSysPreconditionerGPU.cu | 70 ++++++++++++++++++ Common/src/linear_algebra/CSysSolve.cpp | 11 +++ Common/src/linear_algebra/CSysSolveGPU.cu | 68 +++++++++++++++++ Common/src/linear_algebra/CSysVectorGPU.cu | 73 ++++++++++++++++++- Common/src/linear_algebra/meson.build | 2 +- Common/src/meson.build | 30 +++++++- meson.build | 7 +- 15 files changed, 386 insertions(+), 13 deletions(-) create mode 100644 Common/src/linear_algebra/CSysPreconditionerGPU.cu create mode 100644 Common/src/linear_algebra/CSysSolveGPU.cu diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp index e4fc7cf159f..4eff5f947fb 100644 --- a/Common/include/linear_algebra/CPreconditioner.hpp +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -77,6 +77,18 @@ class CPreconditioner { template CPreconditioner::~CPreconditioner() {} +/*! + * \class CIdentityPreconditioner + * \brief No-op preconditioner used when Krylov solvers run without preconditioning. + */ +template +class CIdentityPreconditioner final : public CPreconditioner { + public: + inline void operator()(const CSysVector& u, CSysVector& v) const override { v = u; } + + inline bool IsIdentity() const override { return true; } +}; + /*! * \class CJacobiPreconditioner * \brief Specialization of preconditioner that uses CSysMatrix class. @@ -332,6 +344,9 @@ CPreconditioner* CPreconditioner::Create(ENUM_LINEAR_SOL CPreconditioner* prec = nullptr; switch (kind) { + case NO_PRECONDITIONER: + prec = new CIdentityPreconditioner(); + break; case JACOBI: prec = new CJacobiPreconditioner(jacobian, geometry, config); break; diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 1b8a7565ce5..6368670587e 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -158,6 +158,7 @@ class CSysMatrix { 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(). */ + ScalarType* d_invM = nullptr; /*!< \brief Device inverse diagonal blocks for the Jacobi preconditioner. */ 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. */ mutable CudaSpMVResources* spmv_resources = nullptr; /*!< \brief Cached cuSPARSE resources for GPU SpMV. */ @@ -935,6 +936,13 @@ class CSysMatrix { */ void BuildJacobiPreconditioner(); + /*! + * \brief Build the Jacobi preconditioner on the GPU/device side. + * \note This helper is intended as the implementation hook for GPU-resident Krylov solvers. + * The actual implementation belongs in CSysMatrixGPU.cu. + */ + void BuildJacobiPreconditionerGPU(); + /*! * \brief Multiply CSysVector by the preconditioner * \param[in] vec - CSysVector to be multiplied by the preconditioner. @@ -945,6 +953,14 @@ class CSysMatrix { void ComputeJacobiPreconditioner(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const; + /*! + * \brief Apply the Jacobi preconditioner on the GPU/device side. + * \note This helper is intended as the implementation hook for GPU-resident Krylov solvers. + * The actual implementation belongs in CSysMatrixGPU.cu. + */ + void ComputeJacobiPreconditionerGPU(const CSysVector& vec, CSysVector& prod, + CGeometry* geometry, const CConfig* config) const; + /*! * \brief Build the ILU preconditioner. */ diff --git a/Common/include/linear_algebra/CSysSolve.hpp b/Common/include/linear_algebra/CSysSolve.hpp index 2ea3cbf7df3..67033992372 100644 --- a/Common/include/linear_algebra/CSysSolve.hpp +++ b/Common/include/linear_algebra/CSysSolve.hpp @@ -304,6 +304,16 @@ class CSysSolve { ScalarType& residual, bool monitoring, const CConfig* config, FgcrodrMode mode, unsigned long custom_m) const; + /*! + * \brief CUDA/GPU implementation of the Flexible GMRES solver. + * \note This helper exists so the public solver interface can remain unchanged while the + * implementation is dispatched internally based on the runtime CUDA setting. + * The actual implementation lives in CSysSolveGPU.cu. + */ + unsigned long FGMRES_LinSolver_GPU(const VectorType& b, VectorType& x, const ProductType& mat_vec, + const PrecondType& precond, ScalarType tol, unsigned long m, ScalarType& residual, + bool monitoring, const CConfig* config) const; + /*! * \brief Creates the inner solver for nested preconditioning if the settings allow it. * \returns True if the inner solver can be used. diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 1498b549bbb..decffe26d29 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -240,6 +240,45 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ void GPUSetVal(ScalarType val, bool trigger = true) const; + /*! + * \brief Copy another vector into this vector on the device. + * \note This is an explicit GPU helper for Krylov solver implementations that keep the + * iteration state resident on the device. The implementation belongs in + * CSysVectorGPU.cu. + * \param[in] src - Source vector. + */ + void GPUCopy(const CSysVector& src) const; + + /*! + * \brief Scale this vector on the device. + * \note Explicit GPU helper for solver-side vector operations. + * \param[in] alpha - Scalar multiplier. + */ + void GPUScale(ScalarType alpha) const; + + /*! + * \brief Perform the AXPY operation on the device: this := this + alpha * x. + * \note Explicit GPU helper for solver-side vector operations. + * \param[in] alpha - Scalar multiplier. + * \param[in] x - Input vector. + */ + void GPUAxpy(ScalarType alpha, const CSysVector& x) const; + + /*! + * \brief Dot product between this vector and another vector on the device. + * \note Explicit GPU helper for solver-side reductions. + * \param[in] other - Input vector. + * \return Dot product result. + */ + ScalarType GPUDot(const CSysVector& other) const; + + /*! + * \brief L2 norm of this vector on the device. + * \note Explicit GPU helper for solver-side reductions. + * \return L2 norm result. + */ + ScalarType GPUNorm() const; + /*! * \brief return device pointer that points to the CSysVector values in GPU memory */ diff --git a/Common/include/linear_algebra/GPUComms.cuh b/Common/include/linear_algebra/GPUComms.cuh index 13854381877..eb727477f21 100644 --- a/Common/include/linear_algebra/GPUComms.cuh +++ b/Common/include/linear_algebra/GPUComms.cuh @@ -25,6 +25,11 @@ * License along with SU2. If not, see . */ +#pragma once + +#ifndef SU2_COMMON_LINEAR_ALGEBRA_GPUCOMMS_CUH +#define SU2_COMMON_LINEAR_ALGEBRA_GPUCOMMS_CUH + #include #include @@ -51,3 +56,5 @@ inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=t } #define gpuErrChk(ans) { gpuAssert((ans), __FILE__, __LINE__); } + +#endif // SU2_COMMON_LINEAR_ALGEBRA_GPUCOMMS_CUH diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 4c521c29239..65ba6dfa1ab 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2517,6 +2517,7 @@ enum ENUM_LINEAR_SOLVER_PREC { PASTIX_ILU=10, /*!< \brief PaStiX ILU(k) preconditioner. */ PASTIX_LU_P, /*!< \brief PaStiX LU as preconditioner. */ PASTIX_LDLT_P, /*!< \brief PaStiX LDLT as preconditioner. */ + NO_PRECONDITIONER=20, /*!< \brief No preconditioner. */ }; static const MapType Linear_Solver_Prec_Map = { MakePair("JACOBI", JACOBI) @@ -2526,6 +2527,7 @@ static const MapType Linear_Solver_Prec_Ma MakePair("PASTIX_ILU", PASTIX_ILU) MakePair("PASTIX_LU", PASTIX_LU_P) MakePair("PASTIX_LDLT", PASTIX_LDLT_P) + MakePair("NONE", NO_PRECONDITIONER) }; /*! diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 30c5d52e014..fafa093734e 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -54,6 +54,7 @@ CSysMatrix::CSysMatrix() : rank(SU2_MPI::GetRank()), size(SU2_MPI::G col_ind_ilu = nullptr; invM = nullptr; + d_invM = nullptr; #ifdef USE_MKL MatrixMatrixProductJitter = nullptr; @@ -78,6 +79,7 @@ CSysMatrix::~CSysMatrix() { GPUMemoryAllocation::gpu_free(d_matrix); GPUMemoryAllocation::gpu_free(d_row_ptr); GPUMemoryAllocation::gpu_free(d_col_ind); + GPUMemoryAllocation::gpu_free(d_invM); } #ifdef USE_MKL @@ -202,6 +204,10 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi if (diag_needed) allocAndInit(invM, nPointDomain * nVar * nEqn); + if (useCuda && diag_needed) { + d_invM = GPUMemoryAllocation::gpu_alloc(nPointDomain * nVar * nEqn * sizeof(ScalarType)); + } + /*--- Thread parallel initialization. ---*/ int num_threads = omp_get_max_threads(); @@ -678,6 +684,19 @@ void CSysMatrix::MatrixVectorProduct(const CSysVector& v template void CSysMatrix::BuildJacobiPreconditioner() { SU2_ZONE_SCOPED + + if (useCuda) { +#ifdef HAVE_CUDA + BuildJacobiPreconditionerGPU(); + return; +#else + SU2_MPI::Error( + "\nError in building Jacobi preconditioner\nENABLE_CUDA is set to YES\nPlease compile with CUDA options " + "enabled in Meson to access GPU Functions", + CURRENT_FUNCTION); +#endif + } + /*--- Build Jacobi preconditioner (M = D), compute and store the inverses of the diagonal blocks. ---*/ SU2_OMP_FOR_DYN(omp_heavy_size) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) @@ -690,6 +709,19 @@ void CSysMatrix::ComputeJacobiPreconditioner(const CSysVector& prod, CGeometry* geometry, const CConfig* config) const { SU2_ZONE_SCOPED + + if (config->GetCUDA()) { +#ifdef HAVE_CUDA + ComputeJacobiPreconditionerGPU(vec, prod, geometry, config); + return; +#else + SU2_MPI::Error( + "\nError in applying Jacobi preconditioner\nENABLE_CUDA is set to YES\nPlease compile with CUDA options " + "enabled in Meson to access GPU Functions", + CURRENT_FUNCTION); +#endif + } + /*--- Apply Jacobi preconditioner, y = D^{-1} * x, the inverse of the diagonal is already known. ---*/ SU2_OMP_BARRIER SU2_OMP_FOR_DYN(omp_heavy_size) diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index af31b9d8ec1..b8bb16ce5a9 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -115,8 +115,6 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector ScalarType* d_vec = vec.GetDevicePointer(); ScalarType* d_prod = prod.GetDevicePointer(); - vec.HtDTransfer(); - const auto indexType = GetCusparseIndexType(); const auto valueType = GetCudaDataType(); @@ -173,8 +171,15 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector cusparseErrChk(cusparseDestroyDnVec(vecY)); cusparseErrChk(cusparseDestroyDnVec(vecX)); - - prod.DtHTransfer(); } - -template class CSysMatrix; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits. +template void CSysMatrix::HtDTransfer(bool trigger) const; +template void CSysMatrix::GPUMatrixVectorProduct(const CSysVector& vec, + CSysVector& prod, CGeometry* geometry, + const CConfig* config) const; + +#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION) +template void CSysMatrix::HtDTransfer(bool trigger) const; +template void CSysMatrix::GPUMatrixVectorProduct(const CSysVector& vec, + CSysVector& prod, CGeometry* geometry, + const CConfig* config) const; +#endif diff --git a/Common/src/linear_algebra/CSysPreconditionerGPU.cu b/Common/src/linear_algebra/CSysPreconditionerGPU.cu new file mode 100644 index 00000000000..8bef1c68be9 --- /dev/null +++ b/Common/src/linear_algebra/CSysPreconditionerGPU.cu @@ -0,0 +1,70 @@ +/*! + * \file CSysPreconditionerGPU.cu + * \brief CUDA/GPU skeleton implementations for matrix-based preconditioners. + * \author Jesse Li + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/linear_algebra/CSysMatrix.hpp" +#include "../../include/linear_algebra/GPUComms.cuh" + +template +void CSysMatrix::BuildJacobiPreconditionerGPU() { + /*--- Implementation area for building the GPU/device Jacobi preconditioner. + * This skeleton intentionally leaves the algorithm unimplemented so the + * surrounding dispatch and file structure can be reviewed independently of + * the CUDA preconditioner details. ---*/ + SU2_MPI::Error("CSysMatrix::BuildJacobiPreconditionerGPU skeleton reached without an implementation.", + CURRENT_FUNCTION); +} + +template +void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector& vec, + CSysVector& prod, CGeometry* geometry, + const CConfig* config) const { + (void)vec; + (void)prod; + (void)geometry; + (void)config; + + /*--- Implementation area for applying the GPU/device Jacobi preconditioner. + * This skeleton intentionally leaves the algorithm unimplemented so the + * surrounding dispatch and file structure can be reviewed independently of + * the CUDA preconditioner details. ---*/ + SU2_MPI::Error("CSysMatrix::ComputeJacobiPreconditionerGPU skeleton reached without an implementation.", + CURRENT_FUNCTION); +} + +template void CSysMatrix::BuildJacobiPreconditionerGPU(); +template void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector& vec, + CSysVector& prod, + CGeometry* geometry, + const CConfig* config) const; + +#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION) +template void CSysMatrix::BuildJacobiPreconditionerGPU(); +template void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector& vec, + CSysVector& prod, + CGeometry* geometry, + const CConfig* config) const; +#endif diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index cd100f311ab..3eba623f051 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -421,6 +421,17 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVectorGetCUDA()) { +#ifdef HAVE_CUDA + return FGMRES_LinSolver_GPU(b, x, mat_vec, precond, tol, m, residual, monitoring, config); +#else + SU2_MPI::Error( + "\nError in launching FGMRES solver\nENABLE_CUDA is set to YES\nPlease compile with CUDA options enabled " + "in Meson to access GPU Functions", + CURRENT_FUNCTION); +#endif + } + const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE); const bool flexible = !precond.IsIdentity(); /*--- If we call the solver outside of a parallel region, but the number of threads allows, diff --git a/Common/src/linear_algebra/CSysSolveGPU.cu b/Common/src/linear_algebra/CSysSolveGPU.cu new file mode 100644 index 00000000000..3a23e5eda61 --- /dev/null +++ b/Common/src/linear_algebra/CSysSolveGPU.cu @@ -0,0 +1,68 @@ +/*! + * \file CSysSolveGPU.cu + * \brief CUDA/GPU skeleton implementations for Krylov solvers. + * \author Jesse Li + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/linear_algebra/CSysSolve.hpp" +#include "../../include/linear_algebra/GPUComms.cuh" + +template +unsigned long CSysSolve::FGMRES_LinSolver_GPU(const VectorType& b, VectorType& x, + const ProductType& mat_vec, const PrecondType& precond, + ScalarType tol, unsigned long m, ScalarType& residual, + bool monitoring, const CConfig* config) const { + SU2_ZONE_SCOPED + + (void)b; + (void)x; + (void)mat_vec; + (void)precond; + (void)tol; + (void)m; + (void)residual; + (void)monitoring; + (void)config; + + /*--- Implementation area for the GPU-resident FGMRES solve path. + * This skeleton intentionally leaves the algorithm unimplemented so the + * surrounding dispatch, build integration, and file structure can be + * reviewed independently of the CUDA solver details. ---*/ + SU2_MPI::Error("FGMRES GPU solver skeleton reached without an implementation.", CURRENT_FUNCTION); + return 0; +} + +/*--- Explicit instantiations for the GPU solver helper. + * Keep these aligned with the scalar types instantiated for CSysSolve. ---*/ +template unsigned long CSysSolve::FGMRES_LinSolver_GPU( + const CSysVector& b, CSysVector& x, + const CMatrixVectorProduct& mat_vec, const CPreconditioner& precond, + su2mixedfloat tol, unsigned long m, su2mixedfloat& residual, bool monitoring, const CConfig* config) const; + +#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION) +template unsigned long CSysSolve::FGMRES_LinSolver_GPU( + const CSysVector& b, CSysVector& x, + const CMatrixVectorProduct& mat_vec, const CPreconditioner& precond, + passivedouble tol, unsigned long m, passivedouble& residual, bool monitoring, const CConfig* config) const; +#endif diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 94ec17bb88f..c5964f5de98 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -27,7 +27,6 @@ #include "../../include/linear_algebra/CSysVector.hpp" #include "../../include/linear_algebra/GPUComms.cuh" - template void CSysVector::HtDTransfer(bool trigger) const { @@ -46,4 +45,74 @@ void CSysVector::GPUSetVal(ScalarType val, bool trigger) const if(trigger) gpuErrChk(cudaMemset((void*)(d_vec_val), val, (sizeof(ScalarType)*nElm))); } -template class CSysVector; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits. +template +void CSysVector::GPUCopy(const CSysVector& src) const { + (void)src; + + /*--- Implementation area for GPU-to-GPU vector copy used by device-resident Krylov solvers. ---*/ + SU2_MPI::Error("CSysVector::GPUCopy skeleton reached without an implementation.", CURRENT_FUNCTION); +} + +template +void CSysVector::GPUScale(ScalarType alpha) const { + (void)alpha; + + /*--- Implementation area for GPU vector scaling used by device-resident Krylov solvers. ---*/ + SU2_MPI::Error("CSysVector::GPUScale skeleton reached without an implementation.", CURRENT_FUNCTION); +} + +template +void CSysVector::GPUAxpy(ScalarType alpha, const CSysVector& x) const { + (void)alpha; + (void)x; + + /*--- Implementation area for GPU AXPY used by device-resident Krylov solvers. ---*/ + SU2_MPI::Error("CSysVector::GPUAxpy skeleton reached without an implementation.", CURRENT_FUNCTION); +} + +template +ScalarType CSysVector::GPUDot(const CSysVector& other) const { + (void)other; + + /*--- Implementation area for GPU dot products used by device-resident Krylov solvers. ---*/ + SU2_MPI::Error("CSysVector::GPUDot skeleton reached without an implementation.", CURRENT_FUNCTION); + return ScalarType(0); +} + +template +ScalarType CSysVector::GPUNorm() const { + /*--- Implementation area for GPU norms used by device-resident Krylov solvers. ---*/ + SU2_MPI::Error("CSysVector::GPUNorm skeleton reached without an implementation.", CURRENT_FUNCTION); + return ScalarType(0); +} + +template void CSysVector::HtDTransfer(bool trigger) const; +template void CSysVector::DtHTransfer(bool trigger) const; +template void CSysVector::GPUSetVal(su2mixedfloat val, bool trigger) const; +template void CSysVector::GPUCopy(const CSysVector& src) const; +template void CSysVector::GPUScale(su2mixedfloat alpha) const; +template void CSysVector::GPUAxpy(su2mixedfloat alpha, const CSysVector& x) const; +template su2mixedfloat CSysVector::GPUDot(const CSysVector& other) const; +template su2mixedfloat CSysVector::GPUNorm() const; + +#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION) +template void CSysVector::HtDTransfer(bool trigger) const; +template void CSysVector::DtHTransfer(bool trigger) const; +template void CSysVector::GPUSetVal(passivedouble val, bool trigger) const; +template void CSysVector::GPUCopy(const CSysVector& src) const; +template void CSysVector::GPUScale(passivedouble alpha) const; +template void CSysVector::GPUAxpy(passivedouble alpha, const CSysVector& x) const; +template passivedouble CSysVector::GPUDot(const CSysVector& other) const; +template passivedouble CSysVector::GPUNorm() const; +#endif + +#ifdef CODI_REVERSE_TYPE +template void CSysVector::HtDTransfer(bool trigger) const; +template void CSysVector::DtHTransfer(bool trigger) const; +template void CSysVector::GPUSetVal(su2double val, bool trigger) const; +template void CSysVector::GPUCopy(const CSysVector& src) const; +template void CSysVector::GPUScale(su2double alpha) const; +template void CSysVector::GPUAxpy(su2double alpha, const CSysVector& x) const; +template su2double CSysVector::GPUDot(const CSysVector& other) const; +template su2double CSysVector::GPUNorm() const; +#endif diff --git a/Common/src/linear_algebra/meson.build b/Common/src/linear_algebra/meson.build index 7b880b29c1e..29bb63cd680 100644 --- a/Common/src/linear_algebra/meson.build +++ b/Common/src/linear_algebra/meson.build @@ -6,5 +6,5 @@ common_src += files(['CSysSolve_b.cpp', 'blas_structure.cpp']) if get_option('enable-cuda') - common_src += files(['CSysMatrixGPU.cu', 'CSysVectorGPU.cu',]) + common_src += files(['CSysMatrixGPU.cu', 'CSysVectorGPU.cu', 'CSysSolveGPU.cu', 'CSysPreconditionerGPU.cu',]) endif diff --git a/Common/src/meson.build b/Common/src/meson.build index f385c4a32ed..710e764303f 100644 --- a/Common/src/meson.build +++ b/Common/src/meson.build @@ -6,6 +6,27 @@ common_src =files(['graph_coloring_structure.cpp', '../include/parallelization/mpi_structure.cpp', '../include/parallelization/omp_structure.cpp']) +common_cuda_cpp_args = [] +foreach arg : su2_cpp_args + if arg.startswith('-D') or arg.startswith('-U') + common_cuda_cpp_args += arg + endif +endforeach + +common_cuda_rev_args = common_cuda_cpp_args +foreach arg : codi_rev_args + if arg.startswith('-D') or arg.startswith('-U') + common_cuda_rev_args += arg + endif +endforeach + +common_cuda_for_args = common_cuda_cpp_args +foreach arg : codi_for_args + if arg.startswith('-D') or arg.startswith('-U') + common_cuda_for_args += arg + endif +endforeach + subdir('linear_algebra') subdir('toolboxes') subdir('geometry') @@ -24,7 +45,8 @@ if get_option('enable-normal') common_src, install : false, dependencies : su2_deps, - cpp_args: [default_warning_flags, su2_cpp_args]) + cpp_args: [default_warning_flags, su2_cpp_args], + cuda_args: common_cuda_cpp_args) common_dep = declare_dependency(link_with: common, include_directories : common_include) @@ -36,7 +58,8 @@ if get_option('enable-autodiff') common_src, install : false, dependencies : [su2_deps, codi_dep], - cpp_args: [default_warning_flags, su2_cpp_args, codi_rev_args]) + cpp_args: [default_warning_flags, su2_cpp_args, codi_rev_args], + cuda_args: common_cuda_rev_args) commonAD_dep = declare_dependency(link_with: commonAD, include_directories : common_include) @@ -49,7 +72,8 @@ if get_option('enable-directdiff') common_src, install : false, dependencies : [su2_deps, codi_dep], - cpp_args: [default_warning_flags, su2_cpp_args, codi_for_args]) + cpp_args: [default_warning_flags, su2_cpp_args, codi_for_args], + cuda_args: common_cuda_for_args) commonDD_dep = declare_dependency(link_with: commonDD, include_directories : common_include) diff --git a/meson.build b/meson.build index d43e2b031d8..43786c100e9 100644 --- a/meson.build +++ b/meson.build @@ -20,13 +20,18 @@ python = pymod.find_installation() if get_option('enable-cuda') add_languages('cuda') add_global_arguments('-arch=sm_86', language : 'cuda') - cuda_deps = [meson.get_compiler('cuda').find_library('cusparse', required : true)] + cuda_deps = [ + meson.get_compiler('cuda').find_library('cusparse', required : true), + meson.get_compiler('cuda').find_library('cublas', required : true), + ] else cuda_deps = [] endif su2_cpp_args = [] su2_deps = [declare_dependency(include_directories: 'externals/CLI11')] + cuda_deps +codi_rev_args = [] +codi_for_args = [] default_warning_flags = [] if build_machine.system() != 'windows' From fde2c145cf35a7945c064a67a523ca63b8ac7bb5 Mon Sep 17 00:00:00 2001 From: LwhJesse <256257451+LwhJesse@users.noreply.github.com> Date: Mon, 1 Jun 2026 22:18:31 +0800 Subject: [PATCH 03/10] Implement CUDA vector primitives --- Common/src/linear_algebra/CSysVectorGPU.cu | 139 ++++++++++++++++++--- 1 file changed, 122 insertions(+), 17 deletions(-) diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index c5964f5de98..d69e9fcf234 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -27,6 +27,87 @@ #include "../../include/linear_algebra/CSysVector.hpp" #include "../../include/linear_algebra/GPUComms.cuh" +#include +#include + +namespace { + +constexpr unsigned GPU_VEC_OP_BLOCK_SIZE = 256; + +template +__global__ void GPUCopyKernel(const ScalarType* src, ScalarType* dst, unsigned long nElm) { + const unsigned long idx = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (idx < nElm) dst[idx] = src[idx]; +} + +template +__global__ void GPUScaleKernel(ScalarType alpha, ScalarType* vec, unsigned long nElm) { + const unsigned long idx = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (idx < nElm) vec[idx] *= alpha; +} + +template +__global__ void GPUAxpyKernel(ScalarType alpha, const ScalarType* x, ScalarType* y, unsigned long nElm) { + const unsigned long idx = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (idx < nElm) y[idx] += alpha * x[idx]; +} + +inline dim3 MakeVectorGrid(unsigned long nElm) { + const auto grid = (nElm + GPU_VEC_OP_BLOCK_SIZE - 1) / GPU_VEC_OP_BLOCK_SIZE; + return dim3(static_cast(grid), 1, 1); +} + +thread_local cublasHandle_t active_solver_blas_handle = nullptr; +thread_local unsigned active_solver_blas_depth = 0; + +cublasHandle_t GetBlasHandle(const char* creation_error, bool& owns_handle) { + if (active_solver_blas_handle != nullptr) { + owns_handle = false; + return active_solver_blas_handle; + } + + cublasHandle_t handle = nullptr; + owns_handle = true; + if (cublasCreate(&handle) != CUBLAS_STATUS_SUCCESS) { + SU2_MPI::Error(creation_error, CURRENT_FUNCTION); + } + return handle; +} + +void ReleaseBlasHandle(cublasHandle_t handle, bool owns_handle, const char* destruction_error) { + if (owns_handle && handle != nullptr && cublasDestroy(handle) != CUBLAS_STATUS_SUCCESS) { + SU2_MPI::Error(destruction_error, CURRENT_FUNCTION); + } +} + +} // namespace + +void SU2_GPU_BeginSolverBLASContext() { + if (active_solver_blas_depth == 0) { + cublasHandle_t handle = nullptr; + if (cublasCreate(&handle) != CUBLAS_STATUS_SUCCESS) { + SU2_MPI::Error("cuBLAS handle creation failed for the GPU linear solver context.", CURRENT_FUNCTION); + } + active_solver_blas_handle = handle; + } + ++active_solver_blas_depth; +} + +void SU2_GPU_EndSolverBLASContext() { + if (active_solver_blas_depth == 0) { + SU2_MPI::Error("GPU linear solver BLAS context ended without a matching begin.", CURRENT_FUNCTION); + } + + --active_solver_blas_depth; + if (active_solver_blas_depth == 0) { + auto status = cublasDestroy(active_solver_blas_handle); + active_solver_blas_handle = nullptr; + if (status != CUBLAS_STATUS_SUCCESS) { + SU2_MPI::Error("cuBLAS handle destruction failed for the GPU linear solver context.", CURRENT_FUNCTION); + } + } +} + template void CSysVector::HtDTransfer(bool trigger) const { @@ -47,43 +128,67 @@ void CSysVector::GPUSetVal(ScalarType val, bool trigger) const template void CSysVector::GPUCopy(const CSysVector& src) const { - (void)src; + if (nElm == 0) return; - /*--- Implementation area for GPU-to-GPU vector copy used by device-resident Krylov solvers. ---*/ - SU2_MPI::Error("CSysVector::GPUCopy skeleton reached without an implementation.", CURRENT_FUNCTION); + GPUCopyKernel<<>>(src.GetDevicePointer(), GetDevicePointer(), nElm); + gpuErrChk(cudaPeekAtLastError()); } template void CSysVector::GPUScale(ScalarType alpha) const { - (void)alpha; + if (nElm == 0) return; - /*--- Implementation area for GPU vector scaling used by device-resident Krylov solvers. ---*/ - SU2_MPI::Error("CSysVector::GPUScale skeleton reached without an implementation.", CURRENT_FUNCTION); + GPUScaleKernel<<>>(alpha, GetDevicePointer(), nElm); + gpuErrChk(cudaPeekAtLastError()); } template void CSysVector::GPUAxpy(ScalarType alpha, const CSysVector& x) const { - (void)alpha; - (void)x; + if (nElm == 0) return; - /*--- Implementation area for GPU AXPY used by device-resident Krylov solvers. ---*/ - SU2_MPI::Error("CSysVector::GPUAxpy skeleton reached without an implementation.", CURRENT_FUNCTION); + GPUAxpyKernel<<>>(alpha, x.GetDevicePointer(), GetDevicePointer(), + nElm); + gpuErrChk(cudaPeekAtLastError()); } template ScalarType CSysVector::GPUDot(const CSysVector& other) const { - (void)other; + bool owns_handle = false; + cublasHandle_t handle = GetBlasHandle("cuBLAS handle creation failed in CSysVector::GPUDot.", owns_handle); + cublasStatus_t status = CUBLAS_STATUS_SUCCESS; + + ScalarType local_dot = ScalarType(0); + + if constexpr (std::is_same_v) { + status = cublasSdot(handle, static_cast(nElmDomain), GetDevicePointer(), 1, other.GetDevicePointer(), 1, + &local_dot); + } else if constexpr (std::is_same_v) { + status = cublasDdot(handle, static_cast(nElmDomain), GetDevicePointer(), 1, other.GetDevicePointer(), 1, + &local_dot); + } else { + ReleaseBlasHandle(handle, owns_handle, "cuBLAS handle destruction failed in CSysVector::GPUDot."); + SU2_MPI::Error("Unsupported ScalarType in CSysVector::GPUDot.", CURRENT_FUNCTION); + return ScalarType(0); + } + + if (status != CUBLAS_STATUS_SUCCESS) { + ReleaseBlasHandle(handle, owns_handle, "cuBLAS handle destruction failed in CSysVector::GPUDot."); + SU2_MPI::Error("cuBLAS dot failed in CSysVector::GPUDot.", CURRENT_FUNCTION); + return ScalarType(0); + } + + ReleaseBlasHandle(handle, owns_handle, "cuBLAS handle destruction failed in CSysVector::GPUDot."); + + ScalarType global_dot = ScalarType(0); + const auto mpi_type = (sizeof(ScalarType) < sizeof(double)) ? MPI_FLOAT : MPI_DOUBLE; + SelectMPIWrapper::W::Allreduce(&local_dot, &global_dot, 1, mpi_type, MPI_SUM, SU2_MPI::GetComm()); - /*--- Implementation area for GPU dot products used by device-resident Krylov solvers. ---*/ - SU2_MPI::Error("CSysVector::GPUDot skeleton reached without an implementation.", CURRENT_FUNCTION); - return ScalarType(0); + return global_dot; } template ScalarType CSysVector::GPUNorm() const { - /*--- Implementation area for GPU norms used by device-resident Krylov solvers. ---*/ - SU2_MPI::Error("CSysVector::GPUNorm skeleton reached without an implementation.", CURRENT_FUNCTION); - return ScalarType(0); + return sqrt(GPUDot(*this)); } template void CSysVector::HtDTransfer(bool trigger) const; From 2b4f9d871610b4100288ac9b8e6966f52029df7c Mon Sep 17 00:00:00 2001 From: LwhJesse <256257451+LwhJesse@users.noreply.github.com> Date: Mon, 1 Jun 2026 22:18:42 +0800 Subject: [PATCH 04/10] Implement CUDA FGMRES solve path --- Common/src/linear_algebra/CSysSolveGPU.cu | 206 ++++++++++++++++++++-- 1 file changed, 190 insertions(+), 16 deletions(-) diff --git a/Common/src/linear_algebra/CSysSolveGPU.cu b/Common/src/linear_algebra/CSysSolveGPU.cu index 3a23e5eda61..cb1ead10d89 100644 --- a/Common/src/linear_algebra/CSysSolveGPU.cu +++ b/Common/src/linear_algebra/CSysSolveGPU.cu @@ -26,8 +26,30 @@ */ #include "../../include/linear_algebra/CSysSolve.hpp" +#include "../../include/linear_algebra/CMatrixVectorProduct.hpp" +#include "../../include/linear_algebra/CPreconditioner.hpp" #include "../../include/linear_algebra/GPUComms.cuh" +#include +#include + +void SU2_GPU_BeginSolverBLASContext(); +void SU2_GPU_EndSolverBLASContext(); + +namespace { + +class CGPUSolverBLASContextGuard { + public: + CGPUSolverBLASContextGuard() { SU2_GPU_BeginSolverBLASContext(); } + + ~CGPUSolverBLASContextGuard() { SU2_GPU_EndSolverBLASContext(); } + + CGPUSolverBLASContextGuard(const CGPUSolverBLASContextGuard&) = delete; + CGPUSolverBLASContextGuard& operator=(const CGPUSolverBLASContextGuard&) = delete; +}; + +} // namespace + template unsigned long CSysSolve::FGMRES_LinSolver_GPU(const VectorType& b, VectorType& x, const ProductType& mat_vec, const PrecondType& precond, @@ -35,22 +57,174 @@ unsigned long CSysSolve::FGMRES_LinSolver_GPU(const VectorType& b, V bool monitoring, const CConfig* config) const { SU2_ZONE_SCOPED - (void)b; - (void)x; - (void)mat_vec; - (void)precond; - (void)tol; - (void)m; - (void)residual; - (void)monitoring; - (void)config; - - /*--- Implementation area for the GPU-resident FGMRES solve path. - * This skeleton intentionally leaves the algorithm unimplemented so the - * surrounding dispatch, build integration, and file structure can be - * reviewed independently of the CUDA solver details. ---*/ - SU2_MPI::Error("FGMRES GPU solver skeleton reached without an implementation.", CURRENT_FUNCTION); - return 0; + const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE); + const bool flexible = !precond.IsIdentity(); + + if (m < 1) { + SU2_MPI::Error("Number of linear solver iterations must be greater than 0.", CURRENT_FUNCTION); + } + + if (m > 1000) { + SU2_MPI::Error("FGMRES subspace is too large (>1000).", CURRENT_FUNCTION); + } + + CGPUSolverBLASContextGuard blas_context; + if (V.size() <= m || (flexible && Z.size() <= m)) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + V.resize(m + 1); + for (auto& w : V) w.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); + if (flexible) { + Z.resize(m + 1); + for (auto& z : Z) z.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); + } + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } + + su2vector g(m + 1), sn(m + 1), cs(m + 1), y(m); + g = ScalarType(0); + sn = ScalarType(0); + cs = ScalarType(0); + y = ScalarType(0); + + su2matrix H(m + 1, m); + H = ScalarType(0); + + /*--- b and x enter from the host side. The solver owns all subsequent + * synchronization for vectors it updates with GPU primitives. ---*/ + b.HtDTransfer(); + if (xIsZero) { + x.GPUSetVal(ScalarType(0)); + } else { + x.HtDTransfer(); + } + + ScalarType norm0 = b.GPUNorm(); + + if (!xIsZero) { + mat_vec(x, V[0]); + V[0].GPUAxpy(ScalarType(-1), b); + } else { + V[0].GPUCopy(b); + V[0].GPUScale(ScalarType(-1)); + } + + ScalarType beta = xIsZero ? norm0 : V[0].GPUNorm(); + + SU2_OMP_MASTER + ResetDeflation(); + END_SU2_OMP_MASTER + + if (tol_type == LinearToleranceType::RELATIVE) norm0 = beta; + + if (beta < tol * norm0 || beta < eps) { + if (masterRank) { + SU2_OMP_MASTER + cout << "CSysSolve::FGMRES_GPU(): system solved by initial guess." << endl; + END_SU2_OMP_MASTER + } + residual = beta; + x.DtHTransfer(); + return 0; + } + + V[0].GPUScale(ScalarType(-1) / beta); + g[0] = beta; + + unsigned long i = 0; + if (monitoring && masterRank) { + SU2_OMP_MASTER { + WriteHeader("FGMRES_GPU", tol, beta); + WriteHistory(i, beta / norm0); + } + END_SU2_OMP_MASTER + } + + for (i = 0; i < m; i++) { + if (beta < tol * norm0) break; + + /*--- mat_vec consumes and produces device-resident vectors on the CUDA path. ---*/ + if (flexible) { + precond(V[i], Z[i]); + mat_vec(Z[i], V[i + 1]); + } else { + mat_vec(V[i], V[i + 1]); + } + + /*--- Classical Gram-Schmidt twice, matching the CPU implementation's + * numerical structure while using GPU dot/axpy/norm primitives. ---*/ + for (unsigned long k = 0; k <= i; k++) { + H(k, i) = V[k].GPUDot(V[i + 1]); + V[i + 1].GPUAxpy(-H(k, i), V[k]); + } + + for (unsigned long k = 0; k <= i; k++) { + const ScalarType dh = V[k].GPUDot(V[i + 1]); + H(k, i) += dh; + V[i + 1].GPUAxpy(-dh, V[k]); + } + + const ScalarType nrm = V[i + 1].GPUNorm(); + if (nrm <= ScalarType(0) || nrm != nrm) { + H(i + 1, i) = ScalarType(0); + if (masterRank) { + SU2_OMP_MASTER + cout << "WARNING: FGMRES GPU orthogonalization failed, linear solver diverged." << endl; + END_SU2_OMP_MASTER + } + break; + } + + H(i + 1, i) = nrm; + V[i + 1].GPUScale(ScalarType(1) / nrm); + + for (unsigned long k = 0; k < i; k++) ApplyGivens(sn[k], cs[k], H(k, i), H(k + 1, i)); + GenerateGivens(H(i, i), H(i + 1, i), sn[i], cs[i]); + ApplyGivens(sn[i], cs[i], g[i], g[i + 1]); + + beta = fabs(g[i + 1]); + + if (monitoring && masterRank && ((i + 1) % monitorFreq == 0)) { + SU2_OMP_MASTER + WriteHistory(i + 1, beta / norm0); + END_SU2_OMP_MASTER + } + } + + SolveReduced(i, H, g, y); + + const auto& basis = flexible ? Z : V; + + for (unsigned long k = 0; k < i; k++) { + x.GPUAxpy(y[k], basis[k]); + } + + x.DtHTransfer(); + + if (monitoring && config->GetComm_Level() == COMM_FULL) { + if (masterRank) { + SU2_OMP_MASTER + WriteFinalResidual("FGMRES_GPU", i, beta / norm0); + END_SU2_OMP_MASTER + } + + if (recomputeRes) { + mat_vec(x, V[0]); + V[0].GPUAxpy(ScalarType(-1), b); + const ScalarType res = V[0].GPUNorm(); + + if (fabs(res - beta) > tol * 10) { + if (masterRank) { + SU2_OMP_MASTER + WriteWarning(beta, res, tol); + END_SU2_OMP_MASTER + } + } + } + } + + residual = beta / norm0; + return i; } /*--- Explicit instantiations for the GPU solver helper. From 9c344ee7934fca81a5f50e95057e5ee54e82eda1 Mon Sep 17 00:00:00 2001 From: LwhJesse <256257451+LwhJesse@users.noreply.github.com> Date: Mon, 1 Jun 2026 22:18:52 +0800 Subject: [PATCH 05/10] Implement CUDA Jacobi preconditioner --- .../linear_algebra/CSysPreconditionerGPU.cu | 69 +++++++++++++++---- 1 file changed, 54 insertions(+), 15 deletions(-) diff --git a/Common/src/linear_algebra/CSysPreconditionerGPU.cu b/Common/src/linear_algebra/CSysPreconditionerGPU.cu index 8bef1c68be9..adee93264eb 100644 --- a/Common/src/linear_algebra/CSysPreconditionerGPU.cu +++ b/Common/src/linear_algebra/CSysPreconditionerGPU.cu @@ -25,34 +25,73 @@ * License along with SU2. If not, see . */ -#include "../../include/linear_algebra/CSysMatrix.hpp" +#include "../../include/linear_algebra/CSysMatrix.inl" #include "../../include/linear_algebra/GPUComms.cuh" +namespace { + +template +__global__ void ApplyJacobiPreconditionerKernel(const ScalarType* invM, const ScalarType* vec, ScalarType* prod, + unsigned long nPointDomain, unsigned long nVar) { + const auto iPoint = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (iPoint >= nPointDomain) return; + + const auto block = &invM[iPoint * nVar * nVar]; + const auto rhs = &vec[iPoint * nVar]; + auto out = &prod[iPoint * nVar]; + + for (auto iVar = 0ul; iVar < nVar; ++iVar) { + ScalarType sum = ScalarType(0); + for (auto jVar = 0ul; jVar < nVar; ++jVar) { + sum += block[iVar * nVar + jVar] * rhs[jVar]; + } + out[iVar] = sum; + } +} + +} // namespace + template void CSysMatrix::BuildJacobiPreconditionerGPU() { - /*--- Implementation area for building the GPU/device Jacobi preconditioner. - * This skeleton intentionally leaves the algorithm unimplemented so the - * surrounding dispatch and file structure can be reviewed independently of - * the CUDA preconditioner details. ---*/ - SU2_MPI::Error("CSysMatrix::BuildJacobiPreconditionerGPU skeleton reached without an implementation.", - CURRENT_FUNCTION); + SU2_ZONE_SCOPED + + if (nVar != nEqn) { + SU2_MPI::Error("CUDA Jacobi preconditioner requires square blocks.", CURRENT_FUNCTION); + } + + if (invM == nullptr) { + SU2_MPI::Error("CUDA Jacobi preconditioner was requested without host inverse block storage.", CURRENT_FUNCTION); + } + + if (d_invM == nullptr) { + d_invM = GPUMemoryAllocation::gpu_alloc(nPointDomain * nVar * nEqn * sizeof(ScalarType)); + } + + for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) { + InverseDiagonalBlock(iPoint, &(invM[iPoint * nVar * nVar])); + } + + gpuErrChk(cudaMemcpy(d_invM, invM, nPointDomain * nVar * nVar * sizeof(ScalarType), cudaMemcpyHostToDevice)); } template void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const { - (void)vec; - (void)prod; (void)geometry; (void)config; - /*--- Implementation area for applying the GPU/device Jacobi preconditioner. - * This skeleton intentionally leaves the algorithm unimplemented so the - * surrounding dispatch and file structure can be reviewed independently of - * the CUDA preconditioner details. ---*/ - SU2_MPI::Error("CSysMatrix::ComputeJacobiPreconditionerGPU skeleton reached without an implementation.", - CURRENT_FUNCTION); + SU2_ZONE_SCOPED + + if (d_invM == nullptr) { + SU2_MPI::Error("CUDA Jacobi preconditioner used before BuildJacobiPreconditionerGPU.", CURRENT_FUNCTION); + } + + constexpr unsigned threadsPerBlock = 128; + const auto blocks = static_cast((nPointDomain + threadsPerBlock - 1) / threadsPerBlock); + ApplyJacobiPreconditionerKernel<<>>(d_invM, vec.GetDevicePointer(), prod.GetDevicePointer(), + nPointDomain, nVar); + gpuErrChk(cudaPeekAtLastError()); } template void CSysMatrix::BuildJacobiPreconditionerGPU(); From a875f56767115bfe866bc99b878d20b4319eba1b Mon Sep 17 00:00:00 2001 From: LwhJesse <256257451+LwhJesse@users.noreply.github.com> Date: Thu, 18 Jun 2026 13:58:16 +0800 Subject: [PATCH 06/10] Share FGMRES control flow with CUDA vector dispatch Move the FGMRES iteration into one shared implementation and select host or CUDA vector-operation backends from the existing solver entry point. Keep the CUDA path GPU-resident for SpMV, Jacobi, dot/norm, and vector updates, with only scalar reductions and final solution synchronization crossing back to the host. --- Common/include/linear_algebra/CSysSolve.hpp | 12 +- .../linear_algebra/CSysSolveFGMRES.inl | 267 ++++++++++++++++++ Common/include/linear_algebra/CSysVector.hpp | 151 ++++++++-- .../linear_algebra/vector_expressions.hpp | 41 +-- Common/src/linear_algebra/CSysSolve.cpp | 223 ++------------- Common/src/linear_algebra/CSysSolveGPU.cu | 200 +++---------- Common/src/linear_algebra/CSysVectorGPU.cu | 80 +----- 7 files changed, 494 insertions(+), 480 deletions(-) create mode 100644 Common/include/linear_algebra/CSysSolveFGMRES.inl diff --git a/Common/include/linear_algebra/CSysSolve.hpp b/Common/include/linear_algebra/CSysSolve.hpp index 67033992372..af76e93c07d 100644 --- a/Common/include/linear_algebra/CSysSolve.hpp +++ b/Common/include/linear_algebra/CSysSolve.hpp @@ -304,11 +304,15 @@ class CSysSolve { ScalarType& residual, bool monitoring, const CConfig* config, FgcrodrMode mode, unsigned long custom_m) const; + template + unsigned long FGMRES_LinSolverImpl(const VectorType& b, VectorType& x, const ProductType& mat_vec, + const PrecondType& precond, ScalarType tol, unsigned long m, ScalarType& residual, + bool monitoring, const CConfig* config, const FGMRESOps& ops) const; + /*! - * \brief CUDA/GPU implementation of the Flexible GMRES solver. - * \note This helper exists so the public solver interface can remain unchanged while the - * implementation is dispatched internally based on the runtime CUDA setting. - * The actual implementation lives in CSysSolveGPU.cu. + * \brief CUDA/GPU backend wrapper for the Flexible GMRES solver. + * \note The algorithmic implementation is shared with the host solver. This helper + * only provides GPU vector-operation dispatch and synchronization. */ unsigned long FGMRES_LinSolver_GPU(const VectorType& b, VectorType& x, const ProductType& mat_vec, const PrecondType& precond, ScalarType tol, unsigned long m, ScalarType& residual, diff --git a/Common/include/linear_algebra/CSysSolveFGMRES.inl b/Common/include/linear_algebra/CSysSolveFGMRES.inl new file mode 100644 index 00000000000..cef424a729b --- /dev/null +++ b/Common/include/linear_algebra/CSysSolveFGMRES.inl @@ -0,0 +1,267 @@ +/*! + * \file CSysSolveFGMRES.inl + * \brief Shared implementation of the Flexible GMRES linear solver. + * \author Jesse Li + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +namespace fgmres_impl { + +template +bool Orthogonalize(unsigned long i, su2matrix& H, std::vector>& V, + const FGMRESOps& ops) { + for (unsigned long k = 0; k <= i; k++) { + H(k, i) = ops.Dot(V[k], V[i + 1]); + ops.AddScaled(V[i + 1], -H(k, i), V[k]); + } + + for (unsigned long k = 0; k <= i; k++) { + const ScalarType dh = ops.Dot(V[k], V[i + 1]); + H(k, i) += dh; + ops.AddScaled(V[i + 1], -dh, V[k]); + } + + const ScalarType nrm = ops.Norm(V[i + 1]); + if (nrm <= ScalarType(0) || nrm != nrm) { + H(i + 1, i) = ScalarType(0); + return false; + } + + H(i + 1, i) = nrm; + ops.Divide(V[i + 1], nrm); + return true; +} + +} // namespace fgmres_impl + +template +template +unsigned long CSysSolve::FGMRES_LinSolverImpl(const VectorType& b, VectorType& x, + const ProductType& mat_vec, const PrecondType& precond, + ScalarType tol, unsigned long m, ScalarType& residual, + bool monitoring, const CConfig* config, + const FGMRESOps& ops) const { + SU2_ZONE_SCOPED + + const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE); + const bool flexible = !precond.IsIdentity(); + const bool nestedParallel = ops.NestedParallel(); + + /*--- Check the subspace size. ---*/ + + if (m < 1) { + SU2_MPI::Error("Number of linear solver iterations must be greater than 0.", CURRENT_FUNCTION); + } + + if (m > 1000) { + SU2_MPI::Error("FGMRES subspace is too large (>1000).", CURRENT_FUNCTION); + } + + /*--- Allocate if not allocated yet. ---*/ + + if (V.size() <= m || (flexible && Z.size() <= m)) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + V.resize(m + 1); + for (auto& w : V) w.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); + if (flexible) { + Z.resize(m + 1); + for (auto& z : Z) z.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); + } + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } + + /*--- Define various arrays. In parallel, each thread of each rank has and works + on its own thread, since calculations on these arrays are based on dot products + (reduced across all threads and ranks) all threads do the same computations. ---*/ + + su2vector g(m + 1), sn(m + 1), cs(m + 1), y(m); + g = ScalarType(0); + sn = ScalarType(0); + cs = ScalarType(0); + y = ScalarType(0); + + su2matrix H(m + 1, m); + H = ScalarType(0); + + ops.PrepareInput(b, x, xIsZero); + + /*--- Calculate the norm of the rhs vector. ---*/ + + ScalarType norm0 = ops.Norm(b); + + /*--- Calculate the initial residual (actually the negative residual) and compute its norm. ---*/ + + if (!xIsZero) { + mat_vec(x, V[0]); + ops.Subtract(V[0], b); + } else { + ops.AssignNegative(V[0], b); + } + + ScalarType beta = xIsZero ? norm0 : ops.Norm(V[0]); + + /*--- FGMRES shares V and Z with FGCRODR, reset deflation because we're breaking + * relations between subspaces that FGCRODR relies on for recycling. ---*/ + + SU2_OMP_MASTER + ResetDeflation(); + END_SU2_OMP_MASTER + + /*--- Set the norm to the initial initial residual value. ---*/ + + if (tol_type == LinearToleranceType::RELATIVE) norm0 = beta; + + if (beta < tol * norm0 || beta < eps) { + /*--- System is already solved. ---*/ + + if (masterRank) { + SU2_OMP_MASTER + cout << "CSysSolve::FGMRES(): system solved by initial guess." << endl; + END_SU2_OMP_MASTER + } + residual = beta; + ops.FinalizeOutput(x); + return 0; + } + + /*--- Normalize residual to get w_{0} (the negative sign is because w[0] + holds the negative residual, as mentioned above). ---*/ + + ops.Divide(V[0], -beta); + + /*--- Initialize the RHS of the reduced system. ---*/ + + g[0] = beta; + + /*--- Output header information including initial residual. ---*/ + + unsigned long i = 0; + if (monitoring && masterRank) { + SU2_OMP_MASTER { + WriteHeader("FGMRES", tol, beta); + WriteHistory(i, beta / norm0); + } + END_SU2_OMP_MASTER + } + + /*--- Loop over all search directions. ---*/ + + for (i = 0; i < m; i++) { + /*--- Check if solution has converged. ---*/ + + if (beta < tol * norm0) break; + + if (flexible) { + /*--- Precondition the CSysVector v[i] and store result in z[i]. ---*/ + + precond(V[i], Z[i]); + + /*--- Add to Krylov subspace. ---*/ + + mat_vec(Z[i], V[i + 1]); + } else { + mat_vec(V[i], V[i + 1]); + } + + /*--- Modified Gram-Schmidt orthogonalization. ---*/ + + bool orthog_ok; + if constexpr (FGMRESOps::UseSolverModGramSchmidt) { + if (nestedParallel) { + /*--- "omp parallel if" does not work well here. ---*/ + SU2_OMP_PARALLEL + orthog_ok = ModGramSchmidt(true, i, H, V); + END_SU2_OMP_PARALLEL + } else { + orthog_ok = ModGramSchmidt(false, i, H, V); + } + } else { + orthog_ok = fgmres_impl::Orthogonalize(i, H, V, ops); + } + + if (!orthog_ok) { + if (masterRank) { + SU2_OMP_MASTER + cout << "WARNING: FGMRES orthogonalization failed, linear solver diverged." << endl; + END_SU2_OMP_MASTER + } + break; + } + + /*--- Apply old Givens rotations to new column of the Hessenberg matrix then generate the + new Givens rotation matrix and apply it to the last two elements of H[:][i] and g. ---*/ + + for (unsigned long k = 0; k < i; k++) ApplyGivens(sn[k], cs[k], H(k, i), H(k + 1, i)); + GenerateGivens(H(i, i), H(i + 1, i), sn[i], cs[i]); + ApplyGivens(sn[i], cs[i], g[i], g[i + 1]); + + /*--- Set L2 norm of residual and check if solution has converged. ---*/ + + beta = fabs(g[i + 1]); + + /*--- Output the relative residual if necessary. ---*/ + + if (monitoring && masterRank && ((i + 1) % monitorFreq == 0)) { + SU2_OMP_MASTER + WriteHistory(i + 1, beta / norm0); + END_SU2_OMP_MASTER + } + } + + /*--- Solve the least-squares system and update solution. ---*/ + + SolveReduced(i, H, g, y); + + const auto& basis = flexible ? Z : V; + + ops.UpdateSolution(i, basis, y, x, nestedParallel); + ops.FinalizeOutput(x); + + /*--- Recalculate final (neg.) residual (this should be optional). ---*/ + + if (monitoring && config->GetComm_Level() == COMM_FULL) { + if (masterRank) { + SU2_OMP_MASTER + WriteFinalResidual("FGMRES", i, beta / norm0); + END_SU2_OMP_MASTER + } + + if (recomputeRes) { + mat_vec(x, V[0]); + ops.Subtract(V[0], b); + ScalarType res = ops.Norm(V[0]); + + if (fabs(res - beta) > tol * 10) { + if (masterRank) { + SU2_OMP_MASTER + WriteWarning(beta, res, tol); + END_SU2_OMP_MASTER + } + } + } + } + + residual = beta / norm0; + return i; +} diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index decffe26d29..d67185156b9 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -37,6 +37,10 @@ #include "vector_expressions.hpp" #include "../../include/CConfig.hpp" +#ifdef __CUDACC__ +#include "GPUComms.cuh" +#endif + /*! * \brief OpenMP worksharing construct used in CSysVector for loops. * \note The loop will only run in parallel if methods are called from a @@ -59,6 +63,116 @@ #define END_CSYSVEC_PARFOR #endif +namespace VecExpr { + +enum class DeviceAssignOp { Assign, Add, Subtract, Multiply, Divide }; + +template +class CDeviceVectorView : public CVecExpr, Scalar> { + private: + Scalar* data = nullptr; + unsigned long size = 0; + + public: + static constexpr bool StoreAsRef = false; + + CDeviceVectorView(Scalar* data_, unsigned long size_) : data(data_), size(size_) {} + +#ifdef __CUDACC__ + __device__ const Scalar& operator[](size_t i) const { return data[i]; } + + template + CDeviceVectorView& operator=(const CVecExpr& expr); + + CDeviceVectorView& operator=(Scalar val); + +#define MAKE_DEVICE_COMPOUND(OP) \ + template \ + CDeviceVectorView& operator OP(const CVecExpr& expr); + MAKE_DEVICE_COMPOUND(+=) + MAKE_DEVICE_COMPOUND(-=) + MAKE_DEVICE_COMPOUND(*=) + MAKE_DEVICE_COMPOUND(/=) +#undef MAKE_DEVICE_COMPOUND + +#define MAKE_DEVICE_SCALAR_COMPOUND(OP) CDeviceVectorView& operator OP(Scalar val); + MAKE_DEVICE_SCALAR_COMPOUND(+=) + MAKE_DEVICE_SCALAR_COMPOUND(-=) + MAKE_DEVICE_SCALAR_COMPOUND(*=) + MAKE_DEVICE_SCALAR_COMPOUND(/=) +#undef MAKE_DEVICE_SCALAR_COMPOUND +#endif +}; + +#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; + if (i >= size) return; + + if constexpr (Op == DeviceAssignOp::Assign) { + data[i] = expr[i]; + } else if constexpr (Op == DeviceAssignOp::Add) { + data[i] += expr[i]; + } else if constexpr (Op == DeviceAssignOp::Subtract) { + data[i] -= expr[i]; + } else if constexpr (Op == DeviceAssignOp::Multiply) { + data[i] *= expr[i]; + } else { + data[i] /= expr[i]; + } +} + +template +inline void AssignDeviceExpression(Scalar* data, unsigned long 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); + DeviceAssignKernel<<>>(data, size, expr.derived()); + gpuErrChk(cudaPeekAtLastError()); +} + +template +template +CDeviceVectorView& CDeviceVectorView::operator=(const CVecExpr& expr) { + AssignDeviceExpression(data, size, expr); + return *this; +} + +template +CDeviceVectorView& CDeviceVectorView::operator=(Scalar val) { + AssignDeviceExpression(data, size, Bcast(val)); + return *this; +} + +#define MAKE_DEVICE_COMPOUND(OP, ASSIGN_OP) \ + template \ + template \ + CDeviceVectorView& CDeviceVectorView::operator OP(const CVecExpr& expr) { \ + AssignDeviceExpression(data, size, expr); \ + return *this; \ + } +MAKE_DEVICE_COMPOUND(+=, DeviceAssignOp::Add) +MAKE_DEVICE_COMPOUND(-=, DeviceAssignOp::Subtract) +MAKE_DEVICE_COMPOUND(*=, DeviceAssignOp::Multiply) +MAKE_DEVICE_COMPOUND(/=, DeviceAssignOp::Divide) +#undef MAKE_DEVICE_COMPOUND + +#define MAKE_DEVICE_SCALAR_COMPOUND(OP, ASSIGN_OP) \ + template \ + CDeviceVectorView& CDeviceVectorView::operator OP(Scalar val) { \ + AssignDeviceExpression(data, size, Bcast(val)); \ + return *this; \ + } +MAKE_DEVICE_SCALAR_COMPOUND(+=, DeviceAssignOp::Add) +MAKE_DEVICE_SCALAR_COMPOUND(-=, DeviceAssignOp::Subtract) +MAKE_DEVICE_SCALAR_COMPOUND(*=, DeviceAssignOp::Multiply) +MAKE_DEVICE_SCALAR_COMPOUND(/=, DeviceAssignOp::Divide) +#undef MAKE_DEVICE_SCALAR_COMPOUND +#endif + +} // namespace VecExpr + /*! * \class CSysVector * \ingroup SpLinSys @@ -234,36 +348,6 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ void DtHTransfer(bool trigger = true) const; - /*! - * \brief Sets all the elements of the GPU vector to a certain value - * \param[in] trigger - boolean value that decides whether to conduct the transfer or not. True by default. - */ - void GPUSetVal(ScalarType val, bool trigger = true) const; - - /*! - * \brief Copy another vector into this vector on the device. - * \note This is an explicit GPU helper for Krylov solver implementations that keep the - * iteration state resident on the device. The implementation belongs in - * CSysVectorGPU.cu. - * \param[in] src - Source vector. - */ - void GPUCopy(const CSysVector& src) const; - - /*! - * \brief Scale this vector on the device. - * \note Explicit GPU helper for solver-side vector operations. - * \param[in] alpha - Scalar multiplier. - */ - void GPUScale(ScalarType alpha) const; - - /*! - * \brief Perform the AXPY operation on the device: this := this + alpha * x. - * \note Explicit GPU helper for solver-side vector operations. - * \param[in] alpha - Scalar multiplier. - * \param[in] x - Input vector. - */ - void GPUAxpy(ScalarType alpha, const CSysVector& x) const; - /*! * \brief Dot product between this vector and another vector on the device. * \note Explicit GPU helper for solver-side reductions. @@ -284,6 +368,13 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ inline ScalarType* GetDevicePointer() const { return d_vec_val; } + /*! + * \brief Return an expression-template view of the device values. + */ + inline VecExpr::CDeviceVectorView device() const { + return VecExpr::CDeviceVectorView(d_vec_val, nElm); + } + /*! * \brief return the number of local elements in the CSysVector */ diff --git a/Common/include/linear_algebra/vector_expressions.hpp b/Common/include/linear_algebra/vector_expressions.hpp index a0d0ce28901..b2c6a50f5a7 100644 --- a/Common/include/linear_algebra/vector_expressions.hpp +++ b/Common/include/linear_algebra/vector_expressions.hpp @@ -39,6 +39,12 @@ namespace VecExpr { /// \addtogroup VecExpr /// @{ +#ifdef __CUDACC__ +#define SU2_CUDA_HOST_DEVICE __host__ __device__ +#else +#define SU2_CUDA_HOST_DEVICE +#endif + /*! * \brief Base vector expression class. * \ingroup BLAS @@ -59,7 +65,7 @@ class CVecExpr { /*! * \brief Cast the expression to Derived, usually to allow evaluation via operator[]. */ - FORCEINLINE const Derived& derived() const { return static_cast(*this); } + SU2_CUDA_HOST_DEVICE FORCEINLINE const Derived& derived() const { return static_cast(*this); } // Allowed from C++14, allows nested expression propagation without // manually calling derived() on the expression being evaluated. @@ -76,8 +82,8 @@ class Bcast : public CVecExpr, Scalar> { public: static constexpr bool StoreAsRef = false; - FORCEINLINE Bcast(const Scalar& x_) : x(x_) {} - FORCEINLINE const Scalar& operator[](size_t) const { return x; } + SU2_CUDA_HOST_DEVICE FORCEINLINE Bcast(const Scalar& x_) : x(x_) {} + SU2_CUDA_HOST_DEVICE FORCEINLINE const Scalar& operator[](size_t) const { return x; } }; /*! @@ -120,19 +126,19 @@ namespace math = ::std; /*--- Macro to create expression classes (EXPR) and overloads (FUN) for unary * functions, based on their coefficient-wise implementation (IMPL). ---*/ -#define MAKE_UNARY_FUN(FUN, EXPR, IMPL) \ - /*!--- Expression class. ---*/ \ - template \ - class EXPR : public CVecExpr, Scalar> { \ - store_t u; \ - \ - public: \ - static constexpr bool StoreAsRef = false; \ - FORCEINLINE EXPR(const U& u_) : u(u_) {} \ - FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i])) \ - }; \ - /*!--- Function overload, returns an expression object. ---*/ \ - template \ +#define MAKE_UNARY_FUN(FUN, EXPR, IMPL) \ + /*!--- Expression class. ---*/ \ + template \ + class EXPR : public CVecExpr, Scalar> { \ + store_t u; \ + \ + public: \ + static constexpr bool StoreAsRef = false; \ + FORCEINLINE EXPR(const U& u_) : u(u_) {} \ + SU2_CUDA_HOST_DEVICE FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i])) \ + }; \ + /*!--- Function overload, returns an expression object. ---*/ \ + template \ FORCEINLINE auto FUN(const CVecExpr& u) RETURNS(EXPR(u.derived())) #define sign_impl(x) Scalar(1 - 2 * (x < 0)) @@ -158,7 +164,7 @@ MAKE_UNARY_FUN(sign, sign_, sign_impl) public: \ static constexpr bool StoreAsRef = false; \ FORCEINLINE EXPR(const U& u_, const V& v_) : u(u_), v(v_) {} \ - FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i], v[i])) \ + SU2_CUDA_HOST_DEVICE FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i], v[i])) \ }; \ /*!--- Vector with vector function overload. ---*/ \ template \ @@ -241,4 +247,5 @@ MAKE_BINARY_FUN(operator>, gt_, gt_impl) #undef MAKE_BINARY_FUN /// @} +#undef SU2_CUDA_HOST_DEVICE } // namespace VecExpr diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index 3eba623f051..1045ae68bbc 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -131,8 +131,35 @@ void LinearCombination(bool parallel, Ts&&... args) { } } +template +class CHostFGMRESOps { + public: + static constexpr bool UseSolverModGramSchmidt = true; + + bool NestedParallel() const { return !omp_in_parallel() && omp_get_max_threads() > 1; } + + void PrepareInput(const CSysVector&, CSysVector&, bool) const {} + + void FinalizeOutput(CSysVector&) const {} + + ScalarType Norm(const CSysVector& v) const { return v.norm(); } + + void AssignNegative(CSysVector& dst, const CSysVector& src) const { dst = -src; } + + void Subtract(CSysVector& dst, const CSysVector& src) const { dst -= src; } + + void Divide(CSysVector& dst, ScalarType val) const { dst /= val; } + + void UpdateSolution(unsigned long n, const std::vector>& basis, const su2vector& y, + CSysVector& x, bool nestedParallel) const { + LinearCombination(nestedParallel, n, basis, y, x, true); + } +}; + } // namespace +#include "../../include/linear_algebra/CSysSolveFGMRES.inl" + template CSysSolve::CSysSolve(LINEAR_SOLVER_MODE linear_solver_mode) : eps(linSolEpsilon()), @@ -432,200 +459,8 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector 1; - - /*--- Check the subspace size ---*/ - - if (m < 1) { - SU2_MPI::Error("Number of linear solver iterations must be greater than 0.", CURRENT_FUNCTION); - } - - if (m > 1000) { - SU2_MPI::Error("FGMRES subspace is too large (>1000).", CURRENT_FUNCTION); - } - - /*--- Allocate if not allocated yet ---*/ - - if (V.size() <= m || (flexible && Z.size() <= m)) { - BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - V.resize(m + 1); - for (auto& w : V) w.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); - if (flexible) { - Z.resize(m + 1); - for (auto& z : Z) z.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); - } - } - END_SU2_OMP_SAFE_GLOBAL_ACCESS - } - - /*--- Define various arrays. In parallel, each thread of each rank has and works - on its own thread, since calculations on these arrays are based on dot products - (reduced across all threads and ranks) all threads do the same computations. ---*/ - - su2vector g(m + 1), sn(m + 1), cs(m + 1), y(m); - g = ScalarType(0); - sn = ScalarType(0); - cs = ScalarType(0); - y = ScalarType(0); - su2matrix H(m + 1, m); - H = ScalarType(0); - - /*--- Calculate the norm of the rhs vector. ---*/ - - ScalarType norm0 = b.norm(); - - /*--- Calculate the initial residual (actually the negative residual) and compute its norm. ---*/ - - if (!xIsZero) { - mat_vec(x, V[0]); - V[0] -= b; - } else { - V[0] = -b; - } - - ScalarType beta = xIsZero ? norm0 : V[0].norm(); - - /*--- FGMRES shares V and Z with FGCRODR, reset deflation because we're breaking - * relations between subspaces that FGCRODR relies on for recycling. ---*/ - - SU2_OMP_MASTER - ResetDeflation(); - END_SU2_OMP_MASTER - - /*--- Set the norm to the initial initial residual value ---*/ - - if (tol_type == LinearToleranceType::RELATIVE) norm0 = beta; - - if (beta < tol * norm0 || beta < eps) { - /*--- System is already solved ---*/ - - if (masterRank) { - SU2_OMP_MASTER - cout << "CSysSolve::FGMRES(): system solved by initial guess." << endl; - END_SU2_OMP_MASTER - } - residual = beta; - return 0; - } - - /*--- Normalize residual to get w_{0} (the negative sign is because w[0] - holds the negative residual, as mentioned above). ---*/ - - V[0] /= -beta; - - /*--- Initialize the RHS of the reduced system ---*/ - - g[0] = beta; - - /*--- Output header information including initial residual ---*/ - - unsigned long i = 0; - if (monitoring && masterRank) { - SU2_OMP_MASTER { - WriteHeader("FGMRES", tol, beta); - WriteHistory(i, beta / norm0); - } - END_SU2_OMP_MASTER - } - - /*--- Loop over all search directions ---*/ - - for (i = 0; i < m; i++) { - /*--- Check if solution has converged ---*/ - - if (beta < tol * norm0) break; - - if (flexible) { - /*--- Precondition the CSysVector v[i] and store result in z[i] ---*/ - - precond(V[i], Z[i]); - - /*--- Add to Krylov subspace ---*/ - - mat_vec(Z[i], V[i + 1]); - } else { - mat_vec(V[i], V[i + 1]); - } - - /*--- Modified Gram-Schmidt orthogonalization ---*/ - - bool orthog_ok; - if (nestedParallel) { - /*--- "omp parallel if" does not work well here ---*/ - SU2_OMP_PARALLEL - orthog_ok = ModGramSchmidt(true, i, H, V); - END_SU2_OMP_PARALLEL - } else { - orthog_ok = ModGramSchmidt(false, i, H, V); - } - - if (!orthog_ok) { - if (masterRank) { - SU2_OMP_MASTER - cout << "WARNING: FGMRES orthogonalization failed, linear solver diverged." << endl; - END_SU2_OMP_MASTER - } - break; - } - - /*--- Apply old Givens rotations to new column of the Hessenberg matrix then generate the - new Givens rotation matrix and apply it to the last two elements of H[:][i] and g ---*/ - - for (unsigned long k = 0; k < i; k++) ApplyGivens(sn[k], cs[k], H(k, i), H(k + 1, i)); - GenerateGivens(H(i, i), H(i + 1, i), sn[i], cs[i]); - ApplyGivens(sn[i], cs[i], g[i], g[i + 1]); - - /*--- Set L2 norm of residual and check if solution has converged ---*/ - - beta = fabs(g[i + 1]); - - /*--- Output the relative residual if necessary ---*/ - - if (monitoring && masterRank && ((i + 1) % monitorFreq == 0)) { - SU2_OMP_MASTER - WriteHistory(i + 1, beta / norm0); - END_SU2_OMP_MASTER - } - } - - /*--- Solve the least-squares system and update solution ---*/ - - SolveReduced(i, H, g, y); - - const auto& basis = flexible ? Z : V; - - LinearCombination(nestedParallel, i, basis, y, x, true); - - /*--- Recalculate final (neg.) residual (this should be optional) ---*/ - - if (monitoring && config->GetComm_Level() == COMM_FULL) { - if (masterRank) { - SU2_OMP_MASTER - WriteFinalResidual("FGMRES", i, beta / norm0); - END_SU2_OMP_MASTER - } - - if (recomputeRes) { - mat_vec(x, V[0]); - V[0] -= b; - ScalarType res = V[0].norm(); - - if (fabs(res - beta) > tol * 10) { - if (masterRank) { - SU2_OMP_MASTER - WriteWarning(beta, res, tol); - END_SU2_OMP_MASTER - } - } - } - } - - residual = beta / norm0; - return i; + return FGMRES_LinSolverImpl(b, x, mat_vec, precond, tol, m, residual, monitoring, config, + CHostFGMRESOps()); } template diff --git a/Common/src/linear_algebra/CSysSolveGPU.cu b/Common/src/linear_algebra/CSysSolveGPU.cu index cb1ead10d89..d618cc223b5 100644 --- a/Common/src/linear_algebra/CSysSolveGPU.cu +++ b/Common/src/linear_algebra/CSysSolveGPU.cu @@ -1,6 +1,6 @@ /*! * \file CSysSolveGPU.cu - * \brief CUDA/GPU skeleton implementations for Krylov solvers. + * \brief CUDA/GPU backend dispatch for Krylov solvers. * \author Jesse Li * \version 8.5.0 "Harrier" * @@ -30,9 +30,6 @@ #include "../../include/linear_algebra/CPreconditioner.hpp" #include "../../include/linear_algebra/GPUComms.cuh" -#include -#include - void SU2_GPU_BeginSolverBLASContext(); void SU2_GPU_EndSolverBLASContext(); @@ -48,186 +45,67 @@ class CGPUSolverBLASContextGuard { CGPUSolverBLASContextGuard& operator=(const CGPUSolverBLASContextGuard&) = delete; }; -} // namespace - template -unsigned long CSysSolve::FGMRES_LinSolver_GPU(const VectorType& b, VectorType& x, - const ProductType& mat_vec, const PrecondType& precond, - ScalarType tol, unsigned long m, ScalarType& residual, - bool monitoring, const CConfig* config) const { - SU2_ZONE_SCOPED - - const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE); - const bool flexible = !precond.IsIdentity(); - - if (m < 1) { - SU2_MPI::Error("Number of linear solver iterations must be greater than 0.", CURRENT_FUNCTION); - } +class CDeviceFGMRESOps { + public: + static constexpr bool UseSolverModGramSchmidt = false; - if (m > 1000) { - SU2_MPI::Error("FGMRES subspace is too large (>1000).", CURRENT_FUNCTION); - } + bool NestedParallel() const { return false; } - CGPUSolverBLASContextGuard blas_context; - if (V.size() <= m || (flexible && Z.size() <= m)) { - BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - V.resize(m + 1); - for (auto& w : V) w.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); - if (flexible) { - Z.resize(m + 1); - for (auto& z : Z) z.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); - } + void PrepareInput(const CSysVector& b, CSysVector& x, bool x_is_zero) const { + b.HtDTransfer(); + if (x_is_zero) { + x.device() = ScalarType(0); + } else { + x.HtDTransfer(); } - END_SU2_OMP_SAFE_GLOBAL_ACCESS } - su2vector g(m + 1), sn(m + 1), cs(m + 1), y(m); - g = ScalarType(0); - sn = ScalarType(0); - cs = ScalarType(0); - y = ScalarType(0); - - su2matrix H(m + 1, m); - H = ScalarType(0); - - /*--- b and x enter from the host side. The solver owns all subsequent - * synchronization for vectors it updates with GPU primitives. ---*/ - b.HtDTransfer(); - if (xIsZero) { - x.GPUSetVal(ScalarType(0)); - } else { - x.HtDTransfer(); - } + void FinalizeOutput(CSysVector& x) const { x.DtHTransfer(); } - ScalarType norm0 = b.GPUNorm(); + ScalarType Norm(const CSysVector& v) const { return v.GPUNorm(); } - if (!xIsZero) { - mat_vec(x, V[0]); - V[0].GPUAxpy(ScalarType(-1), b); - } else { - V[0].GPUCopy(b); - V[0].GPUScale(ScalarType(-1)); + void AssignNegative(CSysVector& dst, const CSysVector& src) const { + dst.device() = -src.device(); } - ScalarType beta = xIsZero ? norm0 : V[0].GPUNorm(); - - SU2_OMP_MASTER - ResetDeflation(); - END_SU2_OMP_MASTER - - if (tol_type == LinearToleranceType::RELATIVE) norm0 = beta; - - if (beta < tol * norm0 || beta < eps) { - if (masterRank) { - SU2_OMP_MASTER - cout << "CSysSolve::FGMRES_GPU(): system solved by initial guess." << endl; - END_SU2_OMP_MASTER - } - residual = beta; - x.DtHTransfer(); - return 0; + void Subtract(CSysVector& dst, const CSysVector& src) const { + dst.device() -= src.device(); } - V[0].GPUScale(ScalarType(-1) / beta); - g[0] = beta; + void Divide(CSysVector& dst, ScalarType val) const { dst.device() /= val; } - unsigned long i = 0; - if (monitoring && masterRank) { - SU2_OMP_MASTER { - WriteHeader("FGMRES_GPU", tol, beta); - WriteHistory(i, beta / norm0); - } - END_SU2_OMP_MASTER + ScalarType Dot(const CSysVector& lhs, const CSysVector& rhs) const { + return lhs.GPUDot(rhs); } - for (i = 0; i < m; i++) { - if (beta < tol * norm0) break; - - /*--- mat_vec consumes and produces device-resident vectors on the CUDA path. ---*/ - if (flexible) { - precond(V[i], Z[i]); - mat_vec(Z[i], V[i + 1]); - } else { - mat_vec(V[i], V[i + 1]); - } - - /*--- Classical Gram-Schmidt twice, matching the CPU implementation's - * numerical structure while using GPU dot/axpy/norm primitives. ---*/ - for (unsigned long k = 0; k <= i; k++) { - H(k, i) = V[k].GPUDot(V[i + 1]); - V[i + 1].GPUAxpy(-H(k, i), V[k]); - } - - for (unsigned long k = 0; k <= i; k++) { - const ScalarType dh = V[k].GPUDot(V[i + 1]); - H(k, i) += dh; - V[i + 1].GPUAxpy(-dh, V[k]); - } - - const ScalarType nrm = V[i + 1].GPUNorm(); - if (nrm <= ScalarType(0) || nrm != nrm) { - H(i + 1, i) = ScalarType(0); - if (masterRank) { - SU2_OMP_MASTER - cout << "WARNING: FGMRES GPU orthogonalization failed, linear solver diverged." << endl; - END_SU2_OMP_MASTER - } - break; - } - - H(i + 1, i) = nrm; - V[i + 1].GPUScale(ScalarType(1) / nrm); - - for (unsigned long k = 0; k < i; k++) ApplyGivens(sn[k], cs[k], H(k, i), H(k + 1, i)); - GenerateGivens(H(i, i), H(i + 1, i), sn[i], cs[i]); - ApplyGivens(sn[i], cs[i], g[i], g[i + 1]); - - beta = fabs(g[i + 1]); - - if (monitoring && masterRank && ((i + 1) % monitorFreq == 0)) { - SU2_OMP_MASTER - WriteHistory(i + 1, beta / norm0); - END_SU2_OMP_MASTER - } + void AddScaled(CSysVector& dst, ScalarType alpha, const CSysVector& src) const { + dst.device() += alpha * src.device(); } - SolveReduced(i, H, g, y); - - const auto& basis = flexible ? Z : V; - - for (unsigned long k = 0; k < i; k++) { - x.GPUAxpy(y[k], basis[k]); + void UpdateSolution(unsigned long n, const std::vector>& basis, + const su2vector& y, CSysVector& x, bool) const { + for (unsigned long k = 0; k < n; k++) { + AddScaled(x, y[k], basis[k]); + } } +}; - x.DtHTransfer(); - - if (monitoring && config->GetComm_Level() == COMM_FULL) { - if (masterRank) { - SU2_OMP_MASTER - WriteFinalResidual("FGMRES_GPU", i, beta / norm0); - END_SU2_OMP_MASTER - } +} // namespace - if (recomputeRes) { - mat_vec(x, V[0]); - V[0].GPUAxpy(ScalarType(-1), b); - const ScalarType res = V[0].GPUNorm(); - - if (fabs(res - beta) > tol * 10) { - if (masterRank) { - SU2_OMP_MASTER - WriteWarning(beta, res, tol); - END_SU2_OMP_MASTER - } - } - } - } +#include "../../include/linear_algebra/CSysSolveFGMRES.inl" - residual = beta / norm0; - return i; +template +unsigned long CSysSolve::FGMRES_LinSolver_GPU(const VectorType& b, VectorType& x, + const ProductType& mat_vec, const PrecondType& precond, + ScalarType tol, unsigned long m, ScalarType& residual, + bool monitoring, const CConfig* config) const { + CGPUSolverBLASContextGuard blas_context; + return FGMRES_LinSolverImpl(b, x, mat_vec, precond, tol, m, residual, monitoring, config, + CDeviceFGMRESOps()); } -/*--- Explicit instantiations for the GPU solver helper. +/*--- Explicit instantiations for the GPU backend wrapper. * Keep these aligned with the scalar types instantiated for CSysSolve. ---*/ template unsigned long CSysSolve::FGMRES_LinSolver_GPU( const CSysVector& b, CSysVector& x, diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index d69e9fcf234..dca72c0d76b 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -32,31 +32,6 @@ namespace { -constexpr unsigned GPU_VEC_OP_BLOCK_SIZE = 256; - -template -__global__ void GPUCopyKernel(const ScalarType* src, ScalarType* dst, unsigned long nElm) { - const unsigned long idx = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; - if (idx < nElm) dst[idx] = src[idx]; -} - -template -__global__ void GPUScaleKernel(ScalarType alpha, ScalarType* vec, unsigned long nElm) { - const unsigned long idx = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; - if (idx < nElm) vec[idx] *= alpha; -} - -template -__global__ void GPUAxpyKernel(ScalarType alpha, const ScalarType* x, ScalarType* y, unsigned long nElm) { - const unsigned long idx = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; - if (idx < nElm) y[idx] += alpha * x[idx]; -} - -inline dim3 MakeVectorGrid(unsigned long nElm) { - const auto grid = (nElm + GPU_VEC_OP_BLOCK_SIZE - 1) / GPU_VEC_OP_BLOCK_SIZE; - return dim3(static_cast(grid), 1, 1); -} - thread_local cublasHandle_t active_solver_blas_handle = nullptr; thread_local unsigned active_solver_blas_depth = 0; @@ -108,47 +83,16 @@ void SU2_GPU_EndSolverBLASContext() { } } -template -void CSysVector::HtDTransfer(bool trigger) const -{ - if(trigger) gpuErrChk(cudaMemcpy((void*)(d_vec_val), (void*)&vec_val[0], (sizeof(ScalarType)*nElm), cudaMemcpyHostToDevice)); -} - -template -void CSysVector::DtHTransfer(bool trigger) const -{ - if(trigger) gpuErrChk(cudaMemcpy((void*)(&vec_val[0]), (void*)d_vec_val, (sizeof(ScalarType)*nElm), cudaMemcpyDeviceToHost)); -} - -template -void CSysVector::GPUSetVal(ScalarType val, bool trigger) const -{ - if(trigger) gpuErrChk(cudaMemset((void*)(d_vec_val), val, (sizeof(ScalarType)*nElm))); -} - template -void CSysVector::GPUCopy(const CSysVector& src) const { - if (nElm == 0) return; - - GPUCopyKernel<<>>(src.GetDevicePointer(), GetDevicePointer(), nElm); - gpuErrChk(cudaPeekAtLastError()); -} - -template -void CSysVector::GPUScale(ScalarType alpha) const { - if (nElm == 0) return; - - GPUScaleKernel<<>>(alpha, GetDevicePointer(), nElm); - gpuErrChk(cudaPeekAtLastError()); +void CSysVector::HtDTransfer(bool trigger) const { + if (trigger) + gpuErrChk(cudaMemcpy((void*)(d_vec_val), (void*)&vec_val[0], (sizeof(ScalarType) * nElm), cudaMemcpyHostToDevice)); } template -void CSysVector::GPUAxpy(ScalarType alpha, const CSysVector& x) const { - if (nElm == 0) return; - - GPUAxpyKernel<<>>(alpha, x.GetDevicePointer(), GetDevicePointer(), - nElm); - gpuErrChk(cudaPeekAtLastError()); +void CSysVector::DtHTransfer(bool trigger) const { + if (trigger) + gpuErrChk(cudaMemcpy((void*)(&vec_val[0]), (void*)d_vec_val, (sizeof(ScalarType) * nElm), cudaMemcpyDeviceToHost)); } template @@ -193,20 +137,12 @@ ScalarType CSysVector::GPUNorm() const { template void CSysVector::HtDTransfer(bool trigger) const; template void CSysVector::DtHTransfer(bool trigger) const; -template void CSysVector::GPUSetVal(su2mixedfloat val, bool trigger) const; -template void CSysVector::GPUCopy(const CSysVector& src) const; -template void CSysVector::GPUScale(su2mixedfloat alpha) const; -template void CSysVector::GPUAxpy(su2mixedfloat alpha, const CSysVector& x) const; template su2mixedfloat CSysVector::GPUDot(const CSysVector& other) const; template su2mixedfloat CSysVector::GPUNorm() const; #if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION) template void CSysVector::HtDTransfer(bool trigger) const; template void CSysVector::DtHTransfer(bool trigger) const; -template void CSysVector::GPUSetVal(passivedouble val, bool trigger) const; -template void CSysVector::GPUCopy(const CSysVector& src) const; -template void CSysVector::GPUScale(passivedouble alpha) const; -template void CSysVector::GPUAxpy(passivedouble alpha, const CSysVector& x) const; template passivedouble CSysVector::GPUDot(const CSysVector& other) const; template passivedouble CSysVector::GPUNorm() const; #endif @@ -214,10 +150,6 @@ template passivedouble CSysVector::GPUNorm() const; #ifdef CODI_REVERSE_TYPE template void CSysVector::HtDTransfer(bool trigger) const; template void CSysVector::DtHTransfer(bool trigger) const; -template void CSysVector::GPUSetVal(su2double val, bool trigger) const; -template void CSysVector::GPUCopy(const CSysVector& src) const; -template void CSysVector::GPUScale(su2double alpha) const; -template void CSysVector::GPUAxpy(su2double alpha, const CSysVector& x) const; template su2double CSysVector::GPUDot(const CSysVector& other) const; template su2double CSysVector::GPUNorm() const; #endif From 0ea444401bf902aa61a157812459f77fa8f4ab2f Mon Sep 17 00:00:00 2001 From: LwhJesse <256257451+LwhJesse@users.noreply.github.com> Date: Sat, 27 Jun 2026 23:45:50 +0800 Subject: [PATCH 07/10] Store CSysVector expressions as views for CUDA dispatch --- Common/include/linear_algebra/CSysVector.hpp | 168 +++++++++--------- .../linear_algebra/vector_expressions.hpp | 6 +- Common/src/linear_algebra/CSysSolveGPU.cu | 30 +++- Common/src/linear_algebra/CSysVectorGPU.cu | 12 ++ 4 files changed, 120 insertions(+), 96 deletions(-) diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index d67185156b9..0d279242856 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -39,8 +39,14 @@ #ifdef __CUDACC__ #include "GPUComms.cuh" +#define SU2_CUDA_HOST_DEVICE __host__ __device__ +#else +#define SU2_CUDA_HOST_DEVICE #endif +template +class CSysVector; + /*! * \brief OpenMP worksharing construct used in CSysVector for loops. * \note The loop will only run in parallel if methods are called from a @@ -67,41 +73,35 @@ namespace VecExpr { enum class DeviceAssignOp { Assign, Add, Subtract, Multiply, Divide }; +#ifdef __CUDACC__ +bool DeviceExpressionsEnabled(); +void SetDeviceExpressionsEnabled(bool enabled); +#else +inline bool DeviceExpressionsEnabled() { return false; } +#endif + template -class CDeviceVectorView : public CVecExpr, Scalar> { +class CVectorView : public CVecExpr, Scalar> { private: - Scalar* data = nullptr; - unsigned long size = 0; + const Scalar* data = nullptr; public: static constexpr bool StoreAsRef = false; - CDeviceVectorView(Scalar* data_, unsigned long size_) : data(data_), size(size_) {} + CVectorView(const Scalar* data_) : data(data_) {} + CVectorView(const CSysVector& vector); -#ifdef __CUDACC__ - __device__ const Scalar& operator[](size_t i) const { return data[i]; } + SU2_CUDA_HOST_DEVICE FORCEINLINE const Scalar& operator[](size_t i) const { return data[i]; } +}; - template - CDeviceVectorView& operator=(const CVecExpr& expr); - - CDeviceVectorView& operator=(Scalar val); - -#define MAKE_DEVICE_COMPOUND(OP) \ - template \ - CDeviceVectorView& operator OP(const CVecExpr& expr); - MAKE_DEVICE_COMPOUND(+=) - MAKE_DEVICE_COMPOUND(-=) - MAKE_DEVICE_COMPOUND(*=) - MAKE_DEVICE_COMPOUND(/=) -#undef MAKE_DEVICE_COMPOUND - -#define MAKE_DEVICE_SCALAR_COMPOUND(OP) CDeviceVectorView& operator OP(Scalar val); - MAKE_DEVICE_SCALAR_COMPOUND(+=) - MAKE_DEVICE_SCALAR_COMPOUND(-=) - MAKE_DEVICE_SCALAR_COMPOUND(*=) - MAKE_DEVICE_SCALAR_COMPOUND(/=) -#undef MAKE_DEVICE_SCALAR_COMPOUND -#endif +template +struct store_type> { + using type = CVectorView; +}; + +template +struct store_type> { + using type = CVectorView; }; #ifdef __CUDACC__ @@ -131,44 +131,6 @@ inline void AssignDeviceExpression(Scalar* data, unsigned long size, const CVecE DeviceAssignKernel<<>>(data, size, expr.derived()); gpuErrChk(cudaPeekAtLastError()); } - -template -template -CDeviceVectorView& CDeviceVectorView::operator=(const CVecExpr& expr) { - AssignDeviceExpression(data, size, expr); - return *this; -} - -template -CDeviceVectorView& CDeviceVectorView::operator=(Scalar val) { - AssignDeviceExpression(data, size, Bcast(val)); - return *this; -} - -#define MAKE_DEVICE_COMPOUND(OP, ASSIGN_OP) \ - template \ - template \ - CDeviceVectorView& CDeviceVectorView::operator OP(const CVecExpr& expr) { \ - AssignDeviceExpression(data, size, expr); \ - return *this; \ - } -MAKE_DEVICE_COMPOUND(+=, DeviceAssignOp::Add) -MAKE_DEVICE_COMPOUND(-=, DeviceAssignOp::Subtract) -MAKE_DEVICE_COMPOUND(*=, DeviceAssignOp::Multiply) -MAKE_DEVICE_COMPOUND(/=, DeviceAssignOp::Divide) -#undef MAKE_DEVICE_COMPOUND - -#define MAKE_DEVICE_SCALAR_COMPOUND(OP, ASSIGN_OP) \ - template \ - CDeviceVectorView& CDeviceVectorView::operator OP(Scalar val) { \ - AssignDeviceExpression(data, size, Bcast(val)); \ - return *this; \ - } -MAKE_DEVICE_SCALAR_COMPOUND(+=, DeviceAssignOp::Add) -MAKE_DEVICE_SCALAR_COMPOUND(-=, DeviceAssignOp::Subtract) -MAKE_DEVICE_SCALAR_COMPOUND(*=, DeviceAssignOp::Multiply) -MAKE_DEVICE_SCALAR_COMPOUND(/=, DeviceAssignOp::Divide) -#undef MAKE_DEVICE_SCALAR_COMPOUND #endif } // namespace VecExpr @@ -224,6 +186,23 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> } } + template + CSysVector& AssignDevice(const VecExpr::CVecExpr& expr) { +#ifdef __CUDACC__ + VecExpr::store_t stored_expr(expr.derived()); + VecExpr::AssignDeviceExpression(d_vec_val, nElm, stored_expr); +#endif + return *this; + } + + template + CSysVector& AssignDevice(ScalarType val) { +#ifdef __CUDACC__ + VecExpr::AssignDeviceExpression(d_vec_val, nElm, VecExpr::Bcast(val)); +#endif + return *this; + } + public: static constexpr bool StoreAsRef = true; /*! \brief Required by CVecExpr. */ @@ -369,10 +348,10 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> inline ScalarType* GetDevicePointer() const { return d_vec_val; } /*! - * \brief Return an expression-template view of the device values. + * \brief Return the pointer used by expression-template views. */ - inline VecExpr::CDeviceVectorView device() const { - return VecExpr::CDeviceVectorView(d_vec_val, nElm); + inline const ScalarType* GetExpressionPointer() const { + return VecExpr::DeviceExpressionsEnabled() ? d_vec_val : vec_val; } /*! @@ -431,6 +410,13 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] other - Another vector. */ CSysVector& operator=(const CSysVector& other) { +#ifdef __CUDACC__ + if (VecExpr::DeviceExpressionsEnabled()) { + VecExpr::CVectorView view(other); + VecExpr::AssignDeviceExpression(d_vec_val, nElm, view); + return *this; + } +#endif CSYSVEC_PARFOR for (auto i = 0ul; i < nElm; ++i) vec_val[i] = other.vec_val[i]; END_CSYSVEC_PARFOR @@ -441,25 +427,27 @@ class CSysVector : public VecExpr::CVecExpr, 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 \ - CSysVector& operator OP(const VecExpr::CVecExpr& 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, ASSIGN_OP) \ + CSysVector& operator OP(ScalarType val) { \ + if (VecExpr::DeviceExpressionsEnabled()) return AssignDevice(val); \ + CSYSVEC_PARFOR \ + for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val; \ + END_CSYSVEC_PARFOR \ + return *this; \ + } \ + template \ + CSysVector& operator OP(const VecExpr::CVecExpr& expr) { \ + if (VecExpr::DeviceExpressionsEnabled()) return AssignDevice(expr); \ + CSYSVEC_PARFOR \ + for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP expr.derived()[i]; \ + END_CSYSVEC_PARFOR \ + return *this; \ } - MAKE_COMPOUND(=) - MAKE_COMPOUND(+=) - MAKE_COMPOUND(-=) - MAKE_COMPOUND(*=) - MAKE_COMPOUND(/=) + MAKE_COMPOUND(=, VecExpr::DeviceAssignOp::Assign) + MAKE_COMPOUND(+=, VecExpr::DeviceAssignOp::Add) + MAKE_COMPOUND(-=, VecExpr::DeviceAssignOp::Subtract) + MAKE_COMPOUND(*=, VecExpr::DeviceAssignOp::Multiply) + MAKE_COMPOUND(/=, VecExpr::DeviceAssignOp::Divide) #undef MAKE_COMPOUND /*! @@ -632,5 +620,13 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> } }; +namespace VecExpr { + +template +CVectorView::CVectorView(const CSysVector& vector) : data(vector.GetExpressionPointer()) {} + +} // namespace VecExpr + #undef CSYSVEC_PARFOR #undef END_CSYSVEC_PARFOR +#undef SU2_CUDA_HOST_DEVICE diff --git a/Common/include/linear_algebra/vector_expressions.hpp b/Common/include/linear_algebra/vector_expressions.hpp index b2c6a50f5a7..f9b941a0bdb 100644 --- a/Common/include/linear_algebra/vector_expressions.hpp +++ b/Common/include/linear_algebra/vector_expressions.hpp @@ -107,7 +107,11 @@ struct add_lref_if { using type = remove_reference_t&; }; template -using store_t = typename add_lref_if::type; +struct store_type { + using type = typename add_lref_if::type; +}; +template +using store_t = typename store_type::type; /*--- Namespace from which the math function implementations come. ---*/ diff --git a/Common/src/linear_algebra/CSysSolveGPU.cu b/Common/src/linear_algebra/CSysSolveGPU.cu index d618cc223b5..a637da42eb4 100644 --- a/Common/src/linear_algebra/CSysSolveGPU.cu +++ b/Common/src/linear_algebra/CSysSolveGPU.cu @@ -45,6 +45,21 @@ class CGPUSolverBLASContextGuard { CGPUSolverBLASContextGuard& operator=(const CGPUSolverBLASContextGuard&) = delete; }; +class CDeviceExpressionContextGuard { + private: + bool previous = false; + + public: + CDeviceExpressionContextGuard() : previous(VecExpr::DeviceExpressionsEnabled()) { + VecExpr::SetDeviceExpressionsEnabled(true); + } + + ~CDeviceExpressionContextGuard() { VecExpr::SetDeviceExpressionsEnabled(previous); } + + CDeviceExpressionContextGuard(const CDeviceExpressionContextGuard&) = delete; + CDeviceExpressionContextGuard& operator=(const CDeviceExpressionContextGuard&) = delete; +}; + template class CDeviceFGMRESOps { public: @@ -55,7 +70,7 @@ class CDeviceFGMRESOps { void PrepareInput(const CSysVector& b, CSysVector& x, bool x_is_zero) const { b.HtDTransfer(); if (x_is_zero) { - x.device() = ScalarType(0); + x = ScalarType(0); } else { x.HtDTransfer(); } @@ -65,22 +80,18 @@ class CDeviceFGMRESOps { ScalarType Norm(const CSysVector& v) const { return v.GPUNorm(); } - void AssignNegative(CSysVector& dst, const CSysVector& src) const { - dst.device() = -src.device(); - } + void AssignNegative(CSysVector& dst, const CSysVector& src) const { dst = -src; } - void Subtract(CSysVector& dst, const CSysVector& src) const { - dst.device() -= src.device(); - } + void Subtract(CSysVector& dst, const CSysVector& src) const { dst -= src; } - void Divide(CSysVector& dst, ScalarType val) const { dst.device() /= val; } + void Divide(CSysVector& dst, ScalarType val) const { dst /= val; } ScalarType Dot(const CSysVector& lhs, const CSysVector& rhs) const { return lhs.GPUDot(rhs); } void AddScaled(CSysVector& dst, ScalarType alpha, const CSysVector& src) const { - dst.device() += alpha * src.device(); + dst += alpha * src; } void UpdateSolution(unsigned long n, const std::vector>& basis, @@ -101,6 +112,7 @@ unsigned long CSysSolve::FGMRES_LinSolver_GPU(const VectorType& b, V ScalarType tol, unsigned long m, ScalarType& residual, bool monitoring, const CConfig* config) const { CGPUSolverBLASContextGuard blas_context; + CDeviceExpressionContextGuard device_expression_context; return FGMRES_LinSolverImpl(b, x, mat_vec, precond, tol, m, residual, monitoring, config, CDeviceFGMRESOps()); } diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index dca72c0d76b..64f7f16504b 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -57,6 +57,18 @@ void ReleaseBlasHandle(cublasHandle_t handle, bool owns_handle, const char* dest } // namespace +namespace VecExpr { + +namespace { +thread_local bool device_expressions_enabled = false; +} // namespace + +bool DeviceExpressionsEnabled() { return device_expressions_enabled; } + +void SetDeviceExpressionsEnabled(bool enabled) { device_expressions_enabled = enabled; } + +} // namespace VecExpr + void SU2_GPU_BeginSolverBLASContext() { if (active_solver_blas_depth == 0) { cublasHandle_t handle = nullptr; From b36e87b399635603ae50b089168cb3ece3ab140e Mon Sep 17 00:00:00 2001 From: LwhJesse <256257451+LwhJesse@users.noreply.github.com> Date: Sun, 28 Jun 2026 20:22:29 +0800 Subject: [PATCH 08/10] Move FGMRES implementation back into CSysSolve --- Common/include/linear_algebra/CSysSolve.hpp | 14 - .../linear_algebra/CSysSolveFGMRES.inl | 267 ----------------- Common/include/linear_algebra/CSysVector.hpp | 89 ++++-- .../linear_algebra/vector_expressions.hpp | 3 + Common/src/linear_algebra/CSysSolve.cpp | 271 +++++++++++++++--- Common/src/linear_algebra/CSysSolveGPU.cu | 132 --------- Common/src/linear_algebra/CSysVector.cpp | 17 ++ Common/src/linear_algebra/CSysVectorGPU.cu | 55 ++++ Common/src/linear_algebra/meson.build | 2 +- 9 files changed, 385 insertions(+), 465 deletions(-) delete mode 100644 Common/include/linear_algebra/CSysSolveFGMRES.inl delete mode 100644 Common/src/linear_algebra/CSysSolveGPU.cu diff --git a/Common/include/linear_algebra/CSysSolve.hpp b/Common/include/linear_algebra/CSysSolve.hpp index af76e93c07d..2ea3cbf7df3 100644 --- a/Common/include/linear_algebra/CSysSolve.hpp +++ b/Common/include/linear_algebra/CSysSolve.hpp @@ -304,20 +304,6 @@ class CSysSolve { ScalarType& residual, bool monitoring, const CConfig* config, FgcrodrMode mode, unsigned long custom_m) const; - template - unsigned long FGMRES_LinSolverImpl(const VectorType& b, VectorType& x, const ProductType& mat_vec, - const PrecondType& precond, ScalarType tol, unsigned long m, ScalarType& residual, - bool monitoring, const CConfig* config, const FGMRESOps& ops) const; - - /*! - * \brief CUDA/GPU backend wrapper for the Flexible GMRES solver. - * \note The algorithmic implementation is shared with the host solver. This helper - * only provides GPU vector-operation dispatch and synchronization. - */ - unsigned long FGMRES_LinSolver_GPU(const VectorType& b, VectorType& x, const ProductType& mat_vec, - const PrecondType& precond, ScalarType tol, unsigned long m, ScalarType& residual, - bool monitoring, const CConfig* config) const; - /*! * \brief Creates the inner solver for nested preconditioning if the settings allow it. * \returns True if the inner solver can be used. diff --git a/Common/include/linear_algebra/CSysSolveFGMRES.inl b/Common/include/linear_algebra/CSysSolveFGMRES.inl deleted file mode 100644 index cef424a729b..00000000000 --- a/Common/include/linear_algebra/CSysSolveFGMRES.inl +++ /dev/null @@ -1,267 +0,0 @@ -/*! - * \file CSysSolveFGMRES.inl - * \brief Shared implementation of the Flexible GMRES linear solver. - * \author Jesse Li - * \version 8.5.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) - * - * SU2 is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * SU2 is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with SU2. If not, see . - */ - -namespace fgmres_impl { - -template -bool Orthogonalize(unsigned long i, su2matrix& H, std::vector>& V, - const FGMRESOps& ops) { - for (unsigned long k = 0; k <= i; k++) { - H(k, i) = ops.Dot(V[k], V[i + 1]); - ops.AddScaled(V[i + 1], -H(k, i), V[k]); - } - - for (unsigned long k = 0; k <= i; k++) { - const ScalarType dh = ops.Dot(V[k], V[i + 1]); - H(k, i) += dh; - ops.AddScaled(V[i + 1], -dh, V[k]); - } - - const ScalarType nrm = ops.Norm(V[i + 1]); - if (nrm <= ScalarType(0) || nrm != nrm) { - H(i + 1, i) = ScalarType(0); - return false; - } - - H(i + 1, i) = nrm; - ops.Divide(V[i + 1], nrm); - return true; -} - -} // namespace fgmres_impl - -template -template -unsigned long CSysSolve::FGMRES_LinSolverImpl(const VectorType& b, VectorType& x, - const ProductType& mat_vec, const PrecondType& precond, - ScalarType tol, unsigned long m, ScalarType& residual, - bool monitoring, const CConfig* config, - const FGMRESOps& ops) const { - SU2_ZONE_SCOPED - - const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE); - const bool flexible = !precond.IsIdentity(); - const bool nestedParallel = ops.NestedParallel(); - - /*--- Check the subspace size. ---*/ - - if (m < 1) { - SU2_MPI::Error("Number of linear solver iterations must be greater than 0.", CURRENT_FUNCTION); - } - - if (m > 1000) { - SU2_MPI::Error("FGMRES subspace is too large (>1000).", CURRENT_FUNCTION); - } - - /*--- Allocate if not allocated yet. ---*/ - - if (V.size() <= m || (flexible && Z.size() <= m)) { - BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - V.resize(m + 1); - for (auto& w : V) w.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); - if (flexible) { - Z.resize(m + 1); - for (auto& z : Z) z.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); - } - } - END_SU2_OMP_SAFE_GLOBAL_ACCESS - } - - /*--- Define various arrays. In parallel, each thread of each rank has and works - on its own thread, since calculations on these arrays are based on dot products - (reduced across all threads and ranks) all threads do the same computations. ---*/ - - su2vector g(m + 1), sn(m + 1), cs(m + 1), y(m); - g = ScalarType(0); - sn = ScalarType(0); - cs = ScalarType(0); - y = ScalarType(0); - - su2matrix H(m + 1, m); - H = ScalarType(0); - - ops.PrepareInput(b, x, xIsZero); - - /*--- Calculate the norm of the rhs vector. ---*/ - - ScalarType norm0 = ops.Norm(b); - - /*--- Calculate the initial residual (actually the negative residual) and compute its norm. ---*/ - - if (!xIsZero) { - mat_vec(x, V[0]); - ops.Subtract(V[0], b); - } else { - ops.AssignNegative(V[0], b); - } - - ScalarType beta = xIsZero ? norm0 : ops.Norm(V[0]); - - /*--- FGMRES shares V and Z with FGCRODR, reset deflation because we're breaking - * relations between subspaces that FGCRODR relies on for recycling. ---*/ - - SU2_OMP_MASTER - ResetDeflation(); - END_SU2_OMP_MASTER - - /*--- Set the norm to the initial initial residual value. ---*/ - - if (tol_type == LinearToleranceType::RELATIVE) norm0 = beta; - - if (beta < tol * norm0 || beta < eps) { - /*--- System is already solved. ---*/ - - if (masterRank) { - SU2_OMP_MASTER - cout << "CSysSolve::FGMRES(): system solved by initial guess." << endl; - END_SU2_OMP_MASTER - } - residual = beta; - ops.FinalizeOutput(x); - return 0; - } - - /*--- Normalize residual to get w_{0} (the negative sign is because w[0] - holds the negative residual, as mentioned above). ---*/ - - ops.Divide(V[0], -beta); - - /*--- Initialize the RHS of the reduced system. ---*/ - - g[0] = beta; - - /*--- Output header information including initial residual. ---*/ - - unsigned long i = 0; - if (monitoring && masterRank) { - SU2_OMP_MASTER { - WriteHeader("FGMRES", tol, beta); - WriteHistory(i, beta / norm0); - } - END_SU2_OMP_MASTER - } - - /*--- Loop over all search directions. ---*/ - - for (i = 0; i < m; i++) { - /*--- Check if solution has converged. ---*/ - - if (beta < tol * norm0) break; - - if (flexible) { - /*--- Precondition the CSysVector v[i] and store result in z[i]. ---*/ - - precond(V[i], Z[i]); - - /*--- Add to Krylov subspace. ---*/ - - mat_vec(Z[i], V[i + 1]); - } else { - mat_vec(V[i], V[i + 1]); - } - - /*--- Modified Gram-Schmidt orthogonalization. ---*/ - - bool orthog_ok; - if constexpr (FGMRESOps::UseSolverModGramSchmidt) { - if (nestedParallel) { - /*--- "omp parallel if" does not work well here. ---*/ - SU2_OMP_PARALLEL - orthog_ok = ModGramSchmidt(true, i, H, V); - END_SU2_OMP_PARALLEL - } else { - orthog_ok = ModGramSchmidt(false, i, H, V); - } - } else { - orthog_ok = fgmres_impl::Orthogonalize(i, H, V, ops); - } - - if (!orthog_ok) { - if (masterRank) { - SU2_OMP_MASTER - cout << "WARNING: FGMRES orthogonalization failed, linear solver diverged." << endl; - END_SU2_OMP_MASTER - } - break; - } - - /*--- Apply old Givens rotations to new column of the Hessenberg matrix then generate the - new Givens rotation matrix and apply it to the last two elements of H[:][i] and g. ---*/ - - for (unsigned long k = 0; k < i; k++) ApplyGivens(sn[k], cs[k], H(k, i), H(k + 1, i)); - GenerateGivens(H(i, i), H(i + 1, i), sn[i], cs[i]); - ApplyGivens(sn[i], cs[i], g[i], g[i + 1]); - - /*--- Set L2 norm of residual and check if solution has converged. ---*/ - - beta = fabs(g[i + 1]); - - /*--- Output the relative residual if necessary. ---*/ - - if (monitoring && masterRank && ((i + 1) % monitorFreq == 0)) { - SU2_OMP_MASTER - WriteHistory(i + 1, beta / norm0); - END_SU2_OMP_MASTER - } - } - - /*--- Solve the least-squares system and update solution. ---*/ - - SolveReduced(i, H, g, y); - - const auto& basis = flexible ? Z : V; - - ops.UpdateSolution(i, basis, y, x, nestedParallel); - ops.FinalizeOutput(x); - - /*--- Recalculate final (neg.) residual (this should be optional). ---*/ - - if (monitoring && config->GetComm_Level() == COMM_FULL) { - if (masterRank) { - SU2_OMP_MASTER - WriteFinalResidual("FGMRES", i, beta / norm0); - END_SU2_OMP_MASTER - } - - if (recomputeRes) { - mat_vec(x, V[0]); - ops.Subtract(V[0], b); - ScalarType res = ops.Norm(V[0]); - - if (fabs(res - beta) > tol * 10) { - if (masterRank) { - SU2_OMP_MASTER - WriteWarning(beta, res, tol); - END_SU2_OMP_MASTER - } - } - } - } - - residual = beta / norm0; - return i; -} diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 0d279242856..5caf153565b 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -73,7 +73,7 @@ namespace VecExpr { enum class DeviceAssignOp { Assign, Add, Subtract, Multiply, Divide }; -#ifdef __CUDACC__ +#ifdef HAVE_CUDA bool DeviceExpressionsEnabled(); void SetDeviceExpressionsEnabled(bool enabled); #else @@ -104,6 +104,43 @@ struct store_type> { using type = CVectorView; }; +template +struct is_device_assignable> : std::true_type {}; + +template +struct is_device_assignable> : std::true_type {}; + +template +struct is_device_assignable> : std::true_type {}; + +template +struct is_device_assignable> : std::true_type {}; + +template +struct is_device_assignable, Scalar>> : std::true_type {}; + +template +struct is_device_assignable, Bcast, Scalar>> : std::true_type {}; + +template +struct is_device_assignable< + add_, Bcast, Scalar>, mul_, Bcast, Scalar>, Scalar>> + : std::true_type {}; + +template +struct is_device_assignable< + add_, Bcast, Scalar>, mul_, Bcast, Scalar>, Scalar>, + mul_, Bcast, Scalar>, Scalar>> : std::true_type {}; + +template +struct is_device_assignable, Bcast, Scalar>, mul_, Bcast, Scalar>, Scalar>, + mul_, Bcast, Scalar>, Scalar>, + mul_, Bcast, Scalar>, Scalar>> : std::true_type {}; + +template +void AssignDeviceExpression(Scalar* data, unsigned long size, const CVecExpr& expr); + #ifdef __CUDACC__ template __global__ void DeviceAssignKernel(Scalar* data, unsigned long size, T expr) { @@ -188,7 +225,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> template CSysVector& AssignDevice(const VecExpr::CVecExpr& expr) { -#ifdef __CUDACC__ +#ifdef HAVE_CUDA VecExpr::store_t stored_expr(expr.derived()); VecExpr::AssignDeviceExpression(d_vec_val, nElm, stored_expr); #endif @@ -197,7 +234,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> template CSysVector& AssignDevice(ScalarType val) { -#ifdef __CUDACC__ +#ifdef HAVE_CUDA VecExpr::AssignDeviceExpression(d_vec_val, nElm, VecExpr::Bcast(val)); #endif return *this; @@ -410,7 +447,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \param[in] other - Another vector. */ CSysVector& operator=(const CSysVector& other) { -#ifdef __CUDACC__ +#ifdef HAVE_CUDA if (VecExpr::DeviceExpressionsEnabled()) { VecExpr::CVectorView view(other); VecExpr::AssignDeviceExpression(d_vec_val, nElm, view); @@ -427,21 +464,26 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> * \brief Compound assignement operations with scalars and expressions. * \param[in] val/expr - Scalar value or expression. */ -#define MAKE_COMPOUND(OP, ASSIGN_OP) \ - CSysVector& operator OP(ScalarType val) { \ - if (VecExpr::DeviceExpressionsEnabled()) return AssignDevice(val); \ - CSYSVEC_PARFOR \ - for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val; \ - END_CSYSVEC_PARFOR \ - return *this; \ - } \ - template \ - CSysVector& operator OP(const VecExpr::CVecExpr& expr) { \ - if (VecExpr::DeviceExpressionsEnabled()) return AssignDevice(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, ASSIGN_OP) \ + CSysVector& operator OP(ScalarType val) { \ + if (VecExpr::DeviceExpressionsEnabled()) return AssignDevice(val); \ + CSYSVEC_PARFOR \ + for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val; \ + END_CSYSVEC_PARFOR \ + return *this; \ + } \ + template \ + CSysVector& operator OP(const VecExpr::CVecExpr& expr) { \ + if (VecExpr::DeviceExpressionsEnabled()) { \ + using DeviceExpr = std::remove_cv_t>; \ + if constexpr (VecExpr::is_device_assignable::value) { \ + return AssignDevice(expr); \ + } \ + } \ + CSYSVEC_PARFOR \ + for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP expr.derived()[i]; \ + END_CSYSVEC_PARFOR \ + return *this; \ } MAKE_COMPOUND(=, VecExpr::DeviceAssignOp::Assign) MAKE_COMPOUND(+=, VecExpr::DeviceAssignOp::Add) @@ -462,6 +504,15 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ template ScalarType dot(const VecExpr::CVecExpr& expr) const { +#ifdef HAVE_CUDA + if (VecExpr::DeviceExpressionsEnabled()) { + using DeviceExpr = std::remove_cv_t>; + if constexpr (std::is_same_v) { + return GPUDot(expr.derived()); + } + } +#endif + /*--- All threads get the same "view" of the vectors. ---*/ SU2_OMP_BARRIER diff --git a/Common/include/linear_algebra/vector_expressions.hpp b/Common/include/linear_algebra/vector_expressions.hpp index f9b941a0bdb..f1f2639e22a 100644 --- a/Common/include/linear_algebra/vector_expressions.hpp +++ b/Common/include/linear_algebra/vector_expressions.hpp @@ -113,6 +113,9 @@ struct store_type { template using store_t = typename store_type::type; +template +struct is_device_assignable : std::false_type {}; + /*--- Namespace from which the math function implementations come. ---*/ #if defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE) diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index 1045ae68bbc..d001930aeaf 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -44,6 +44,11 @@ SU2_RESTORE_WARNING #include #include +#ifdef HAVE_CUDA +void SU2_GPU_BeginSolverBLASContext(); +void SU2_GPU_EndSolverBLASContext(); +#endif + namespace { /*! * \brief Epsilon used in CSysSolve depending on datatype to @@ -132,34 +137,56 @@ void LinearCombination(bool parallel, Ts&&... args) { } template -class CHostFGMRESOps { - public: - static constexpr bool UseSolverModGramSchmidt = true; - - bool NestedParallel() const { return !omp_in_parallel() && omp_get_max_threads() > 1; } - - void PrepareInput(const CSysVector&, CSysVector&, bool) const {} - - void FinalizeOutput(CSysVector&) const {} - - ScalarType Norm(const CSysVector& v) const { return v.norm(); } - - void AssignNegative(CSysVector& dst, const CSysVector& src) const { dst = -src; } - - void Subtract(CSysVector& dst, const CSysVector& src) const { dst -= src; } +class CDeviceSolveContextGuard { + private: + bool enabled = false; +#ifdef HAVE_CUDA + bool previous = false; + CSysVector* solution = nullptr; +#endif - void Divide(CSysVector& dst, ScalarType val) const { dst /= val; } + public: + CDeviceSolveContextGuard(bool enabled_, const CSysVector& b, CSysVector& x, bool x_is_zero) + : enabled(enabled_) { +#ifdef HAVE_CUDA + if (enabled) { + solution = &x; + SU2_GPU_BeginSolverBLASContext(); + previous = VecExpr::DeviceExpressionsEnabled(); + VecExpr::SetDeviceExpressionsEnabled(true); + b.HtDTransfer(); + if (x_is_zero) { + x = ScalarType(0); + } else { + x.HtDTransfer(); + } + } +#else + if (enabled) { + SU2_MPI::Error( + "\nError in launching FGMRES solver\nENABLE_CUDA is set to YES\nPlease compile with CUDA options enabled " + "in Meson to access GPU Functions", + CURRENT_FUNCTION); + } +#endif + } - void UpdateSolution(unsigned long n, const std::vector>& basis, const su2vector& y, - CSysVector& x, bool nestedParallel) const { - LinearCombination(nestedParallel, n, basis, y, x, true); + ~CDeviceSolveContextGuard() { +#ifdef HAVE_CUDA + if (enabled) { + solution->DtHTransfer(); + VecExpr::SetDeviceExpressionsEnabled(previous); + SU2_GPU_EndSolverBLASContext(); + } +#endif } + + CDeviceSolveContextGuard(const CDeviceSolveContextGuard&) = delete; + CDeviceSolveContextGuard& operator=(const CDeviceSolveContextGuard&) = delete; }; } // namespace -#include "../../include/linear_algebra/CSysSolveFGMRES.inl" - template CSysSolve::CSysSolve(LINEAR_SOLVER_MODE linear_solver_mode) : eps(linSolEpsilon()), @@ -448,19 +475,199 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVectorGetCUDA()) { -#ifdef HAVE_CUDA - return FGMRES_LinSolver_GPU(b, x, mat_vec, precond, tol, m, residual, monitoring, config); -#else - SU2_MPI::Error( - "\nError in launching FGMRES solver\nENABLE_CUDA is set to YES\nPlease compile with CUDA options enabled " - "in Meson to access GPU Functions", - CURRENT_FUNCTION); -#endif + const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE); + const bool flexible = !precond.IsIdentity(); + const bool nestedParallel = !omp_in_parallel() && omp_get_max_threads() > 1; + CDeviceSolveContextGuard device_context(config->GetCUDA(), b, x, xIsZero); + + /*--- Check the subspace size. ---*/ + + if (m < 1) { + SU2_MPI::Error("Number of linear solver iterations must be greater than 0.", CURRENT_FUNCTION); + } + + if (m > 1000) { + SU2_MPI::Error("FGMRES subspace is too large (>1000).", CURRENT_FUNCTION); + } + + /*--- Allocate if not allocated yet. ---*/ + + if (V.size() <= m || (flexible && Z.size() <= m)) { + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + V.resize(m + 1); + for (auto& w : V) w.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); + if (flexible) { + Z.resize(m + 1); + for (auto& z : Z) z.Initialize(x.GetNBlk(), x.GetNBlkDomain(), x.GetNVar(), nullptr); + } + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + } + + /*--- Define various arrays. In parallel, each thread of each rank has and works + on its own thread, since calculations on these arrays are based on dot products + (reduced across all threads and ranks) all threads do the same computations. ---*/ + + su2vector g(m + 1), sn(m + 1), cs(m + 1), y(m); + g = ScalarType(0); + sn = ScalarType(0); + cs = ScalarType(0); + y = ScalarType(0); + + su2matrix H(m + 1, m); + H = ScalarType(0); + + /*--- Calculate the norm of the rhs vector. ---*/ + + ScalarType norm0 = b.norm(); + + /*--- Calculate the initial residual (actually the negative residual) and compute its norm. ---*/ + + if (!xIsZero) { + mat_vec(x, V[0]); + V[0] -= b; + } else { + V[0] = -b; + } + + ScalarType beta = xIsZero ? norm0 : V[0].norm(); + + /*--- FGMRES shares V and Z with FGCRODR, reset deflation because we're breaking + * relations between subspaces that FGCRODR relies on for recycling. ---*/ + + SU2_OMP_MASTER + ResetDeflation(); + END_SU2_OMP_MASTER + + /*--- Set the norm to the initial initial residual value. ---*/ + + if (tol_type == LinearToleranceType::RELATIVE) norm0 = beta; + + if (beta < tol * norm0 || beta < eps) { + /*--- System is already solved. ---*/ + + if (masterRank) { + SU2_OMP_MASTER + cout << "CSysSolve::FGMRES(): system solved by initial guess." << endl; + END_SU2_OMP_MASTER + } + residual = beta; + return 0; + } + + /*--- Normalize residual to get w_{0} (the negative sign is because w[0] + holds the negative residual, as mentioned above). ---*/ + + V[0] /= -beta; + + /*--- Initialize the RHS of the reduced system. ---*/ + + g[0] = beta; + + /*--- Output header information including initial residual. ---*/ + + unsigned long i = 0; + if (monitoring && masterRank) { + SU2_OMP_MASTER { + WriteHeader("FGMRES", tol, beta); + WriteHistory(i, beta / norm0); + } + END_SU2_OMP_MASTER + } + + /*--- Loop over all search directions. ---*/ + + for (i = 0; i < m; i++) { + /*--- Check if solution has converged. ---*/ + + if (beta < tol * norm0) break; + + if (flexible) { + /*--- Precondition the CSysVector v[i] and store result in z[i]. ---*/ + + precond(V[i], Z[i]); + + /*--- Add to Krylov subspace. ---*/ + + mat_vec(Z[i], V[i + 1]); + } else { + mat_vec(V[i], V[i + 1]); + } + + /*--- Modified Gram-Schmidt orthogonalization. ---*/ + + bool orthog_ok; + if (nestedParallel) { + /*--- "omp parallel if" does not work well here. ---*/ + SU2_OMP_PARALLEL + orthog_ok = ModGramSchmidt(true, i, H, V); + END_SU2_OMP_PARALLEL + } else { + orthog_ok = ModGramSchmidt(false, i, H, V); + } + + if (!orthog_ok) { + if (masterRank) { + SU2_OMP_MASTER + cout << "WARNING: FGMRES orthogonalization failed, linear solver diverged." << endl; + END_SU2_OMP_MASTER + } + break; + } + + /*--- Apply old Givens rotations to new column of the Hessenberg matrix then generate the + new Givens rotation matrix and apply it to the last two elements of H[:][i] and g. ---*/ + + for (unsigned long k = 0; k < i; k++) ApplyGivens(sn[k], cs[k], H(k, i), H(k + 1, i)); + GenerateGivens(H(i, i), H(i + 1, i), sn[i], cs[i]); + ApplyGivens(sn[i], cs[i], g[i], g[i + 1]); + + /*--- Set L2 norm of residual and check if solution has converged. ---*/ + + beta = fabs(g[i + 1]); + + /*--- Output the relative residual if necessary. ---*/ + + if (monitoring && masterRank && ((i + 1) % monitorFreq == 0)) { + SU2_OMP_MASTER + WriteHistory(i + 1, beta / norm0); + END_SU2_OMP_MASTER + } } - return FGMRES_LinSolverImpl(b, x, mat_vec, precond, tol, m, residual, monitoring, config, - CHostFGMRESOps()); + /*--- Solve the least-squares system and update solution. ---*/ + + SolveReduced(i, H, g, y); + + const auto& basis = flexible ? Z : V; + LinearCombination(nestedParallel, i, basis, y, x, true); + + /*--- Recalculate final (neg.) residual (this should be optional). ---*/ + + if (monitoring && config->GetComm_Level() == COMM_FULL) { + if (masterRank) { + SU2_OMP_MASTER + WriteFinalResidual("FGMRES", i, beta / norm0); + END_SU2_OMP_MASTER + } + + if (recomputeRes) { + mat_vec(x, V[0]); + V[0] -= b; + ScalarType res = V[0].norm(); + + if (fabs(res - beta) > tol * 10) { + if (masterRank) { + SU2_OMP_MASTER + WriteWarning(beta, res, tol); + END_SU2_OMP_MASTER + } + } + } + } + + residual = beta / norm0; + return i; } template diff --git a/Common/src/linear_algebra/CSysSolveGPU.cu b/Common/src/linear_algebra/CSysSolveGPU.cu deleted file mode 100644 index a637da42eb4..00000000000 --- a/Common/src/linear_algebra/CSysSolveGPU.cu +++ /dev/null @@ -1,132 +0,0 @@ -/*! - * \file CSysSolveGPU.cu - * \brief CUDA/GPU backend dispatch for Krylov solvers. - * \author Jesse Li - * \version 8.5.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) - * - * SU2 is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * SU2 is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with SU2. If not, see . - */ - -#include "../../include/linear_algebra/CSysSolve.hpp" -#include "../../include/linear_algebra/CMatrixVectorProduct.hpp" -#include "../../include/linear_algebra/CPreconditioner.hpp" -#include "../../include/linear_algebra/GPUComms.cuh" - -void SU2_GPU_BeginSolverBLASContext(); -void SU2_GPU_EndSolverBLASContext(); - -namespace { - -class CGPUSolverBLASContextGuard { - public: - CGPUSolverBLASContextGuard() { SU2_GPU_BeginSolverBLASContext(); } - - ~CGPUSolverBLASContextGuard() { SU2_GPU_EndSolverBLASContext(); } - - CGPUSolverBLASContextGuard(const CGPUSolverBLASContextGuard&) = delete; - CGPUSolverBLASContextGuard& operator=(const CGPUSolverBLASContextGuard&) = delete; -}; - -class CDeviceExpressionContextGuard { - private: - bool previous = false; - - public: - CDeviceExpressionContextGuard() : previous(VecExpr::DeviceExpressionsEnabled()) { - VecExpr::SetDeviceExpressionsEnabled(true); - } - - ~CDeviceExpressionContextGuard() { VecExpr::SetDeviceExpressionsEnabled(previous); } - - CDeviceExpressionContextGuard(const CDeviceExpressionContextGuard&) = delete; - CDeviceExpressionContextGuard& operator=(const CDeviceExpressionContextGuard&) = delete; -}; - -template -class CDeviceFGMRESOps { - public: - static constexpr bool UseSolverModGramSchmidt = false; - - bool NestedParallel() const { return false; } - - void PrepareInput(const CSysVector& b, CSysVector& x, bool x_is_zero) const { - b.HtDTransfer(); - if (x_is_zero) { - x = ScalarType(0); - } else { - x.HtDTransfer(); - } - } - - void FinalizeOutput(CSysVector& x) const { x.DtHTransfer(); } - - ScalarType Norm(const CSysVector& v) const { return v.GPUNorm(); } - - void AssignNegative(CSysVector& dst, const CSysVector& src) const { dst = -src; } - - void Subtract(CSysVector& dst, const CSysVector& src) const { dst -= src; } - - void Divide(CSysVector& dst, ScalarType val) const { dst /= val; } - - ScalarType Dot(const CSysVector& lhs, const CSysVector& rhs) const { - return lhs.GPUDot(rhs); - } - - void AddScaled(CSysVector& dst, ScalarType alpha, const CSysVector& src) const { - dst += alpha * src; - } - - void UpdateSolution(unsigned long n, const std::vector>& basis, - const su2vector& y, CSysVector& x, bool) const { - for (unsigned long k = 0; k < n; k++) { - AddScaled(x, y[k], basis[k]); - } - } -}; - -} // namespace - -#include "../../include/linear_algebra/CSysSolveFGMRES.inl" - -template -unsigned long CSysSolve::FGMRES_LinSolver_GPU(const VectorType& b, VectorType& x, - const ProductType& mat_vec, const PrecondType& precond, - ScalarType tol, unsigned long m, ScalarType& residual, - bool monitoring, const CConfig* config) const { - CGPUSolverBLASContextGuard blas_context; - CDeviceExpressionContextGuard device_expression_context; - return FGMRES_LinSolverImpl(b, x, mat_vec, precond, tol, m, residual, monitoring, config, - CDeviceFGMRESOps()); -} - -/*--- Explicit instantiations for the GPU backend wrapper. - * Keep these aligned with the scalar types instantiated for CSysSolve. ---*/ -template unsigned long CSysSolve::FGMRES_LinSolver_GPU( - const CSysVector& b, CSysVector& x, - const CMatrixVectorProduct& mat_vec, const CPreconditioner& precond, - su2mixedfloat tol, unsigned long m, su2mixedfloat& residual, bool monitoring, const CConfig* config) const; - -#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION) -template unsigned long CSysSolve::FGMRES_LinSolver_GPU( - const CSysVector& b, CSysVector& x, - const CMatrixVectorProduct& mat_vec, const CPreconditioner& precond, - passivedouble tol, unsigned long m, passivedouble& residual, bool monitoring, const CConfig* config) const; -#endif diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp index f7df34633a0..27f07ce68b2 100644 --- a/Common/src/linear_algebra/CSysVector.cpp +++ b/Common/src/linear_algebra/CSysVector.cpp @@ -78,6 +78,23 @@ const su2matrix& CSysVector::multiDot(const std::vector< if (n == 0 || m == 0) return shared; +#ifdef HAVE_CUDA + if (VecExpr::DeviceExpressionsEnabled()) { + SU2_OMP_MASTER { + shared.resize(n, m); + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < m; ++j) { + shared(i, j) = V[i0 + i].GPUDot(W[j]); + } + } + } + END_SU2_OMP_MASTER + + SU2_OMP_BARRIER + return shared; + } +#endif + SU2_OMP_BARRIER const size_t size = V[0].nElmDomain; diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 64f7f16504b..7b06efbb7a0 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -147,6 +147,61 @@ ScalarType CSysVector::GPUNorm() const { return sqrt(GPUDot(*this)); } +namespace { + +template +using DeviceView = VecExpr::CVectorView; + +template +using DeviceBcast = VecExpr::Bcast; + +template +using DeviceNeg = VecExpr::minus_, Scalar>; + +template +using DeviceScale = VecExpr::mul_, DeviceBcast, Scalar>; + +template +using DeviceScale2 = VecExpr::add_, DeviceScale, Scalar>; + +template +using DeviceScale3 = VecExpr::add_, DeviceScale, Scalar>; + +template +using DeviceScale4 = VecExpr::add_, DeviceScale, Scalar>; + +} // namespace + +#define INSTANTIATE_DEVICE_ASSIGN(SCALAR, OP, EXPR) \ + template void VecExpr::AssignDeviceExpression>( \ + SCALAR*, unsigned long, const VecExpr::CVecExpr, SCALAR>&) + +#define INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, EXPR) \ + INSTANTIATE_DEVICE_ASSIGN(SCALAR, Assign, EXPR); \ + INSTANTIATE_DEVICE_ASSIGN(SCALAR, Add, EXPR); \ + INSTANTIATE_DEVICE_ASSIGN(SCALAR, Subtract, EXPR); \ + INSTANTIATE_DEVICE_ASSIGN(SCALAR, Multiply, EXPR); \ + INSTANTIATE_DEVICE_ASSIGN(SCALAR, Divide, EXPR) + +#define INSTANTIATE_DEVICE_ASSIGN_SCALAR(SCALAR) \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceBcast); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceView); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceNeg); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale2); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale3); \ + INSTANTIATE_DEVICE_ASSIGN_EXPR(SCALAR, DeviceScale4) + +INSTANTIATE_DEVICE_ASSIGN_SCALAR(su2mixedfloat); + +#if defined(USE_MIXED_PRECISION) && !defined(USE_SINGLE_PRECISION) +INSTANTIATE_DEVICE_ASSIGN_SCALAR(passivedouble); +#endif + +#undef INSTANTIATE_DEVICE_ASSIGN_SCALAR +#undef INSTANTIATE_DEVICE_ASSIGN_EXPR +#undef INSTANTIATE_DEVICE_ASSIGN + template void CSysVector::HtDTransfer(bool trigger) const; template void CSysVector::DtHTransfer(bool trigger) const; template su2mixedfloat CSysVector::GPUDot(const CSysVector& other) const; diff --git a/Common/src/linear_algebra/meson.build b/Common/src/linear_algebra/meson.build index 29bb63cd680..0142e24ec4d 100644 --- a/Common/src/linear_algebra/meson.build +++ b/Common/src/linear_algebra/meson.build @@ -6,5 +6,5 @@ common_src += files(['CSysSolve_b.cpp', 'blas_structure.cpp']) if get_option('enable-cuda') - common_src += files(['CSysMatrixGPU.cu', 'CSysVectorGPU.cu', 'CSysSolveGPU.cu', 'CSysPreconditionerGPU.cu',]) + common_src += files(['CSysMatrixGPU.cu', 'CSysVectorGPU.cu', 'CSysPreconditionerGPU.cu',]) endif From 059c029e003ad7fc752a7a73c039bb4f43562062 Mon Sep 17 00:00:00 2001 From: LwhJesse <256257451+LwhJesse@users.noreply.github.com> Date: Sun, 28 Jun 2026 20:52:09 +0800 Subject: [PATCH 09/10] Keep FGMRES solver agnostic to GPU backend --- .../linear_algebra/CMatrixVectorProduct.hpp | 3 +- Common/include/linear_algebra/CSysVector.hpp | 89 +++++++- Common/src/linear_algebra/CSysMatrixGPU.cu | 4 + .../linear_algebra/CSysPreconditionerGPU.cu | 4 + Common/src/linear_algebra/CSysSolve.cpp | 196 +++++++----------- Common/src/linear_algebra/CSysVector.cpp | 5 + Common/src/linear_algebra/CSysVectorGPU.cu | 99 +++++++-- 7 files changed, 253 insertions(+), 147 deletions(-) diff --git a/Common/include/linear_algebra/CMatrixVectorProduct.hpp b/Common/include/linear_algebra/CMatrixVectorProduct.hpp index a0cecaa63d7..e6829461f47 100644 --- a/Common/include/linear_algebra/CMatrixVectorProduct.hpp +++ b/Common/include/linear_algebra/CMatrixVectorProduct.hpp @@ -72,6 +72,7 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct { const CSysMatrix& matrix; /*!< \brief pointer to matrix that defines the product. */ CGeometry* geometry; /*!< \brief geometry associated with the matrix. */ const CConfig* config; /*!< \brief config of the problem. */ + VecExpr::CDeviceExpressionContext device_context; mutable bool matrix_uploaded = false; /*!< \brief Upload the matrix lazily on the first actual GPU matvec. */ public: @@ -83,7 +84,7 @@ class CSysMatrixVectorProduct final : public CMatrixVectorProduct { */ inline CSysMatrixVectorProduct(const CSysMatrix& matrix_ref, CGeometry* geometry_ref, const CConfig* config_ref) - : matrix(matrix_ref), geometry(geometry_ref), config(config_ref) {} + : matrix(matrix_ref), geometry(geometry_ref), config(config_ref), device_context(config_ref->GetCUDA()) {} /*! * \note This class cannot be default constructed as that would leave us with invalid pointers. diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index 5caf153565b..e8fcd624cbd 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -75,11 +75,34 @@ enum class DeviceAssignOp { Assign, Add, Subtract, Multiply, Divide }; #ifdef HAVE_CUDA bool DeviceExpressionsEnabled(); -void SetDeviceExpressionsEnabled(bool enabled); +void BeginDeviceExpressionContext(bool enabled); +void EndDeviceExpressionContext(); +void FlushDeviceExpressionContext(); +unsigned long CurrentDeviceExpressionContextId(); +void RegisterDeviceModifiedVector(const void* vector, void (*sync)(const void*)); +void UnregisterDeviceModifiedVector(const void* vector); #else inline bool DeviceExpressionsEnabled() { return false; } +inline void BeginDeviceExpressionContext(bool) {} +inline void EndDeviceExpressionContext() {} +inline void FlushDeviceExpressionContext() {} +inline unsigned long CurrentDeviceExpressionContextId() { return 0; } +inline void RegisterDeviceModifiedVector(const void*, void (*)(const void*)) {} +inline void UnregisterDeviceModifiedVector(const void*) {} #endif +class CDeviceExpressionContext { + private: + bool enabled = false; + + public: + explicit CDeviceExpressionContext(bool enabled_) : enabled(enabled_) { BeginDeviceExpressionContext(enabled); } + ~CDeviceExpressionContext() { EndDeviceExpressionContext(); } + + CDeviceExpressionContext(const CDeviceExpressionContext&) = delete; + CDeviceExpressionContext& operator=(const CDeviceExpressionContext&) = delete; +}; + template class CVectorView : public CVecExpr, Scalar> { private: @@ -190,6 +213,9 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> unsigned long 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. */ + mutable bool device_data_valid = false; + mutable bool host_data_valid = true; + mutable unsigned long device_context_id = 0; #ifdef HAVE_OMP mutable std::unique_ptr @@ -226,8 +252,10 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> template CSysVector& AssignDevice(const VecExpr::CVecExpr& expr) { #ifdef HAVE_CUDA + if constexpr (Op != VecExpr::DeviceAssignOp::Assign) EnsureDeviceData(); VecExpr::store_t stored_expr(expr.derived()); VecExpr::AssignDeviceExpression(d_vec_val, nElm, stored_expr); + MarkDeviceDataModified(); #endif return *this; } @@ -235,11 +263,18 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> template CSysVector& AssignDevice(ScalarType val) { #ifdef HAVE_CUDA + if constexpr (Op != VecExpr::DeviceAssignOp::Assign) EnsureDeviceData(); VecExpr::AssignDeviceExpression(d_vec_val, nElm, VecExpr::Bcast(val)); + MarkDeviceDataModified(); #endif return *this; } + void MarkHostDataModified() const { + host_data_valid = true; + device_data_valid = false; + } + public: static constexpr bool StoreAsRef = true; /*! \brief Required by CVecExpr. */ @@ -303,6 +338,9 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> std::swap(omp_chunk_size, other.omp_chunk_size); std::swap(vec_val, other.vec_val); std::swap(d_vec_val, other.d_vec_val); + std::swap(device_data_valid, other.device_data_valid); + std::swap(host_data_valid, other.host_data_valid); + std::swap(device_context_id, other.device_context_id); std::swap(nElm, other.nElm); std::swap(nElmDomain, other.nElmDomain); std::swap(nVar, other.nVar); @@ -350,6 +388,8 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> CSYSVEC_PARFOR for (auto i = 0ul; i < nElm; i++) vec_val[i] = SU2_TYPE::GetValue(other[i]); END_CSYSVEC_PARFOR + + MarkHostDataModified(); } /*! @@ -364,6 +404,39 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> */ void DtHTransfer(bool trigger = true) const; + void EnsureDeviceData() const { +#ifdef HAVE_CUDA + const auto context_id = VecExpr::CurrentDeviceExpressionContextId(); + if (VecExpr::DeviceExpressionsEnabled() && (!device_data_valid || device_context_id != context_id)) { + HtDTransfer(); + device_data_valid = true; + host_data_valid = true; + device_context_id = context_id; + } +#endif + } + + void MarkDeviceDataModified() const { +#ifdef HAVE_CUDA + if (VecExpr::DeviceExpressionsEnabled()) { + device_data_valid = true; + host_data_valid = false; + device_context_id = VecExpr::CurrentDeviceExpressionContextId(); + VecExpr::RegisterDeviceModifiedVector( + this, [](const void* vector) { static_cast(vector)->SyncHostFromDevice(); }); + } +#endif + } + + void SyncHostFromDevice() const { +#ifdef HAVE_CUDA + if (device_data_valid && !host_data_valid) { + DtHTransfer(); + host_data_valid = true; + } +#endif + } + /*! * \brief Dot product between this vector and another vector on the device. * \note Explicit GPU helper for solver-side reductions. @@ -449,14 +522,17 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> CSysVector& operator=(const CSysVector& other) { #ifdef HAVE_CUDA if (VecExpr::DeviceExpressionsEnabled()) { + other.EnsureDeviceData(); VecExpr::CVectorView view(other); VecExpr::AssignDeviceExpression(d_vec_val, nElm, view); + MarkDeviceDataModified(); return *this; } #endif CSYSVEC_PARFOR for (auto i = 0ul; i < nElm; ++i) vec_val[i] = other.vec_val[i]; END_CSYSVEC_PARFOR + MarkHostDataModified(); return *this; } @@ -470,6 +546,7 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> CSYSVEC_PARFOR \ for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val; \ END_CSYSVEC_PARFOR \ + MarkHostDataModified(); \ return *this; \ } \ template \ @@ -479,10 +556,12 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> if constexpr (VecExpr::is_device_assignable::value) { \ return AssignDevice(expr); \ } \ + VecExpr::FlushDeviceExpressionContext(); \ } \ CSYSVEC_PARFOR \ for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP expr.derived()[i]; \ END_CSYSVEC_PARFOR \ + MarkHostDataModified(); \ return *this; \ } MAKE_COMPOUND(=, VecExpr::DeviceAssignOp::Assign) @@ -508,8 +587,11 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> if (VecExpr::DeviceExpressionsEnabled()) { using DeviceExpr = std::remove_cv_t>; if constexpr (std::is_same_v) { + EnsureDeviceData(); + expr.derived().EnsureDeviceData(); return GPUDot(expr.derived()); } + VecExpr::FlushDeviceExpressionContext(); } #endif @@ -674,7 +756,10 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> namespace VecExpr { template -CVectorView::CVectorView(const CSysVector& vector) : data(vector.GetExpressionPointer()) {} +CVectorView::CVectorView(const CSysVector& vector) { + vector.EnsureDeviceData(); + data = vector.GetExpressionPointer(); +} } // namespace VecExpr diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu index b8bb16ce5a9..35e3c644496 100644 --- a/Common/src/linear_algebra/CSysMatrixGPU.cu +++ b/Common/src/linear_algebra/CSysMatrixGPU.cu @@ -108,6 +108,8 @@ void CSysMatrix::HtDTransfer(bool trigger) const template void CSysMatrix::GPUMatrixVectorProduct(const CSysVector& vec, CSysVector& prod, CGeometry* geometry, const CConfig* config) const { + vec.EnsureDeviceData(); + if (nVar != nEqn) { SU2_MPI::Error("CUDA CSysMatrix matvec with cuSPARSE BSR requires square blocks.", CURRENT_FUNCTION); } @@ -171,6 +173,8 @@ void CSysMatrix::GPUMatrixVectorProduct(const CSysVector cusparseErrChk(cusparseDestroyDnVec(vecY)); cusparseErrChk(cusparseDestroyDnVec(vecX)); + + prod.MarkDeviceDataModified(); } template void CSysMatrix::HtDTransfer(bool trigger) const; template void CSysMatrix::GPUMatrixVectorProduct(const CSysVector& vec, diff --git a/Common/src/linear_algebra/CSysPreconditionerGPU.cu b/Common/src/linear_algebra/CSysPreconditionerGPU.cu index adee93264eb..4241a0fd538 100644 --- a/Common/src/linear_algebra/CSysPreconditionerGPU.cu +++ b/Common/src/linear_algebra/CSysPreconditionerGPU.cu @@ -83,6 +83,8 @@ void CSysMatrix::ComputeJacobiPreconditionerGPU(const CSysVector::ComputeJacobiPreconditionerGPU(const CSysVector>>(d_invM, vec.GetDevicePointer(), prod.GetDevicePointer(), nPointDomain, nVar); gpuErrChk(cudaPeekAtLastError()); + + prod.MarkDeviceDataModified(); } template void CSysMatrix::BuildJacobiPreconditionerGPU(); diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index d001930aeaf..f6deffea4df 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -41,14 +41,8 @@ SU2_RESTORE_WARNING #include #include #include -#include #include -#ifdef HAVE_CUDA -void SU2_GPU_BeginSolverBLASContext(); -void SU2_GPU_EndSolverBLASContext(); -#endif - namespace { /*! * \brief Epsilon used in CSysSolve depending on datatype to @@ -136,55 +130,6 @@ void LinearCombination(bool parallel, Ts&&... args) { } } -template -class CDeviceSolveContextGuard { - private: - bool enabled = false; -#ifdef HAVE_CUDA - bool previous = false; - CSysVector* solution = nullptr; -#endif - - public: - CDeviceSolveContextGuard(bool enabled_, const CSysVector& b, CSysVector& x, bool x_is_zero) - : enabled(enabled_) { -#ifdef HAVE_CUDA - if (enabled) { - solution = &x; - SU2_GPU_BeginSolverBLASContext(); - previous = VecExpr::DeviceExpressionsEnabled(); - VecExpr::SetDeviceExpressionsEnabled(true); - b.HtDTransfer(); - if (x_is_zero) { - x = ScalarType(0); - } else { - x.HtDTransfer(); - } - } -#else - if (enabled) { - SU2_MPI::Error( - "\nError in launching FGMRES solver\nENABLE_CUDA is set to YES\nPlease compile with CUDA options enabled " - "in Meson to access GPU Functions", - CURRENT_FUNCTION); - } -#endif - } - - ~CDeviceSolveContextGuard() { -#ifdef HAVE_CUDA - if (enabled) { - solution->DtHTransfer(); - VecExpr::SetDeviceExpressionsEnabled(previous); - SU2_GPU_EndSolverBLASContext(); - } -#endif - } - - CDeviceSolveContextGuard(const CDeviceSolveContextGuard&) = delete; - CDeviceSolveContextGuard& operator=(const CDeviceSolveContextGuard&) = delete; -}; - } // namespace template @@ -478,7 +423,6 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector 1; - CDeviceSolveContextGuard device_context(config->GetCUDA(), b, x, xIsZero); /*--- Check the subspace size. ---*/ @@ -1517,84 +1461,86 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con HandleTemporariesIn(LinSysRes, LinSysSol); - auto mat_vec = CSysMatrixVectorProduct(Jacobian, geometry, config); - - /*--- Build preconditioner. ---*/ - - const auto kindPrec = static_cast(KindPrecond); - auto* normal_prec = CPreconditioner::Create(kindPrec, Jacobian, geometry, config); - normal_prec->Build(); - - CPreconditioner* nested_prec = nullptr; - if (nested) { - auto f = [&](const CSysVector& u, CSysVector& v) { - /*--- Initialize to 0 to be safe. ---*/ - v = ScalarType{}; - ScalarType res{}; - /*--- Handle other types here if desired but do not call Solve because - * that will create issues with the AD external function. ---*/ - if (config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::BCGSTAB) { - inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, sqrt(SolverTol), MaxIter, res, false, config); - } else { - const auto smooth_iter = static_cast(std::round(fmax(2, sqrt(MaxIter)))); - inner_solver->Smoother_LinSolver(u, v, mat_vec, *normal_prec, 0, smooth_iter, res, false, config); - } - }; - nested_prec = new CAbstractPreconditioner(f); - } - const auto* precond = nested ? nested_prec : normal_prec; + { + auto mat_vec = CSysMatrixVectorProduct(Jacobian, geometry, config); + + /*--- Build preconditioner. ---*/ + + const auto kindPrec = static_cast(KindPrecond); + auto* normal_prec = CPreconditioner::Create(kindPrec, Jacobian, geometry, config); + normal_prec->Build(); + + CPreconditioner* nested_prec = nullptr; + if (nested) { + auto f = [&](const CSysVector& u, CSysVector& v) { + /*--- Initialize to 0 to be safe. ---*/ + v = ScalarType{}; + ScalarType res{}; + /*--- Handle other types here if desired but do not call Solve because + * that will create issues with the AD external function. ---*/ + if (config->GetKind_Linear_Solver_Inner() == LINEAR_SOLVER_INNER::BCGSTAB) { + inner_solver->BCGSTAB_LinSolver(u, v, mat_vec, *normal_prec, sqrt(SolverTol), MaxIter, res, false, config); + } else { + const auto smooth_iter = static_cast(std::round(fmax(2, sqrt(MaxIter)))); + inner_solver->Smoother_LinSolver(u, v, mat_vec, *normal_prec, 0, smooth_iter, res, false, config); + } + }; + nested_prec = new CAbstractPreconditioner(f); + } + const auto* precond = nested ? nested_prec : normal_prec; - /*--- Solve system. ---*/ + /*--- Solve system. ---*/ - ScalarType residual = 0.0; + ScalarType residual = 0.0; - switch (KindSolver) { - case BCGSTAB: - IterLinSol = BCGSTAB_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual, - ScreenOutput, config); - break; - case FGMRES: - IterLinSol = FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual, - ScreenOutput, config); - break; - case FGCRODR: - IterLinSol = FGCRODR_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual, - ScreenOutput, config); - break; - case RESTARTED_FGMRES: - IterLinSol = RFGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual, - ScreenOutput, config); - break; - case CONJUGATE_GRADIENT: - IterLinSol = CG_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual, - ScreenOutput, config); - break; - case SMOOTHER: - IterLinSol = Smoother_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual, + switch (KindSolver) { + case BCGSTAB: + IterLinSol = BCGSTAB_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, + residual, ScreenOutput, config); + break; + case FGMRES: + IterLinSol = FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual, ScreenOutput, config); - break; - case PASTIX_LDLT: - case PASTIX_LU: - Jacobian.BuildPastixPreconditioner(geometry, config, KindSolver); - Jacobian.ComputePastixPreconditioner(*LinSysRes_ptr, *LinSysSol_ptr, geometry, config); - IterLinSol = 1; - residual = 1e-20; - break; - default: - SU2_MPI::Error("Unknown type of linear solver.", CURRENT_FUNCTION); - } + break; + case FGCRODR: + IterLinSol = FGCRODR_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, + residual, ScreenOutput, config); + break; + case RESTARTED_FGMRES: + IterLinSol = RFGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, + residual, ScreenOutput, config); + break; + case CONJUGATE_GRADIENT: + IterLinSol = CG_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, residual, + ScreenOutput, config); + break; + case SMOOTHER: + IterLinSol = Smoother_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, + residual, ScreenOutput, config); + break; + case PASTIX_LDLT: + case PASTIX_LU: + Jacobian.BuildPastixPreconditioner(geometry, config, KindSolver); + Jacobian.ComputePastixPreconditioner(*LinSysRes_ptr, *LinSysSol_ptr, geometry, config); + IterLinSol = 1; + residual = 1e-20; + break; + default: + SU2_MPI::Error("Unknown type of linear solver.", CURRENT_FUNCTION); + } - SU2_OMP_MASTER { - Residual = residual; - Iterations = IterLinSol; + SU2_OMP_MASTER { + Residual = residual; + Iterations = IterLinSol; + } + END_SU2_OMP_MASTER + + delete normal_prec; + delete nested_prec; } - END_SU2_OMP_MASTER HandleTemporariesOut(LinSysSol); - delete normal_prec; - delete nested_prec; - if (TapeActive) { /*--- To keep the behavior of SU2_DOT, but not strictly required since jacobian is symmetric(?). ---*/ const bool RequiresTranspose = diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp index 27f07ce68b2..c5c4dc9a3cb 100644 --- a/Common/src/linear_algebra/CSysVector.cpp +++ b/Common/src/linear_algebra/CSysVector.cpp @@ -65,6 +65,10 @@ void CSysVector::Initialize(unsigned long numBlk, unsigned long numB for (auto i = 0ul; i < nElm; i++) vec_val[i] = val[i]; } } + + host_data_valid = true; + device_data_valid = false; + device_context_id = 0; } template @@ -151,6 +155,7 @@ const su2matrix& CSysVector::multiDot(const std::vector< template CSysVector::~CSysVector() { + VecExpr::UnregisterDeviceModifiedVector(this); if constexpr (!std::is_trivial_v) { for (auto i = 0ul; i < nElm; i++) vec_val[i].~ScalarType(); } diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu index 7b06efbb7a0..e6429fcb351 100644 --- a/Common/src/linear_algebra/CSysVectorGPU.cu +++ b/Common/src/linear_algebra/CSysVectorGPU.cu @@ -28,7 +28,9 @@ #include "../../include/linear_algebra/CSysVector.hpp" #include "../../include/linear_algebra/GPUComms.cuh" #include +#include #include +#include namespace { @@ -55,34 +57,20 @@ void ReleaseBlasHandle(cublasHandle_t handle, bool owns_handle, const char* dest } } -} // namespace - -namespace VecExpr { - -namespace { -thread_local bool device_expressions_enabled = false; -} // namespace - -bool DeviceExpressionsEnabled() { return device_expressions_enabled; } - -void SetDeviceExpressionsEnabled(bool enabled) { device_expressions_enabled = enabled; } - -} // namespace VecExpr - -void SU2_GPU_BeginSolverBLASContext() { +void BeginDeviceBlasContext() { if (active_solver_blas_depth == 0) { cublasHandle_t handle = nullptr; if (cublasCreate(&handle) != CUBLAS_STATUS_SUCCESS) { - SU2_MPI::Error("cuBLAS handle creation failed for the GPU linear solver context.", CURRENT_FUNCTION); + SU2_MPI::Error("cuBLAS handle creation failed for the GPU linear algebra context.", CURRENT_FUNCTION); } active_solver_blas_handle = handle; } ++active_solver_blas_depth; } -void SU2_GPU_EndSolverBLASContext() { +void EndDeviceBlasContext() { if (active_solver_blas_depth == 0) { - SU2_MPI::Error("GPU linear solver BLAS context ended without a matching begin.", CURRENT_FUNCTION); + SU2_MPI::Error("GPU linear algebra BLAS context ended without a matching begin.", CURRENT_FUNCTION); } --active_solver_blas_depth; @@ -90,11 +78,81 @@ void SU2_GPU_EndSolverBLASContext() { auto status = cublasDestroy(active_solver_blas_handle); active_solver_blas_handle = nullptr; if (status != CUBLAS_STATUS_SUCCESS) { - SU2_MPI::Error("cuBLAS handle destruction failed for the GPU linear solver context.", CURRENT_FUNCTION); + SU2_MPI::Error("cuBLAS handle destruction failed for the GPU linear algebra context.", CURRENT_FUNCTION); } } } +} // namespace + +namespace VecExpr { + +namespace { +thread_local bool device_expressions_enabled = false; +struct DeviceExpressionContextState { + bool enabled = false; + bool previous = false; + unsigned long previous_context_id = 0; +}; +struct ModifiedVector { + const void* vector = nullptr; + void (*sync)(const void*) = nullptr; +}; +thread_local std::vector device_context_stack; +thread_local std::vector modified_vectors; +thread_local unsigned long next_device_context_id = 0; +thread_local unsigned long current_device_context_id = 0; +} // namespace + +bool DeviceExpressionsEnabled() { return device_expressions_enabled; } + +void BeginDeviceExpressionContext(bool enabled) { + device_context_stack.push_back({enabled, device_expressions_enabled, current_device_context_id}); + if (!enabled) return; + + BeginDeviceBlasContext(); + device_expressions_enabled = true; + current_device_context_id = ++next_device_context_id; +} + +void FlushDeviceExpressionContext() { + for (const auto& entry : modified_vectors) entry.sync(entry.vector); + modified_vectors.clear(); +} + +void EndDeviceExpressionContext() { + if (device_context_stack.empty()) { + SU2_MPI::Error("Device expression context ended without a matching begin.", CURRENT_FUNCTION); + } + + const auto state = device_context_stack.back(); + device_context_stack.pop_back(); + + if (state.enabled) { + FlushDeviceExpressionContext(); + device_expressions_enabled = state.previous; + current_device_context_id = state.previous_context_id; + EndDeviceBlasContext(); + } +} + +unsigned long CurrentDeviceExpressionContextId() { return current_device_context_id; } + +void RegisterDeviceModifiedVector(const void* vector, void (*sync)(const void*)) { + if (!DeviceExpressionsEnabled()) return; + const auto exists = std::any_of(modified_vectors.begin(), modified_vectors.end(), + [vector](const auto& entry) { return entry.vector == vector; }); + if (!exists) modified_vectors.push_back({vector, sync}); +} + +void UnregisterDeviceModifiedVector(const void* vector) { + modified_vectors.erase(std::remove_if(modified_vectors.begin(), modified_vectors.end(), + [vector](const auto& entry) { return entry.vector == vector; }), + modified_vectors.end()); +} + +} // namespace VecExpr + template void CSysVector::HtDTransfer(bool trigger) const { if (trigger) @@ -109,6 +167,9 @@ void CSysVector::DtHTransfer(bool trigger) const { template ScalarType CSysVector::GPUDot(const CSysVector& other) const { + EnsureDeviceData(); + other.EnsureDeviceData(); + bool owns_handle = false; cublasHandle_t handle = GetBlasHandle("cuBLAS handle creation failed in CSysVector::GPUDot.", owns_handle); cublasStatus_t status = CUBLAS_STATUS_SUCCESS; From 5e16f4d9d8895c78a1e2a20fcb1cb5e6a3ba7a4b Mon Sep 17 00:00:00 2001 From: LwhJesse <256257451+LwhJesse@users.noreply.github.com> Date: Sun, 28 Jun 2026 22:11:06 +0800 Subject: [PATCH 10/10] Fix CI ThreadSanitizer data race ThreadSanitizer hybrid_regression jobs reported a data race in CSysVector::MarkHostDataModified() on host/device validity flags. Avoid updating those CUDA-only validity flags in non-CUDA builds. --- Common/include/linear_algebra/CSysVector.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp index e8fcd624cbd..aabb1aa69a4 100644 --- a/Common/include/linear_algebra/CSysVector.hpp +++ b/Common/include/linear_algebra/CSysVector.hpp @@ -271,8 +271,10 @@ class CSysVector : public VecExpr::CVecExpr, ScalarType> } void MarkHostDataModified() const { +#ifdef HAVE_CUDA host_data_valid = true; device_data_valid = false; +#endif } public: