From 9d098128778d13efef85ffe4efa99401999032fe Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Thu, 29 May 2025 16:34:17 +0200 Subject: [PATCH] formatting --- opm/simulators/linalg/gpuistl/GpuDILU.cpp | 183 +++++++++-------- opm/simulators/linalg/gpuistl/GpuDILU.hpp | 29 ++- .../preconditionerKernels/DILUKernels.cu | 128 ++++++------ .../preconditionerKernels/DILUKernels.hpp | 22 +- .../ILU_variants_helper_kernels.cu | 189 +++++++++--------- .../ILU_variants_helper_kernels.hpp | 52 ++--- 6 files changed, 306 insertions(+), 297 deletions(-) diff --git a/opm/simulators/linalg/gpuistl/GpuDILU.cpp b/opm/simulators/linalg/gpuistl/GpuDILU.cpp index 07ffdfa4e7..3cc141787e 100644 --- a/opm/simulators/linalg/gpuistl/GpuDILU.cpp +++ b/opm/simulators/linalg/gpuistl/GpuDILU.cpp @@ -21,28 +21,33 @@ #include #include #include +#include #include #include #include #include -#include #include #include #include +#include #include #include #include #include -#include -#include -#include #include +#include +#include namespace Opm::gpuistl { template -GpuDILU::GpuDILU(const typename GpuDILU::GPUMatrix& gpuMatrix, const M& cpuMatrix, bool splitMatrix, bool tuneKernels, int mixedPrecisionScheme, bool reorder) +GpuDILU::GpuDILU(const typename GpuDILU::GPUMatrix& gpuMatrix, + const M& cpuMatrix, + bool splitMatrix, + bool tuneKernels, + int mixedPrecisionScheme, + bool reorder) : m_levelSets(Opm::getMatrixRowColoring(cpuMatrix, Opm::ColoringType::LOWER)) , m_reorderedToNatural(detail::createReorderedToNatural(m_levelSets)) , m_naturalToReordered(detail::createNaturalToReordered(m_levelSets)) @@ -59,8 +64,9 @@ GpuDILU::GpuDILU(const typename GpuDILU::GPUMatrix& gpuM { // TODO: Should in some way verify that this matrix is symmetric, only do it debug mode? // Some sanity check - OPM_ERROR_IF(cpuMatrix.N() != m_gpuMatrix.N(), - fmt::format("CuSparse matrix not same size as DUNE matrix. {} vs {}.", m_gpuMatrix.N(), cpuMatrix.N())); + OPM_ERROR_IF( + cpuMatrix.N() != m_gpuMatrix.N(), + fmt::format("CuSparse matrix not same size as DUNE matrix. {} vs {}.", m_gpuMatrix.N(), cpuMatrix.N())); OPM_ERROR_IF(cpuMatrix[0][0].N() != m_gpuMatrix.blockSize(), fmt::format("CuSparse matrix not same blocksize as DUNE matrix. {} vs {}.", m_gpuMatrix.blockSize(), @@ -74,57 +80,56 @@ GpuDILU::GpuDILU(const typename GpuDILU::GPUMatrix& gpuM m_gpuMatrix.nonzeroes(), cpuMatrix.nonzeroes())); if (m_reorder && (m_splitMatrix || m_mixedPrecisionScheme != MatrixStorageMPScheme::DOUBLE_DIAG_DOUBLE_OFFDIAG)) { - OPM_THROW(std::runtime_error, "Reordering is only supported for full precision matrices without matrix splitting."); + OPM_THROW(std::runtime_error, + "Reordering is only supported for full precision matrices without matrix splitting."); } if (m_splitMatrix) { m_gpuMatrixReorderedDiag = std::make_unique>(blocksize_ * blocksize_ * cpuMatrix.N()); std::tie(m_gpuMatrixReorderedLower, m_gpuMatrixReorderedUpper) = detail::extractLowerAndUpperMatrices>(cpuMatrix, - m_reorderedToNatural); - } - else { + m_reorderedToNatural); + } else { m_gpuMatrixReordered = detail::createReorderedMatrix>( cpuMatrix, m_reorderedToNatural); } if (m_mixedPrecisionScheme != MatrixStorageMPScheme::DOUBLE_DIAG_DOUBLE_OFFDIAG) { - if (!m_splitMatrix){ + if (!m_splitMatrix) { OPM_THROW(std::runtime_error, "Matrix must be split when storing as float."); } - m_gpuMatrixReorderedLowerFloat = std::make_unique(m_gpuMatrixReorderedLower->getRowIndices(), m_gpuMatrixReorderedLower->getColumnIndices(), blocksize_); - m_gpuMatrixReorderedUpperFloat = std::make_unique(m_gpuMatrixReorderedUpper->getRowIndices(), m_gpuMatrixReorderedUpper->getColumnIndices(), blocksize_); + m_gpuMatrixReorderedLowerFloat = std::make_unique( + m_gpuMatrixReorderedLower->getRowIndices(), m_gpuMatrixReorderedLower->getColumnIndices(), blocksize_); + m_gpuMatrixReorderedUpperFloat = std::make_unique( + m_gpuMatrixReorderedUpper->getRowIndices(), m_gpuMatrixReorderedUpper->getColumnIndices(), blocksize_); if (m_mixedPrecisionScheme == MatrixStorageMPScheme::FLOAT_DIAG_FLOAT_OFFDIAG) { - m_gpuDInvFloat = std::make_unique(m_gpuMatrix.N() * m_gpuMatrix.blockSize() * m_gpuMatrix.blockSize()); + m_gpuDInvFloat + = std::make_unique(m_gpuMatrix.N() * m_gpuMatrix.blockSize() * m_gpuMatrix.blockSize()); } } // Precompute diagaon indices and handle the actual reordering if that is enabled. - if (m_reorder){ - detail::computeDiagIndices( - m_gpuMatrixReordered->getNonZeroValues().data(), - m_gpuMatrixReordered->getRowIndices().data(), - m_gpuMatrixReordered->getColumnIndices().data(), - m_gpuReorderToNatural.data(), - m_gpuMatrix.N(), - m_diagIdxs.data()); + if (m_reorder) { + detail::computeDiagIndices(m_gpuMatrixReordered->getNonZeroValues().data(), + m_gpuMatrixReordered->getRowIndices().data(), + m_gpuMatrixReordered->getColumnIndices().data(), + m_gpuReorderToNatural.data(), + m_gpuMatrix.N(), + m_diagIdxs.data()); reorderAndSplitMatrix(m_moveThreadBlockSize); - } - else{ - detail::computeDiagIndicesNoReorder( - m_gpuMatrix.getNonZeroValues().data(), - m_gpuMatrix.getRowIndices().data(), - m_gpuMatrix.getColumnIndices().data(), - m_gpuLevelSets.data(), - m_gpuMatrix.N(), - m_diagIdxs.data()); + } else { + detail::computeDiagIndicesNoReorder(m_gpuMatrix.getNonZeroValues().data(), + m_gpuMatrix.getRowIndices().data(), + m_gpuMatrix.getColumnIndices().data(), + m_gpuLevelSets.data(), + m_gpuMatrix.N(), + m_diagIdxs.data()); } computeDiagonal(m_DILUFactorizationThreadBlockSize); if (m_tuneThreadBlockSizes) { tuneThreadBlockSizes(); } - } template @@ -154,7 +159,8 @@ GpuDILU::apply(X& v, const Y& d) apply(v, d, m_lowerSolveThreadBlockSize, m_upperSolveThreadBlockSize); OPM_GPU_SAFE_CALL(cudaStreamEndCapture(m_stream.get(), &m_apply_graphs[ptrs].get())); - OPM_GPU_SAFE_CALL(cudaGraphInstantiate(&m_executableGraphs[ptrs].get(), m_apply_graphs[ptrs].get(), nullptr, nullptr, 0)); + OPM_GPU_SAFE_CALL( + cudaGraphInstantiate(&m_executableGraphs[ptrs].get(), m_apply_graphs[ptrs].get(), nullptr, nullptr, 0)); } OPM_GPU_SAFE_CALL(cudaGraphLaunch(m_executableGraphs[ptrs].get(), 0)); } @@ -227,18 +233,17 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int lowerSolveThreadBlockSize, m_stream.get()); } else { - detail::DILU::solveLowerLevelSetNoReorder( - m_gpuMatrix.getNonZeroValues().data(), - m_gpuMatrix.getRowIndices().data(), - m_gpuMatrix.getColumnIndices().data(), - m_gpuLevelSets.data(), - levelStartIdx, - numOfRowsInLevel, - m_gpuDInv.data(), - d.data(), - v.data(), - lowerSolveThreadBlockSize, - m_stream.get()); + detail::DILU::solveLowerLevelSetNoReorder(m_gpuMatrix.getNonZeroValues().data(), + m_gpuMatrix.getRowIndices().data(), + m_gpuMatrix.getColumnIndices().data(), + m_gpuLevelSets.data(), + levelStartIdx, + numOfRowsInLevel, + m_gpuDInv.data(), + d.data(), + v.data(), + lowerSolveThreadBlockSize, + m_stream.get()); } } levelStartIdx += numOfRowsInLevel; @@ -250,7 +255,7 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int const int numOfRowsInLevel = m_levelSets[level].size(); levelStartIdx -= numOfRowsInLevel; if (m_splitMatrix) { - if (m_mixedPrecisionScheme == MatrixStorageMPScheme::FLOAT_DIAG_FLOAT_OFFDIAG){ + if (m_mixedPrecisionScheme == MatrixStorageMPScheme::FLOAT_DIAG_FLOAT_OFFDIAG) { detail::DILU::solveUpperLevelSetSplit( m_gpuMatrixReorderedUpperFloat->getNonZeroValues().data(), m_gpuMatrixReorderedUpperFloat->getRowIndices().data(), @@ -262,7 +267,7 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int v.data(), upperSolveThreadBlockSize, m_stream.get()); - } else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG){ + } else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG) { detail::DILU::solveUpperLevelSetSplit( m_gpuMatrixReorderedUpperFloat->getNonZeroValues().data(), m_gpuMatrixReorderedUpperFloat->getRowIndices().data(), @@ -301,17 +306,16 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int upperSolveThreadBlockSize, m_stream.get()); } else { - detail::DILU::solveUpperLevelSetNoReorder( - m_gpuMatrix.getNonZeroValues().data(), - m_gpuMatrix.getRowIndices().data(), - m_gpuMatrix.getColumnIndices().data(), - m_gpuLevelSets.data(), - levelStartIdx, - numOfRowsInLevel, - m_gpuDInv.data(), - v.data(), - upperSolveThreadBlockSize, - m_stream.get()); + detail::DILU::solveUpperLevelSetNoReorder(m_gpuMatrix.getNonZeroValues().data(), + m_gpuMatrix.getRowIndices().data(), + m_gpuMatrix.getColumnIndices().data(), + m_gpuLevelSets.data(), + levelStartIdx, + numOfRowsInLevel, + m_gpuDInv.data(), + v.data(), + upperSolveThreadBlockSize, + m_stream.get()); } } } @@ -381,12 +385,12 @@ GpuDILU::reorderAndSplitMatrix(int moveThreadBlockSize) moveThreadBlockSize); } else { detail::copyMatDataToReordered(m_gpuMatrix.getNonZeroValues().data(), - m_gpuMatrix.getRowIndices().data(), - m_gpuMatrixReordered->getNonZeroValues().data(), - m_gpuMatrixReordered->getRowIndices().data(), - m_gpuNaturalToReorder.data(), - m_gpuMatrixReordered->N(), - moveThreadBlockSize); + m_gpuMatrix.getRowIndices().data(), + m_gpuMatrixReordered->getNonZeroValues().data(), + m_gpuMatrixReordered->getRowIndices().data(), + m_gpuNaturalToReorder.data(), + m_gpuMatrixReordered->N(), + moveThreadBlockSize); } } @@ -399,7 +403,10 @@ GpuDILU::computeDiagonal(int factorizationBlockSize) const int numOfRowsInLevel = m_levelSets[level].size(); if (m_splitMatrix) { if (m_mixedPrecisionScheme == MatrixStorageMPScheme::FLOAT_DIAG_FLOAT_OFFDIAG) { - detail::DILU::computeDiluDiagonalSplit( + detail::DILU::computeDiluDiagonalSplit( m_gpuMatrixReorderedLower->getNonZeroValues().data(), m_gpuMatrixReorderedLower->getRowIndices().data(), m_gpuMatrixReorderedLower->getColumnIndices().data(), @@ -417,7 +424,10 @@ GpuDILU::computeDiagonal(int factorizationBlockSize) m_gpuMatrixReorderedUpperFloat->getNonZeroValues().data(), factorizationBlockSize); } else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG) { - detail::DILU::computeDiluDiagonalSplit( + detail::DILU::computeDiluDiagonalSplit( m_gpuMatrixReorderedLower->getNonZeroValues().data(), m_gpuMatrixReorderedLower->getRowIndices().data(), m_gpuMatrixReorderedLower->getColumnIndices().data(), @@ -436,7 +446,10 @@ GpuDILU::computeDiagonal(int factorizationBlockSize) factorizationBlockSize); } else { // TODO: should this be field type twice or field type then float in the template? - detail::DILU::computeDiluDiagonalSplit( + detail::DILU::computeDiluDiagonalSplit( m_gpuMatrixReorderedLower->getNonZeroValues().data(), m_gpuMatrixReorderedLower->getRowIndices().data(), m_gpuMatrixReorderedLower->getColumnIndices().data(), @@ -490,41 +503,45 @@ GpuDILU::tuneThreadBlockSizes() { if (m_reorder) { // tune the thread-block size of the update function - auto tuneMoveThreadBlockSizeInUpdate = [this](int moveThreadBlockSize){ + auto tuneMoveThreadBlockSizeInUpdate = [this](int moveThreadBlockSize) { this->update(moveThreadBlockSize, m_DILUFactorizationThreadBlockSize); }; - m_moveThreadBlockSize = detail::tuneThreadBlockSize(tuneMoveThreadBlockSizeInUpdate, "(in DILU update) Move data to reordered matrix"); + m_moveThreadBlockSize = detail::tuneThreadBlockSize(tuneMoveThreadBlockSizeInUpdate, + "(in DILU update) Move data to reordered matrix"); } - auto tuneFactorizationThreadBlockSizeInUpdate = [this](int factorizationThreadBlockSize){ + auto tuneFactorizationThreadBlockSizeInUpdate = [this](int factorizationThreadBlockSize) { this->update(m_moveThreadBlockSize, factorizationThreadBlockSize); }; - m_DILUFactorizationThreadBlockSize = detail::tuneThreadBlockSize(tuneFactorizationThreadBlockSizeInUpdate, "(in DILU update) DILU factorization"); + m_DILUFactorizationThreadBlockSize + = detail::tuneThreadBlockSize(tuneFactorizationThreadBlockSizeInUpdate, "(in DILU update) DILU factorization"); // tune the thread-block size of the apply GpuVector tmpV(m_gpuMatrix.N() * m_gpuMatrix.blockSize()); GpuVector tmpD(m_gpuMatrix.N() * m_gpuMatrix.blockSize()); tmpD = 1; - auto tuneLowerSolveThreadBlockSizeInApply = [this, &tmpV, &tmpD](int lowerSolveThreadBlockSize){ + auto tuneLowerSolveThreadBlockSizeInApply = [this, &tmpV, &tmpD](int lowerSolveThreadBlockSize) { this->apply(tmpV, tmpD, lowerSolveThreadBlockSize, m_DILUFactorizationThreadBlockSize); }; - m_lowerSolveThreadBlockSize = detail::tuneThreadBlockSize(tuneLowerSolveThreadBlockSizeInApply, "(in DILU apply) Triangular lower solve"); + m_lowerSolveThreadBlockSize + = detail::tuneThreadBlockSize(tuneLowerSolveThreadBlockSizeInApply, "(in DILU apply) Triangular lower solve"); - auto tuneUpperSolveThreadBlockSizeInApply = [this, &tmpV, &tmpD](int upperSolveThreadBlockSize){ + auto tuneUpperSolveThreadBlockSizeInApply = [this, &tmpV, &tmpD](int upperSolveThreadBlockSize) { this->apply(tmpV, tmpD, m_lowerSolveThreadBlockSize, upperSolveThreadBlockSize); }; - m_upperSolveThreadBlockSize = detail::tuneThreadBlockSize(tuneUpperSolveThreadBlockSizeInApply, "(in DILU apply) Triangular upper solve"); + m_upperSolveThreadBlockSize + = detail::tuneThreadBlockSize(tuneUpperSolveThreadBlockSizeInApply, "(in DILU apply) Triangular upper solve"); } } // namespace Opm::gpuistl -#define INSTANTIATE_CUDILU_DUNE(realtype, blockdim) \ - template class ::Opm::gpuistl::GpuDILU>, \ - ::Opm::gpuistl::GpuVector, \ - ::Opm::gpuistl::GpuVector>; \ - template class ::Opm::gpuistl::GpuDILU>, \ - ::Opm::gpuistl::GpuVector, \ - ::Opm::gpuistl::GpuVector> +#define INSTANTIATE_CUDILU_DUNE(realtype, blockdim) \ + template class ::Opm::gpuistl::GpuDILU>, \ + ::Opm::gpuistl::GpuVector, \ + ::Opm::gpuistl::GpuVector>; \ + template class ::Opm::gpuistl::GpuDILU>, \ + ::Opm::gpuistl::GpuVector, \ + ::Opm::gpuistl::GpuVector> INSTANTIATE_CUDILU_DUNE(double, 1); INSTANTIATE_CUDILU_DUNE(double, 2); diff --git a/opm/simulators/linalg/gpuistl/GpuDILU.hpp b/opm/simulators/linalg/gpuistl/GpuDILU.hpp index e69addcea2..5dce74606b 100644 --- a/opm/simulators/linalg/gpuistl/GpuDILU.hpp +++ b/opm/simulators/linalg/gpuistl/GpuDILU.hpp @@ -19,18 +19,18 @@ #ifndef OPM_GPUDILU_HPP #define OPM_GPUDILU_HPP +#include #include #include #include -#include -#include #include +#include #include #include +#include #include -#include -#include #include +#include namespace Opm::gpuistl @@ -66,11 +66,16 @@ public: //! \brief The matrix type the preconditioner is for. using matrix_type = GPUMatrix; - + //! \brief Constructor. //! //! Constructor gets all parameters to operate the prec. - explicit GpuDILU(const GPUMatrix& gpuMatrix, const CPUMatrixT& cpuMatrix, bool splitMatrix, bool tuneKernels, int mixedPrecisionScheme, bool reorder); + explicit GpuDILU(const GPUMatrix& gpuMatrix, + const CPUMatrixT& cpuMatrix, + bool splitMatrix, + bool tuneKernels, + int mixedPrecisionScheme, + bool reorder); //! \brief Prepare the preconditioner. //! \note Does nothing at the time being. @@ -111,7 +116,8 @@ public: return false; } - virtual bool hasPerfectUpdate() const override { + virtual bool hasPerfectUpdate() const override + { return true; } @@ -148,7 +154,7 @@ private: //! row conversion from reordered to natural matrix indices stored on the GPU GpuVector m_gpuReorderToNatural; //! \brief Potentially stores the level sets on the GPU - GpuBuffer m_gpuLevelSets; //TODO: make this an optional to not spend time on the copy if not needed + GpuBuffer m_gpuLevelSets; // TODO: make this an optional to not spend time on the copy if not needed //! \brief Stores the inverted diagonal that we use in DILU GpuVector m_gpuDInv; //! \brief Bool storing whether or not we should store matrices in a split format @@ -166,6 +172,7 @@ private: //! \brief The reorder parameter can be set to false to shorten the update time, but slows the apply time bool m_reorder; + //! \brief The indices of the blocks that are the diagonal entries GpuBuffer m_diagIdxs; // Graphs for Apply @@ -173,10 +180,10 @@ private: std::map, GPUGraphExec> m_executableGraphs; // Stream for the DILU operations on the GPU - GPUStream m_stream{}; + GPUStream m_stream {}; // Events for synchronization with main stream - GPUEvent m_before{}; - GPUEvent m_after{}; + GPUEvent m_before {}; + GPUEvent m_after {}; }; } // end namespace Opm::gpuistl diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu index 48b7932b32..4a7780a783 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu @@ -61,14 +61,14 @@ namespace template __global__ void cuSolveLowerLevelSetNoReorder(const T* mat, - const int* rowIndices, - const int* colIndices, - const size_t* indexConversion, - int startIdx, - int rowsInLevelSet, - const T* dInv, - const T* d, - T* v) + const int* rowIndices, + const int* colIndices, + const size_t* indexConversion, + int startIdx, + int rowsInLevelSet, + const T* dInv, + const T* d, + T* v) { const auto rawIdx = startIdx + (blockDim.x * blockIdx.x + threadIdx.x); if (rawIdx < rowsInLevelSet + startIdx) { @@ -153,13 +153,13 @@ namespace template __global__ void cuSolveUpperLevelSetNoReorder(const T* mat, - const int* rowIndices, - const int* colIndices, - const size_t* indexConversion, - int startIdx, - int rowsInLevelSet, - const T* dInv, - T* v) + const int* rowIndices, + const int* colIndices, + const size_t* indexConversion, + int startIdx, + int rowsInLevelSet, + const T* dInv, + T* v) { const auto rawIdx = startIdx + (blockDim.x * blockIdx.x + threadIdx.x); if (rawIdx < rowsInLevelSet + startIdx) { @@ -270,13 +270,13 @@ namespace template __global__ void cuComputeDiluDiagonalNoReorder(const T* mat, - const int* rowIndices, - const int* colIndices, - const size_t* indexConversion, - const size_t* diagIndices, - const int startIdx, - int rowsInLevelSet, - T* dInv) + const int* rowIndices, + const int* colIndices, + const size_t* indexConversion, + const size_t* diagIndices, + const int startIdx, + int rowsInLevelSet, + T* dInv) { const auto rawIdx = startIdx + blockDim.x * blockIdx.x + threadIdx.x; if (rawIdx < rowsInLevelSet + startIdx) { @@ -428,19 +428,19 @@ solveLowerLevelSet(T* reorderedMat, template void solveLowerLevelSetNoReorder(const T* reorderedMat, - const int* rowIndices, - const int* colIndices, - const size_t* indexConversion, - int startIdx, - int rowsInLevelSet, - const T* dInv, - const T* d, - T* v, - int thrBlockSize, - cudaStream_t stream) + const int* rowIndices, + const int* colIndices, + const size_t* indexConversion, + int startIdx, + int rowsInLevelSet, + const T* dInv, + const T* d, + T* v, + int thrBlockSize, + cudaStream_t stream) { - int threadBlockSize - = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveLowerLevelSetNoReorder, thrBlockSize); + int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( + cuSolveLowerLevelSetNoReorder, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); cuSolveLowerLevelSetNoReorder<<>>( reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, d, v); @@ -491,18 +491,18 @@ solveUpperLevelSet(T* reorderedMat, template void solveUpperLevelSetNoReorder(const T* reorderedMat, - const int* rowIndices, - const int* colIndices, - const size_t* indexConversion, - int startIdx, - int rowsInLevelSet, - const T* dInv, - T* v, - int thrBlockSize, - cudaStream_t stream) + const int* rowIndices, + const int* colIndices, + const size_t* indexConversion, + int startIdx, + int rowsInLevelSet, + const T* dInv, + T* v, + int thrBlockSize, + cudaStream_t stream) { - int threadBlockSize - = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveUpperLevelSetNoReorder, thrBlockSize); + int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( + cuSolveUpperLevelSetNoReorder, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); cuSolveUpperLevelSetNoReorder<<>>( reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v); @@ -563,27 +563,21 @@ computeDiluDiagonal(T* reorderedMat, template void computeDiluDiagonalNoReorder(const T* reorderedMat, - const int* rowIndices, - const int* colIndices, - const size_t* indexConversion, - const size_t* diagIndices, - int startIdx, - int rowsInLevelSet, - T* dInv, - int thrBlockSize) + const int* rowIndices, + const int* colIndices, + const size_t* indexConversion, + const size_t* diagIndices, + int startIdx, + int rowsInLevelSet, + T* dInv, + int thrBlockSize) { if (blocksize <= 3) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( cuComputeDiluDiagonalNoReorder, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuComputeDiluDiagonalNoReorder<<>>(reorderedMat, - rowIndices, - colIndices, - indexConversion, - diagIndices, - startIdx, - rowsInLevelSet, - dInv); + cuComputeDiluDiagonalNoReorder<<>>( + reorderedMat, rowIndices, colIndices, indexConversion, diagIndices, startIdx, rowsInLevelSet, dInv); } else { OPM_THROW(std::invalid_argument, "Inverting diagonal is not implemented for blocksizes > 3"); } @@ -634,8 +628,9 @@ computeDiluDiagonalSplit(const InputScalar* srcReorderedLowerMat, } #define INSTANTIATE_KERNEL_WRAPPERS(T, blocksize) \ - template void computeDiluDiagonal(T*, int*, int*, int*, int*, size_t*, const int, int, T*, int); \ - template void computeDiluDiagonalNoReorder(const T*, const int*, const int*, const size_t*, const size_t*, int, int, T*, int); \ + template void computeDiluDiagonal(T*, int*, int*, int*, int*, size_t*, const int, int, T*, int); \ + template void computeDiluDiagonalNoReorder( \ + const T*, const int*, const int*, const size_t*, const size_t*, int, int, T*, int); \ template void computeDiluDiagonalSplit( \ const T*, \ int*, \ @@ -739,10 +734,11 @@ computeDiluDiagonalSplit(const InputScalar* srcReorderedLowerMat, double*, \ int); \ template void solveUpperLevelSet(T*, int*, int*, int*, int, int, const T*, T*, int, cudaStream_t); \ - template void solveUpperLevelSetNoReorder(const T*, const int*, const int*, const size_t*, int, int, const T*, T*, int, cudaStream_t); \ + template void solveUpperLevelSetNoReorder( \ + const T*, const int*, const int*, const size_t*, int, int, const T*, T*, int, cudaStream_t); \ template void solveLowerLevelSet( \ - T*, int*, int*, int*, int, int, const T*, const T*, T*, int, cudaStream_t); \ - template void solveLowerLevelSetNoReorder( \ + T*, int*, int*, int*, int, int, const T*, const T*, T*, int, cudaStream_t); \ + template void solveLowerLevelSetNoReorder( \ const T*, const int*, const int*, const size_t*, int, int, const T*, const T*, T*, int, cudaStream_t); INSTANTIATE_KERNEL_WRAPPERS(float, 1); diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp index 3752ad49d5..42e2d0c81d 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp @@ -59,16 +59,16 @@ void solveLowerLevelSet(T* reorderedMat, template void solveLowerLevelSetNoReorder(const T* reorderedMat, - const int* rowIndices, - const int* colIndices, - const size_t* indexConversion, - int startIdx, - int rowsInLevelSet, - const T* dInv, - const T* d, - T* v, - int threadBlockSize, - cudaStream_t stream); + const int* rowIndices, + const int* colIndices, + const size_t* indexConversion, + int startIdx, + int rowsInLevelSet, + const T* dInv, + const T* d, + T* v, + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform a lower solve on certain rows in a matrix that can safely be computed in parallel @@ -193,7 +193,7 @@ void computeDiluDiagonal(T* reorderedMat, T* dInv, int threadBlockSize); - /** +/** * @brief Computes the DILU of a BCSR matrix and stores it in a vector containing the diagonal blocks * @param mat pointer to GPU memory containing nonzerovalues of the sparse matrix. The matrix reordered such * that rows in the same level sets are contiguous diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.cu b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.cu index 11be1d72bd..dc53ddd471 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.cu +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.cu @@ -21,126 +21,115 @@ namespace Opm::gpuistl::detail { - namespace { - // precompute the diagonal indices to speedup the update - // The matrix here is not reordered, so the i'th thread computes - // the diagonal element index not the i'th row of the matrix, - // That is what we would have done if there was no coloring. - // Instead it represents what the j'th rows in the k'th levelset - // correspons to via the indexConversion array. - template - __global__ void cuComputeDiagIndicesNoReorder(const T* mat, - const int* rowIndices, - const int* colIndices, - const size_t* indexConversion, - int rows, - size_t* diagIndices) - { - const auto rawIdx = blockDim.x * blockIdx.x + threadIdx.x; - if (rawIdx < rows) { - const size_t realRowIdx = indexConversion[rawIdx]; - const size_t nnzIdx = rowIndices[realRowIdx]; - - size_t diagIdx = nnzIdx; - while (colIndices[diagIdx] != realRowIdx) { - ++diagIdx; - } - - diagIndices[realRowIdx] = diagIdx; - } - } - - // In this case we have a reordered matrix, so the firsth i rows will - // in the same levelset and contiguous in the matrix. For this reason - // we convert from the reordered index (i) to the natural index (j) which - // will give us the diagonal element index as the row and column indices match there - template - __global__ void cuComputeDiagIndices(const T* mat, - const int* rowIndices, - const int* colIndices, - const int* reorderedToNatural, - int rows, - size_t* diagIndices) - { - const auto rawIdx = blockDim.x * blockIdx.x + threadIdx.x; - if (rawIdx < rows) { - const int naturalRowIdx = reorderedToNatural[rawIdx]; - const size_t nnzIdx = rowIndices[rawIdx]; - - size_t diagIdx = nnzIdx; - while (colIndices[diagIdx] != naturalRowIdx) { - ++diagIdx; - } - - // I choose to have this indexed in the same order as the matrix in hopes - // of memory coalescing. - diagIndices[rawIdx] = diagIdx; - } - } - - } // empty namespace - - +namespace +{ + // precompute the diagonal indices to speedup the update + // The matrix here is not reordered, so the i'th thread computes + // the diagonal element index not the i'th row of the matrix, + // That is what we would have done if there was no coloring. + // Instead it represents what the j'th rows in the k'th levelset + // correspons to via the indexConversion array. template - void computeDiagIndicesNoReorder(const T* mat, + __global__ void cuComputeDiagIndicesNoReorder(const T* mat, + const int* rowIndices, + const int* colIndices, + const size_t* indexConversion, + int rows, + size_t* diagIndices) + { + const auto rawIdx = blockDim.x * blockIdx.x + threadIdx.x; + if (rawIdx < rows) { + const size_t realRowIdx = indexConversion[rawIdx]; + const size_t nnzIdx = rowIndices[realRowIdx]; + + size_t diagIdx = nnzIdx; + while (colIndices[diagIdx] != realRowIdx) { + ++diagIdx; + } + + diagIndices[realRowIdx] = diagIdx; + } + } + + // In this case we have a reordered matrix, so the firsth i rows will + // in the same levelset and contiguous in the matrix. For this reason + // we convert from the reordered index (i) to the natural index (j) which + // will give us the diagonal element index as the row and column indices match there + template + __global__ void cuComputeDiagIndices(const T* mat, const int* rowIndices, const int* colIndices, - const size_t* indexConversion, + const int* reorderedToNatural, int rows, size_t* diagIndices) { - // Have cuda/hip automatically pick a reasonable block size for this kernel based on static analysis - int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( - cuComputeDiagIndicesNoReorder - ); + const auto rawIdx = blockDim.x * blockIdx.x + threadIdx.x; + if (rawIdx < rows) { + const int naturalRowIdx = reorderedToNatural[rawIdx]; + const size_t nnzIdx = rowIndices[rawIdx]; - // Calculate the number of blocks needed - int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rows, threadBlockSize); + size_t diagIdx = nnzIdx; + while (colIndices[diagIdx] != naturalRowIdx) { + ++diagIdx; + } - // launch kernel - cuComputeDiagIndicesNoReorder<<>>( - mat, rowIndices, colIndices, indexConversion, rows, diagIndices - ); + // I choose to have this indexed in the same order as the matrix in hopes + // of memory coalescing. + diagIndices[rawIdx] = diagIdx; + } } - template - void computeDiagIndices(const T* mat, +} // namespace + + +template +void +computeDiagIndicesNoReorder(const T* mat, const int* rowIndices, const int* colIndices, - const int* reorderedToNatural, + const size_t* indexConversion, int rows, size_t* diagIndices) - { - // Have cuda/hip automatically pick a reasonable block size for this kernel based on static analysis - int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( - cuComputeDiagIndices - ); +{ + // Have cuda/hip automatically pick a reasonable block size for this kernel based on static analysis + int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuComputeDiagIndicesNoReorder); - // Calculate the number of blocks needed - int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rows, threadBlockSize); + // Calculate the number of blocks needed + int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rows, threadBlockSize); - // launch kernel - cuComputeDiagIndices<<>>( - mat, rowIndices, colIndices, reorderedToNatural, rows, diagIndices - ); - } + // launch kernel + cuComputeDiagIndicesNoReorder + <<>>(mat, rowIndices, colIndices, indexConversion, rows, diagIndices); +} + +template +void +computeDiagIndices(const T* mat, + const int* rowIndices, + const int* colIndices, + const int* reorderedToNatural, + int rows, + size_t* diagIndices) +{ + // Have cuda/hip automatically pick a reasonable block size for this kernel based on static analysis + int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuComputeDiagIndices); + + // Calculate the number of blocks needed + int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rows, threadBlockSize); + + // launch kernel + cuComputeDiagIndices + <<>>(mat, rowIndices, colIndices, reorderedToNatural, rows, diagIndices); +} } // namespace Opm::gpuistl::detail -#define DEFINE_KERNEL_WRAPPERS_FOR_TYPES(T) \ - template void Opm::gpuistl::detail::computeDiagIndicesNoReorder(const T*, \ - const int*, \ - const int*, \ - const size_t*, \ - int, \ - size_t*); \ - template void Opm::gpuistl::detail::computeDiagIndices(const T*, \ - const int*, \ - const int*, \ - const int*, \ - int, \ - size_t*); +#define DEFINE_KERNEL_WRAPPERS_FOR_TYPES(T) \ + template void Opm::gpuistl::detail::computeDiagIndicesNoReorder( \ + const T*, const int*, const int*, const size_t*, int, size_t*); \ + template void Opm::gpuistl::detail::computeDiagIndices( \ + const T*, const int*, const int*, const int*, int, size_t*); DEFINE_KERNEL_WRAPPERS_FOR_TYPES(float) DEFINE_KERNEL_WRAPPERS_FOR_TYPES(double) diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.hpp b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.hpp index ddeaf43498..0a3aa417a7 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.hpp +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.hpp @@ -43,32 +43,32 @@ namespace Opm::gpuistl::detail * @param[in] rows Number of rows in the matrix * @param[out] diagIndices Array to store the indices of diagonal elements */ - template - void computeDiagIndicesNoReorder(const T* mat, - const int* rowIndices, - const int* colIndices, - const size_t* indexConversion, - int rows, - size_t* diagIndices); +template +void computeDiagIndicesNoReorder(const T* mat, + const int* rowIndices, + const int* colIndices, + const size_t* indexConversion, + int rows, + size_t* diagIndices); - /** - * @brief Computes indices of diagonal elements for non-reordered GPU preconditioner. - * - * @tparam T Matrix element type - * @param[in] mat The matrix in BCSR format - * @param[in] rowIndices Row indices of the BCSR matrix - * @param[in] colIndices Column indices of the BCSR matrix - * @param[in] indexConversion Index conversion mapping - * @param[in] rows Number of rows in the matrix - * @param[out] diagIndices Array to store the indices of diagonal elements - */ - template - void computeDiagIndices(const T* mat, - const int* rowIndices, - const int* colIndices, - const int* reorderedToNatural, - int rows, - size_t* diagIndices); +/** + * @brief Computes indices of diagonal elements for non-reordered GPU preconditioner. + * + * @tparam T Matrix element type + * @param[in] mat The matrix in BCSR format + * @param[in] rowIndices Row indices of the BCSR matrix + * @param[in] colIndices Column indices of the BCSR matrix + * @param[in] indexConversion Index conversion mapping + * @param[in] rows Number of rows in the matrix + * @param[out] diagIndices Array to store the indices of diagonal elements + */ +template +void computeDiagIndices(const T* mat, + const int* rowIndices, + const int* colIndices, + const int* reorderedToNatural, + int rows, + size_t* diagIndices); -} +} // namespace Opm::gpuistl::detail #endif