Loading...
Searching...
No Matches
GpuSparseMatrixGeneric.hpp
/*
Copyright 2025 Equinor ASA
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_GPUSPARSEMATRIXGENERIC_HPP
#define OPM_GPUSPARSEMATRIXGENERIC_HPP
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/linalg/gpuistl/GpuBuffer.hpp>
#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
#include <opm/simulators/linalg/gpuistl/detail/CuSparseHandle.hpp>
#include <opm/simulators/linalg/gpuistl/detail/cusparse_safe_call.hpp>
#include <opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_utilities.hpp>
#include <opm/simulators/linalg/gpuistl/detail/safe_conversion.hpp>
#include <cusparse.h>
#include <cstddef>
namespace Opm::gpuistl
{
template <typename T>
class GpuSparseMatrixGeneric
{
public:
using field_type = T;
static constexpr int max_block_size = 6;
GpuSparseMatrixGeneric(const T* nonZeroElements,
const int* rowIndices,
const int* columnIndices,
size_t numberOfNonzeroBlocks,
size_t blockSize,
size_t numberOfRows);
GpuSparseMatrixGeneric(const GpuVector<int>& rowIndices, const GpuVector<int>& columnIndices, size_t blockSize);
GpuSparseMatrixGeneric(const GpuSparseMatrixGeneric&);
void preprocessSpMV();
// We want to have this as non-mutable as possible, that is we do not want
// to deal with changing matrix sizes and sparsity patterns.
GpuSparseMatrixGeneric& operator=(const GpuSparseMatrixGeneric&) = delete;
virtual ~GpuSparseMatrixGeneric();
template <class MatrixType>
static GpuSparseMatrixGeneric<T> fromMatrix(const MatrixType& matrix, bool copyNonZeroElementsDirectly = false);
size_t N() const
{
return detail::to_size_t(m_numberOfRows);
}
size_t nonzeroes() const
{
// Technically this safe conversion is not needed since we enforce these to be
// non-negative in the constructor, but keeping them for added sanity for now.
//
// We don't believe this will yield any performance penality (it's used too far away from the inner loop),
// but should that be false, they can be removed.
return detail::to_size_t(m_numberOfNonzeroBlocks);
}
GpuVector<T>& getNonZeroValues()
{
return m_nonZeroElements;
}
const GpuVector<T>& getNonZeroValues() const
{
return m_nonZeroElements;
}
GpuVector<int>& getRowIndices()
{
return m_rowIndices;
}
const GpuVector<int>& getRowIndices() const
{
return m_rowIndices;
}
GpuVector<int>& getColumnIndices()
{
return m_columnIndices;
}
const GpuVector<int>& getColumnIndices() const
{
return m_columnIndices;
}
size_t dim() const
{
// Technically this safe conversion is not needed since we enforce these to be
// non-negative in the constructor, but keeping them for added sanity for now.
//
// We don't believe this will yield any performance penality (it's used too far away from the inner loop),
// but should that be false, they can be removed.
return detail::to_size_t(m_blockSize) * detail::to_size_t(m_numberOfRows);
}
size_t blockSize() const
{
// Technically this safe conversion is not needed since we enforce these to be
// non-negative in the constructor, but keeping them for added sanity for now.
//
// We don't believe this will yield any performance penality (it's used too far away from the inner loop),
// but should that be false, they can be removed.
return detail::to_size_t(m_blockSize);
}
virtual void mv(const GpuVector<T>& x, GpuVector<T>& y) const;
virtual void umv(const GpuVector<T>& x, GpuVector<T>& y) const;
virtual void usmv(T alpha, const GpuVector<T>& x, GpuVector<T>& y) const;
template <class MatrixType>
void updateNonzeroValues(const MatrixType& matrix, bool copyNonZeroElementsDirectly = false);
void updateNonzeroValues(const GpuSparseMatrixGeneric<T>& matrix);
template<class FunctionType>
auto dispatchOnBlocksize(FunctionType function) const
{
return dispatchOnBlocksizeImpl<max_block_size>(function);
}
private:
GpuVector<T> m_nonZeroElements;
GpuVector<int> m_columnIndices;
GpuVector<int> m_rowIndices;
// Notice that we store these three as int to make sure we are cusparse compatible.
//
// This gives the added benefit of checking the size constraints at construction of the matrix
// rather than in some call to cusparse.
const int m_numberOfNonzeroBlocks;
const int m_numberOfRows;
const int m_blockSize;
// Generic API descriptors
decltype(detail::makeSafeMatrixDescriptor()) m_matrixDescriptor;
detail::CuSparseHandle& m_cusparseHandle;
// Cached buffer for operations
mutable GpuBuffer<std::byte> m_buffer;
// Helper methods for SpMV operations
void spMV(T alpha, const GpuVector<T>& x, T beta, GpuVector<T>& y) const;
// Initialize matrix descriptor based on block size
void initializeMatrixDescriptor();
template <class VectorType>
void assertSameSize(const VectorType& vector) const;
// Helper to get cuSPARSE data type from C++ type
constexpr cudaDataType getDataType() const
{
if constexpr (std::is_same_v<T, float>) {
return CUDA_R_32F;
} else if constexpr (std::is_same_v<T, double>) {
return CUDA_R_64F;
} else {
static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>, "Only float and double are supported");
return CUDA_R_32F; // Unreachable, but needed to compile
}
}
template<int blockSizeCompileTime, class FunctionType>
auto dispatchOnBlocksizeImpl(FunctionType function) const
{
if (blockSizeCompileTime == m_blockSize) {
return function(std::integral_constant<int, blockSizeCompileTime>());
}
if constexpr (blockSizeCompileTime > 1) {
return dispatchOnBlocksizeImpl<blockSizeCompileTime - 1>(function);
} else {
OPM_THROW(std::runtime_error, fmt::format("Unsupported block size: {}", m_blockSize));
}
}
};
} // namespace Opm::gpuistl
#endif
Opm::gpuistl::GpuSparseMatrixGeneric
The GpuSparseMatrixGeneric class uses cuSPARSE Generic API for sparse matrix operations.
Definition GpuSparseMatrixGeneric.hpp:51
Opm::gpuistl::GpuSparseMatrixGeneric::dim
size_t dim() const
dim returns the dimension of the vector space on which this matrix acts
Definition GpuSparseMatrixGeneric.hpp:204
Opm::gpuistl::GpuSparseMatrixGeneric::getColumnIndices
const GpuVector< int > & getColumnIndices() const
getColumnIndices returns the column indices used to represent the BSR structure.
Definition GpuSparseMatrixGeneric.hpp:193
Opm::gpuistl::GpuSparseMatrixGeneric::blockSize
size_t blockSize() const
blockSize size of the blocks
Definition GpuSparseMatrixGeneric.hpp:217
Opm::gpuistl::GpuSparseMatrixGeneric::updateNonzeroValues
void updateNonzeroValues(const MatrixType &matrix, bool copyNonZeroElementsDirectly=false)
updateNonzeroValues updates the non-zero values by using the non-zero values of the supplied matrix
Definition GpuSparseMatrixGeneric.cpp:231
Opm::gpuistl::GpuSparseMatrixGeneric::getColumnIndices
GpuVector< int > & getColumnIndices()
getColumnIndices returns the column indices used to represent the BSR structure.
Definition GpuSparseMatrixGeneric.hpp:183
Opm::gpuistl::GpuSparseMatrixGeneric::getNonZeroValues
const GpuVector< T > & getNonZeroValues() const
getNonZeroValues returns the GPU vector containing the non-zero values (ordered by block)
Definition GpuSparseMatrixGeneric.hpp:153
Opm::gpuistl::GpuSparseMatrixGeneric::getNonZeroValues
GpuVector< T > & getNonZeroValues()
getNonZeroValues returns the GPU vector containing the non-zero values (ordered by block)
Definition GpuSparseMatrixGeneric.hpp:143
Opm::gpuistl::GpuSparseMatrixGeneric::getRowIndices
GpuVector< int > & getRowIndices()
getRowIndices returns the row indices used to represent the BSR structure.
Definition GpuSparseMatrixGeneric.hpp:163
Opm::gpuistl::GpuSparseMatrixGeneric::nonzeroes
size_t nonzeroes() const
nonzeroes behaves as the Dune::BCRSMatrix::nonzeros() function and returns the number of non zero blo...
Definition GpuSparseMatrixGeneric.hpp:128
Opm::gpuistl::GpuSparseMatrixGeneric::umv
virtual void umv(const GpuVector< T > &x, GpuVector< T > &y) const
umv computes y=Ax+y
Definition GpuSparseMatrixGeneric.cpp:299
Opm::gpuistl::GpuSparseMatrixGeneric::mv
virtual void mv(const GpuVector< T > &x, GpuVector< T > &y) const
mv performs matrix vector multiply y = Ax
Definition GpuSparseMatrixGeneric.cpp:290
Opm::gpuistl::GpuSparseMatrixGeneric::GpuSparseMatrixGeneric
GpuSparseMatrixGeneric(const T *nonZeroElements, const int *rowIndices, const int *columnIndices, size_t numberOfNonzeroBlocks, size_t blockSize, size_t numberOfRows)
Create the sparse matrix specified by the raw data.
Definition GpuSparseMatrixGeneric.cpp:42
Opm::gpuistl::GpuSparseMatrixGeneric::preprocessSpMV
void preprocessSpMV()
Preprocess SpMV operation to optimize for sparsity pattern.
Definition GpuSparseMatrixGeneric.cpp:140
Opm::gpuistl::GpuSparseMatrixGeneric::max_block_size
static constexpr int max_block_size
Maximum block size supported by this implementation.
Definition GpuSparseMatrixGeneric.hpp:63
Opm::gpuistl::GpuSparseMatrixGeneric::N
size_t N() const
N returns the number of rows (which is equal to the number of columns).
Definition GpuSparseMatrixGeneric.hpp:119
Opm::gpuistl::GpuSparseMatrixGeneric::dispatchOnBlocksize
auto dispatchOnBlocksize(FunctionType function) const
Dispatches a function based on the block size of the matrix.
Definition GpuSparseMatrixGeneric.hpp:291
Opm::gpuistl::GpuSparseMatrixGeneric::usmv
virtual void usmv(T alpha, const GpuVector< T > &x, GpuVector< T > &y) const
umv computes y=alpha * Ax + y
Definition GpuSparseMatrixGeneric.cpp:308
Opm::gpuistl::GpuSparseMatrixGeneric::fromMatrix
static GpuSparseMatrixGeneric< T > fromMatrix(const MatrixType &matrix, bool copyNonZeroElementsDirectly=false)
fromMatrix creates a new matrix with the same block size and values as the given matrix
Definition GpuSparseMatrixGeneric.cpp:202
Opm::gpuistl::GpuSparseMatrixGeneric::getRowIndices
const GpuVector< int > & getRowIndices() const
getRowIndices returns the row indices used to represent the BSR structure.
Definition GpuSparseMatrixGeneric.hpp:173
Opm::gpuistl::GpuVector
Definition gpu_type_detection.hpp:30
Opm::gpuistl::detail::CuSparseHandle
The CuSparseHandle class provides a singleton for the simulator universal cuSparseHandle.
Definition CuSparseHandle.hpp:41
Opm::gpuistl::detail::makeSafeMatrixDescriptor
auto makeSafeMatrixDescriptor()
Create RAII-managed cuSPARSE sparse matrix descriptor.
Definition gpusparse_matrix_utilities.hpp:56
Opm::gpuistl::detail::to_size_t
__host__ __device__ std::size_t to_size_t(int i)
to_size_t converts a (on most relevant platforms) a 32 bit signed int to a 64 bits unsigned int
Definition safe_conversion.hpp:86