OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-11-03 12:47:42 -06:00
Code Issues Packages Projects Releases Wiki Activity

formatting

Browse Source
This commit is contained in:
Tobias Meyer Andersen committed 2025-06-26 13:47:48 +02:00
1 parent a232398dd3
commit 9d09812877
6 files changed
+306 -297
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

opm/simulators/linalg/gpuistl/GpuDILU.cpp
+100 -83
View File
@@ -21,28 +21,33 @@
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <fmt/core.h>
#include <functional>
#include <limits>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/TimingMacros.hpp>
#include <opm/simulators/linalg/GraphColoring.hpp>
#include <opm/simulators/linalg/gpuistl/detail/autotuner.hpp>
#include <opm/simulators/linalg/gpuistl/GpuDILU.hpp>
#include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
#include <opm/simulators/linalg/gpuistl/detail/autotuner.hpp>
#include <opm/simulators/linalg/gpuistl/detail/coloringAndReorderingUtils.hpp>
#include <opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.hpp>
#include <opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp>
#include <opm/simulators/linalg/matrixblock.hh>
#include <tuple>
#include <functional>
#include <utility>
#include <string>
#include <tuple>
#include <utility>
namespace Opm::gpuistl
{
template <class M, class X, class Y, int l>
GpuDILU<M, X, Y, l>::GpuDILU(const typename GpuDILU<M, X, Y, l>::GPUMatrix& gpuMatrix, const M& cpuMatrix, bool splitMatrix, bool tuneKernels, int mixedPrecisionScheme, bool reorder)
GpuDILU<M, X, Y, l>::GpuDILU(const typename GpuDILU<M, X, Y, l>::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<M, X, Y, l>::GpuDILU(const typename GpuDILU<M, X, Y, l>::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<M, X, Y, l>::GpuDILU(const typename GpuDILU<M, X, Y, l>::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<GpuVector<field_type>>(blocksize_ * blocksize_ * cpuMatrix.N());
std::tie(m_gpuMatrixReorderedLower, m_gpuMatrixReorderedUpper)
= detail::extractLowerAndUpperMatrices<M, field_type, GpuSparseMatrix<field_type>>(cpuMatrix,
m_reorderedToNatural);
}
else {
m_reorderedToNatural);
} else {
m_gpuMatrixReordered = detail::createReorderedMatrix<M, field_type, GpuSparseMatrix<field_type>>(
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<FloatMat>(m_gpuMatrixReorderedLower->getRowIndices(), m_gpuMatrixReorderedLower->getColumnIndices(), blocksize_);
m_gpuMatrixReorderedUpperFloat = std::make_unique<FloatMat>(m_gpuMatrixReorderedUpper->getRowIndices(), m_gpuMatrixReorderedUpper->getColumnIndices(), blocksize_);
m_gpuMatrixReorderedLowerFloat = std::make_unique<FloatMat>(
m_gpuMatrixReorderedLower->getRowIndices(), m_gpuMatrixReorderedLower->getColumnIndices(), blocksize_);
m_gpuMatrixReorderedUpperFloat = std::make_unique<FloatMat>(
m_gpuMatrixReorderedUpper->getRowIndices(), m_gpuMatrixReorderedUpper->getColumnIndices(), blocksize_);
if (m_mixedPrecisionScheme == MatrixStorageMPScheme::FLOAT_DIAG_FLOAT_OFFDIAG) {
m_gpuDInvFloat = std::make_unique<FloatVec>(m_gpuMatrix.N() * m_gpuMatrix.blockSize() * m_gpuMatrix.blockSize());
m_gpuDInvFloat
= std::make_unique<FloatVec>(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<field_type>(
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<field_type>(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<field_type>(
m_gpuMatrix.getNonZeroValues().data(),
m_gpuMatrix.getRowIndices().data(),
m_gpuMatrix.getColumnIndices().data(),
m_gpuLevelSets.data(),
m_gpuMatrix.N(),
m_diagIdxs.data());
} else {
detail::computeDiagIndicesNoReorder<field_type>(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 <class M, class X, class Y, int l>
@@ -154,7 +159,8 @@ GpuDILU<M, X, Y, l>::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<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int
lowerSolveThreadBlockSize,
m_stream.get());
} else {
detail::DILU::solveLowerLevelSetNoReorder<field_type, blocksize_>(
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<field_type, blocksize_>(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<M, X, Y, l>::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<blocksize_, field_type, float>(
m_gpuMatrixReorderedUpperFloat->getNonZeroValues().data(),
m_gpuMatrixReorderedUpperFloat->getRowIndices().data(),
@@ -262,7 +267,7 @@ GpuDILU<M, X, Y, l>::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<blocksize_, field_type, float>(
m_gpuMatrixReorderedUpperFloat->getNonZeroValues().data(),
m_gpuMatrixReorderedUpperFloat->getRowIndices().data(),
@@ -301,17 +306,16 @@ GpuDILU<M, X, Y, l>::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int
upperSolveThreadBlockSize,
m_stream.get());
} else {
detail::DILU::solveUpperLevelSetNoReorder<field_type, blocksize_>(
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<field_type, blocksize_>(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<M, X, Y, l>::reorderAndSplitMatrix(int moveThreadBlockSize)
moveThreadBlockSize);
} else {
detail::copyMatDataToReordered<field_type, blocksize_>(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<M, X, Y, l>::computeDiagonal(int factorizationBlockSize)
const int numOfRowsInLevel = m_levelSets[level].size();
if (m_splitMatrix) {
if (m_mixedPrecisionScheme == MatrixStorageMPScheme::FLOAT_DIAG_FLOAT_OFFDIAG) {
detail::DILU::computeDiluDiagonalSplit<blocksize_, field_type, float, MatrixStorageMPScheme::FLOAT_DIAG_FLOAT_OFFDIAG>(
detail::DILU::computeDiluDiagonalSplit<blocksize_,
field_type,
float,
MatrixStorageMPScheme::FLOAT_DIAG_FLOAT_OFFDIAG>(
m_gpuMatrixReorderedLower->getNonZeroValues().data(),
m_gpuMatrixReorderedLower->getRowIndices().data(),
m_gpuMatrixReorderedLower->getColumnIndices().data(),
@@ -417,7 +424,10 @@ GpuDILU<M, X, Y, l>::computeDiagonal(int factorizationBlockSize)
m_gpuMatrixReorderedUpperFloat->getNonZeroValues().data(),
factorizationBlockSize);
} else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG) {
detail::DILU::computeDiluDiagonalSplit<blocksize_, field_type, float, MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG>(
detail::DILU::computeDiluDiagonalSplit<blocksize_,
field_type,
float,
MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG>(
m_gpuMatrixReorderedLower->getNonZeroValues().data(),
m_gpuMatrixReorderedLower->getRowIndices().data(),
m_gpuMatrixReorderedLower->getColumnIndices().data(),
@@ -436,7 +446,10 @@ GpuDILU<M, X, Y, l>::computeDiagonal(int factorizationBlockSize)
factorizationBlockSize);
} else {
// TODO: should this be field type twice or field type then float in the template?
detail::DILU::computeDiluDiagonalSplit<blocksize_, field_type, float, MatrixStorageMPScheme::DOUBLE_DIAG_DOUBLE_OFFDIAG>(
detail::DILU::computeDiluDiagonalSplit<blocksize_,
field_type,
float,
MatrixStorageMPScheme::DOUBLE_DIAG_DOUBLE_OFFDIAG>(
m_gpuMatrixReorderedLower->getNonZeroValues().data(),
m_gpuMatrixReorderedLower->getRowIndices().data(),
m_gpuMatrixReorderedLower->getColumnIndices().data(),
@@ -490,41 +503,45 @@ GpuDILU<M, X, Y, l>::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<field_type> tmpV(m_gpuMatrix.N() * m_gpuMatrix.blockSize());
GpuVector<field_type> 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<Dune::BCRSMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>, \
::Opm::gpuistl::GpuVector<realtype>, \
::Opm::gpuistl::GpuVector<realtype>>; \
template class ::Opm::gpuistl::GpuDILU<Dune::BCRSMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>, \
::Opm::gpuistl::GpuVector<realtype>, \
::Opm::gpuistl::GpuVector<realtype>>
#define INSTANTIATE_CUDILU_DUNE(realtype, blockdim) \
template class ::Opm::gpuistl::GpuDILU<Dune::BCRSMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>, \
::Opm::gpuistl::GpuVector<realtype>, \
::Opm::gpuistl::GpuVector<realtype>>; \
template class ::Opm::gpuistl::GpuDILU<Dune::BCRSMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>, \
::Opm::gpuistl::GpuVector<realtype>, \
::Opm::gpuistl::GpuVector<realtype>>
INSTANTIATE_CUDILU_DUNE(double, 1);
INSTANTIATE_CUDILU_DUNE(double, 2);
opm/simulators/linalg/gpuistl/GpuDILU.hpp
+18 -11
View File
@@ -19,18 +19,18 @@
#ifndef OPM_GPUDILU_HPP
#define OPM_GPUDILU_HPP
#include <map>
#include <memory>
#include <opm/grid/utility/SparseTable.hpp>
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.hpp>
#include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
#include <opm/simulators/linalg/gpuistl/GpuBuffer.hpp>
#include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
#include <opm/simulators/linalg/gpuistl/detail/kernel_enums.hpp>
#include <opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.hpp>
#include <opm/simulators/linalg/gpuistl/gpu_resources.hpp>
#include <vector>
#include <map>
#include <utility>
#include <vector>
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<int> m_gpuReorderToNatural;
//! \brief Potentially stores the level sets on the GPU
GpuBuffer<size_t> m_gpuLevelSets; //TODO: make this an optional to not spend time on the copy if not needed
GpuBuffer<size_t> 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<field_type> 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<size_t> m_diagIdxs;
// Graphs for Apply
@@ -173,10 +180,10 @@ private:
std::map<std::pair<field_type*, const field_type*>, 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
opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu
+62 -66
View File
@@ -61,14 +61,14 @@ namespace
template <class T, int blocksize>
__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 <class T, int blocksize>
__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 <class T, int blocksize>
__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 <class T, int blocksize>
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<T, blocksize>, thrBlockSize);
int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(
cuSolveLowerLevelSetNoReorder<T, blocksize>, thrBlockSize);
int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
cuSolveLowerLevelSetNoReorder<T, blocksize><<<nThreadBlocks, threadBlockSize, 0, stream>>>(
reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, d, v);
@@ -491,18 +491,18 @@ solveUpperLevelSet(T* reorderedMat,
template <class T, int blocksize>
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<T, blocksize>, thrBlockSize);
int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(
cuSolveUpperLevelSetNoReorder<T, blocksize>, thrBlockSize);
int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
cuSolveUpperLevelSetNoReorder<T, blocksize><<<nThreadBlocks, threadBlockSize, 0, stream>>>(
reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v);
@@ -563,27 +563,21 @@ computeDiluDiagonal(T* reorderedMat,
template <class T, int blocksize>
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<T, blocksize>, thrBlockSize);
int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize);
cuComputeDiluDiagonalNoReorder<T, blocksize><<<nThreadBlocks, threadBlockSize>>>(reorderedMat,
rowIndices,
colIndices,
indexConversion,
diagIndices,
startIdx,
rowsInLevelSet,
dInv);
cuComputeDiluDiagonalNoReorder<T, blocksize><<<nThreadBlocks, threadBlockSize>>>(
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, blocksize>(T*, int*, int*, int*, int*, size_t*, const int, int, T*, int); \
template void computeDiluDiagonalNoReorder<T, blocksize>(const T*, const int*, const int*, const size_t*, const size_t*, int, int, T*, int); \
template void computeDiluDiagonal<T, blocksize>(T*, int*, int*, int*, int*, size_t*, const int, int, T*, int); \
template void computeDiluDiagonalNoReorder<T, blocksize>( \
const T*, const int*, const int*, const size_t*, const size_t*, int, int, T*, int); \
template void computeDiluDiagonalSplit<blocksize, T, double, MatrixStorageMPScheme::DOUBLE_DIAG_DOUBLE_OFFDIAG>( \
const T*, \
int*, \
@@ -739,10 +734,11 @@ computeDiluDiagonalSplit(const InputScalar* srcReorderedLowerMat,
double*, \
int); \
template void solveUpperLevelSet<T, blocksize>(T*, int*, int*, int*, int, int, const T*, T*, int, cudaStream_t); \
template void solveUpperLevelSetNoReorder<T, blocksize>(const T*, const int*, const int*, const size_t*, int, int, const T*, T*, int, cudaStream_t); \
template void solveUpperLevelSetNoReorder<T, blocksize>( \
const T*, const int*, const int*, const size_t*, int, int, const T*, T*, int, cudaStream_t); \
template void solveLowerLevelSet<T, blocksize>( \
T*, int*, int*, int*, int, int, const T*, const T*, T*, int, cudaStream_t); \
template void solveLowerLevelSetNoReorder<T, blocksize>( \
T*, int*, int*, int*, int, int, const T*, const T*, T*, int, cudaStream_t); \
template void solveLowerLevelSetNoReorder<T, blocksize>( \
const T*, const int*, const int*, const size_t*, int, int, const T*, const T*, T*, int, cudaStream_t);
INSTANTIATE_KERNEL_WRAPPERS(float, 1);
opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp
+11 -11
View File
@@ -59,16 +59,16 @@ void solveLowerLevelSet(T* reorderedMat,
template <class T, int blocksize>
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
opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.cu
+89 -100
View File
@@ -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 <class T>
__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 <class T>
__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 <class T>
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 <class T>
__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<T>
);
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<T><<<threadBlockSize, nThreadBlocks>>>(
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 <class T>
void computeDiagIndices(const T* mat,
} // namespace
template <class T>
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<T>
);
{
// Have cuda/hip automatically pick a reasonable block size for this kernel based on static analysis
int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuComputeDiagIndicesNoReorder<T>);
// 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<T><<<threadBlockSize, nThreadBlocks>>>(
mat, rowIndices, colIndices, reorderedToNatural, rows, diagIndices
);
}
// launch kernel
cuComputeDiagIndicesNoReorder<T>
<<<threadBlockSize, nThreadBlocks>>>(mat, rowIndices, colIndices, indexConversion, rows, diagIndices);
}
template <class T>
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<T>);
// Calculate the number of blocks needed
int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rows, threadBlockSize);
// launch kernel
cuComputeDiagIndices<T>
<<<threadBlockSize, nThreadBlocks>>>(mat, rowIndices, colIndices, reorderedToNatural, rows, diagIndices);
}
} // namespace Opm::gpuistl::detail
#define DEFINE_KERNEL_WRAPPERS_FOR_TYPES(T) \
template void Opm::gpuistl::detail::computeDiagIndicesNoReorder<T>(const T*, \
const int*, \
const int*, \
const size_t*, \
int, \
size_t*); \
template void Opm::gpuistl::detail::computeDiagIndices<T>(const T*, \
const int*, \
const int*, \
const int*, \
int, \
size_t*);
#define DEFINE_KERNEL_WRAPPERS_FOR_TYPES(T) \
template void Opm::gpuistl::detail::computeDiagIndicesNoReorder<T>( \
const T*, const int*, const int*, const size_t*, int, size_t*); \
template void Opm::gpuistl::detail::computeDiagIndices<T>( \
const T*, const int*, const int*, const int*, int, size_t*);
DEFINE_KERNEL_WRAPPERS_FOR_TYPES(float)
DEFINE_KERNEL_WRAPPERS_FOR_TYPES(double)
opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU_variants_helper_kernels.hpp
+26 -26
View File
@@ -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 <class T>
void computeDiagIndicesNoReorder(const T* mat,
const int* rowIndices,
const int* colIndices,
const size_t* indexConversion,
int rows,
size_t* diagIndices);
template <class T>
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 <class T>
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 <class T>
void computeDiagIndices(const T* mat,
const int* rowIndices,
const int* colIndices,
const int* reorderedToNatural,
int rows,
size_t* diagIndices);
}
} // namespace Opm::gpuistl::detail
#endif