Loading...
Searching...
No Matches
GpuSparseMatrixWrapper.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_GPUSPARSEMATRIXWRAPPER_HPP
#define OPM_GPUSPARSEMATRIXWRAPPER_HPP
#include <cusparse.h>
#include <iostream>
#include <memory>
#include <stdexcept>
#include <type_traits>
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
#include <opm/simulators/linalg/gpuistl/GpuSparseMatrix.hpp>
#include <opm/simulators/linalg/gpuistl/GpuSparseMatrixGeneric.hpp>
#include <opm/simulators/linalg/gpuistl/detail/CuMatrixDescription.hpp>
#include <opm/simulators/linalg/gpuistl/detail/CuSparseHandle.hpp>
#include <opm/simulators/linalg/gpuistl/detail/safe_conversion.hpp>
#include <vector>
namespace Opm::gpuistl
{
template <typename T>
class GpuSparseMatrixWrapper
{
public:
using field_type = T;
/*
Here is the primary function of this class.
Since the generic API for CUDA/HIP is primaryly supported on CUDA 13 (and not yet HIP) for blocked
matrices, places wanting to use blocked matrices can invoke this class which handles which API to use.
Basically we just check if HIP is present, or if we are using cuda and a version prior to 13.
*/
#if USE_HIP || (!USE_HIP && CUDA_VERSION < 13000)
using matrix_type = GpuSparseMatrix<T>;
#else
using matrix_type = GpuSparseMatrixGeneric<T>;
#endif
// Arrow operator overloads for direct access to the underlying matrix
matrix_type* operator->() {
if (!m_matrix) {
throw std::runtime_error("GpuSparseMatrixWrapper: underlying matrix is nullptr.");
}
return m_matrix.get();
}
const matrix_type* operator->() const {
if (!m_matrix) {
throw std::runtime_error("GpuSparseMatrixWrapper: underlying matrix is nullptr.");
}
return m_matrix.get();
}
static constexpr int max_block_size = 6;
GpuSparseMatrixWrapper(const T* nonZeroElements,
const int* rowIndices,
const int* columnIndices,
size_t numberOfNonzeroBlocks,
size_t blockSize,
size_t numberOfRows)
{
m_matrix = std::make_unique<matrix_type>(nonZeroElements,
rowIndices,
columnIndices,
numberOfNonzeroBlocks,
blockSize,
numberOfRows);
}
GpuSparseMatrixWrapper(const GpuVector<int>& rowIndices,
const GpuVector<int>& columnIndices,
size_t blockSize)
{
m_matrix = std::make_unique<matrix_type>(rowIndices, columnIndices, blockSize);
}
GpuSparseMatrixWrapper(const GpuSparseMatrixWrapper& other)
{
if (!other.m_matrix) {
throw std::runtime_error("Internal error, other.m_matrix is a nullptr.");
}
m_matrix = std::make_unique<matrix_type>(*other.m_matrix);
}
// 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.
GpuSparseMatrixWrapper& operator=(const GpuSparseMatrixWrapper&) = delete;
~GpuSparseMatrixWrapper() = default;
GpuSparseMatrixWrapper() = default;
const matrix_type& get() const {
if (!m_matrix) {
throw std::runtime_error("GpuSparseMatrixWrapper: underlying matrix is nullptr.");
}
return *m_matrix;
}
template <class MatrixType>
static GpuSparseMatrixWrapper<T> fromMatrix(const MatrixType& matrix, bool copyNonZeroElementsDirectly = false)
{
GpuSparseMatrixWrapper<T> gpuSparseMatrixWrapper;
gpuSparseMatrixWrapper.m_matrix = std::make_unique<matrix_type>(
matrix_type::fromMatrix(matrix, copyNonZeroElementsDirectly));
return gpuSparseMatrixWrapper;
}
// Only participates in overload resolution when matrix_type == GpuSparseMatrix<T>
template <class M = matrix_type,
typename = std::enable_if_t<std::is_same_v<M, GpuSparseMatrix<T>>>>
void setUpperTriangular()
{
m_matrix->setUpperTriangular();
}
template <class M = matrix_type,
typename = std::enable_if_t<std::is_same_v<M, GpuSparseMatrix<T>>>>
void setLowerTriangular()
{
m_matrix->setLowerTriangular();
}
template <class M = matrix_type,
typename = std::enable_if_t<std::is_same_v<M, GpuSparseMatrix<T>>>>
void setUnitDiagonal()
{
m_matrix->setUnitDiagonal();
}
template <class M = matrix_type,
typename = std::enable_if_t<std::is_same_v<M, GpuSparseMatrix<T>>>>
void setNonUnitDiagonal()
{
m_matrix->setNonUnitDiagonal();
}
size_t N() const
{
return m_matrix->N();
}
size_t nonzeroes() const
{
return m_matrix->nonzeroes();
}
GpuVector<T>& getNonZeroValues()
{
return m_matrix->getNonZeroValues();
}
const GpuVector<T>& getNonZeroValues() const
{
return m_matrix->getNonZeroValues();
}
GpuVector<int>& getRowIndices()
{
return m_matrix->getRowIndices();
}
const GpuVector<int>& getRowIndices() const
{
return m_matrix->getRowIndices();
}
GpuVector<int>& getColumnIndices()
{
return m_matrix->getColumnIndices();
}
const GpuVector<int>& getColumnIndices() const
{
return m_matrix->getColumnIndices();
}
size_t dim() const
{
return m_matrix->dim();
}
size_t blockSize() const
{
return m_matrix->blockSize();
}
detail::GpuSparseMatrixDescription& getDescription()
{
return m_matrix->getDescription();
}
virtual void mv(const GpuVector<T>& x, GpuVector<T>& y) const
{
m_matrix->mv(x, y);
}
virtual void umv(const GpuVector<T>& x, GpuVector<T>& y) const
{
m_matrix->umv(x, y);
}
virtual void usmv(T alpha, const GpuVector<T>& x, GpuVector<T>& y) const
{
m_matrix->usmv(alpha, x, y);
}
template <class MatrixType>
void updateNonzeroValues(const MatrixType& matrix, bool copyNonZeroElementsDirectly = false)
{
m_matrix->updateNonzeroValues(matrix, copyNonZeroElementsDirectly);
}
void updateNonzeroValues(const GpuSparseMatrixWrapper<T>& matrix)
{
m_matrix->updateNonzeroValues(matrix.get());
}
template<class FunctionType>
auto dispatchOnBlocksize(FunctionType function) const
{
return m_matrix->dispatchOnBlocksize(function);
}
private:
std::unique_ptr<matrix_type> m_matrix = nullptr;
};
} // namespace Opm::gpuistl
#endif
Opm::gpuistl::GpuSparseMatrixGeneric
The GpuSparseMatrixGeneric class uses cuSPARSE Generic API for sparse matrix operations.
Definition GpuSparseMatrixGeneric.hpp:51
Opm::gpuistl::GpuSparseMatrixWrapper
The GpuSparseMatrixWrapper Checks CUDA/HIP version and dispatches a version either using the old or t...
Definition GpuSparseMatrixWrapper.hpp:61
Opm::gpuistl::GpuSparseMatrixWrapper::getRowIndices
GpuVector< int > & getRowIndices()
getRowIndices returns the row indices used to represent the BSR structure.
Definition GpuSparseMatrixWrapper.hpp:267
Opm::gpuistl::GpuSparseMatrixWrapper::usmv
virtual void usmv(T alpha, const GpuVector< T > &x, GpuVector< T > &y) const
umv computes y=alpha * Ax + y
Definition GpuSparseMatrixWrapper.hpp:358
Opm::gpuistl::GpuSparseMatrixWrapper::dispatchOnBlocksize
auto dispatchOnBlocksize(FunctionType function) const
Dispatches a function based on the block size of the matrix.
Definition GpuSparseMatrixWrapper.hpp:411
Opm::gpuistl::GpuSparseMatrixWrapper::updateNonzeroValues
void updateNonzeroValues(const GpuSparseMatrixWrapper< T > &matrix)
updateNonzeroValues updates the non-zero values by using the non-zero values of the supplied matrix
Definition GpuSparseMatrixWrapper.hpp:384
Opm::gpuistl::GpuSparseMatrixWrapper::getDescription
detail::GpuSparseMatrixDescription & getDescription()
getDescription the cusparse matrix description.
Definition GpuSparseMatrixWrapper.hpp:326
Opm::gpuistl::GpuSparseMatrixWrapper::fromMatrix
static GpuSparseMatrixWrapper< 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 GpuSparseMatrixWrapper.hpp:179
Opm::gpuistl::GpuSparseMatrixWrapper::getColumnIndices
GpuVector< int > & getColumnIndices()
getColumnIndices returns the column indices used to represent the BSR structure.
Definition GpuSparseMatrixWrapper.hpp:287
Opm::gpuistl::GpuSparseMatrixWrapper::mv
virtual void mv(const GpuVector< T > &x, GpuVector< T > &y) const
mv performs matrix vector multiply y = Ax
Definition GpuSparseMatrixWrapper.hpp:336
Opm::gpuistl::GpuSparseMatrixWrapper::setUnitDiagonal
void setUnitDiagonal()
setUnitDiagonal sets the CuSparse flag that this has unit diagional.
Definition GpuSparseMatrixWrapper.hpp:210
Opm::gpuistl::GpuSparseMatrixWrapper::getNonZeroValues
GpuVector< T > & getNonZeroValues()
getNonZeroValues returns the GPU vector containing the non-zero values (ordered by block)
Definition GpuSparseMatrixWrapper.hpp:247
Opm::gpuistl::GpuSparseMatrixWrapper::setLowerTriangular
void setLowerTriangular()
setLowerTriangular sets the CuSparse flag that this is an lower diagonal (with non-unit diagonal) mat...
Definition GpuSparseMatrixWrapper.hpp:200
Opm::gpuistl::GpuSparseMatrixWrapper< field_type >::blockSize
size_t blockSize() const
Definition GpuSparseMatrixWrapper.hpp:316
Opm::gpuistl::GpuSparseMatrixWrapper::getRowIndices
const GpuVector< int > & getRowIndices() const
getRowIndices returns the row indices used to represent the BSR structure.
Definition GpuSparseMatrixWrapper.hpp:277
Opm::gpuistl::GpuSparseMatrixWrapper::nonzeroes
size_t nonzeroes() const
nonzeroes behaves as the Dune::BCRSMatrix::nonzeros() function and returns the number of non zero blo...
Definition GpuSparseMatrixWrapper.hpp:237
Opm::gpuistl::GpuSparseMatrixWrapper::getNonZeroValues
const GpuVector< T > & getNonZeroValues() const
getNonZeroValues returns the GPU vector containing the non-zero values (ordered by block)
Definition GpuSparseMatrixWrapper.hpp:257
Opm::gpuistl::GpuSparseMatrixWrapper::GpuSparseMatrixWrapper
GpuSparseMatrixWrapper(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 GpuSparseMatrixWrapper.hpp:114
Opm::gpuistl::GpuSparseMatrixWrapper::getColumnIndices
const GpuVector< int > & getColumnIndices() const
getColumnIndices returns the column indices used to represent the BSR structure.
Definition GpuSparseMatrixWrapper.hpp:297
Opm::gpuistl::GpuSparseMatrixWrapper::umv
virtual void umv(const GpuVector< T > &x, GpuVector< T > &y) const
umv computes y=Ax+y
Definition GpuSparseMatrixWrapper.hpp:346
Opm::gpuistl::GpuSparseMatrixWrapper::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 GpuSparseMatrixWrapper.hpp:374
Opm::gpuistl::GpuSparseMatrixWrapper::setNonUnitDiagonal
void setNonUnitDiagonal()
setNonUnitDiagonal sets the CuSparse flag that this has non-unit diagional.
Definition GpuSparseMatrixWrapper.hpp:220
Opm::gpuistl::GpuSparseMatrixWrapper< field_type >::max_block_size
static constexpr int max_block_size
Definition GpuSparseMatrixWrapper.hpp:99
Opm::gpuistl::GpuSparseMatrixWrapper::N
size_t N() const
N returns the number of rows (which is equal to the number of columns).
Definition GpuSparseMatrixWrapper.hpp:228
Opm::gpuistl::GpuSparseMatrixWrapper::GpuSparseMatrixWrapper
GpuSparseMatrixWrapper(const GpuVector< int > &rowIndices, const GpuVector< int > &columnIndices, size_t blockSize)
Create a sparse matrix by copying the sparsity structure of another matrix, not filling in the values...
Definition GpuSparseMatrixWrapper.hpp:139
Opm::gpuistl::GpuSparseMatrixWrapper::dim
size_t dim() const
dim returns the dimension of the vector space on which this matrix acts
Definition GpuSparseMatrixWrapper.hpp:308
Opm::gpuistl::GpuSparseMatrix
The GpuSparseMatrix class simple wrapper class for a CuSparse matrix.
Definition GpuSparseMatrix.hpp:60
Opm::gpuistl::GpuSparseMatrix::fromMatrix
static GpuSparseMatrix< 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 GpuSparseMatrix.cpp:115
Opm::gpuistl::GpuVector
Definition gpu_type_detection.hpp:30
Opm::gpuistl::detail::GpuSparseMatrixDescription
CuSparseResource< cusparseMatDescr_t > GpuSparseMatrixDescription
GpuSparseMatrixDescription holder.
Definition CuMatrixDescription.hpp:30