GpuSparseMatrixWrapper.hpp
Go to the documentation of this file.
1/*
2 Copyright 2025 Equinor ASA
3
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#ifndef OPM_GPUSPARSEMATRIXWRAPPER_HPP
20#define OPM_GPUSPARSEMATRIXWRAPPER_HPP
21#include <cusparse.h>
22#include <iostream>
23#include <memory>
24#include <stdexcept>
25#include <type_traits>
26#include <opm/common/ErrorMacros.hpp>
33#include <vector>
34
35namespace Opm::gpuistl
36{
37
58template <typename T>
60
61{
62public:
63 using field_type = T;
64
65 /*
66 Here is the primary function of this class.
67 Since the generic API for CUDA/HIP is primaryly supported on CUDA 13 (and not yet HIP) for blocked
68 matrices, places wanting to use blocked matrices can invoke this class which handles which API to use.
69 Basically we just check if HIP is present, or if we are using cuda and a version prior to 13.
70 */
71#if USE_HIP || (!USE_HIP && CUDA_VERSION < 13000)
73#else
75#endif
76
77 // Arrow operator overloads for direct access to the underlying matrix
79 if (!m_matrix) {
80 throw std::runtime_error("GpuSparseMatrixWrapper: underlying matrix is nullptr.");
81 }
82 return m_matrix.get();
83 }
84 const matrix_type* operator->() const {
85 if (!m_matrix) {
86 throw std::runtime_error("GpuSparseMatrixWrapper: underlying matrix is nullptr.");
87 }
88 return m_matrix.get();
89 }
90
99 static constexpr int max_block_size = 6;
100
114 GpuSparseMatrixWrapper(const T* nonZeroElements,
115 const int* rowIndices,
116 const int* columnIndices,
117 size_t numberOfNonzeroBlocks,
118 size_t blockSize,
119 size_t numberOfRows)
120 {
121 m_matrix = std::make_unique<matrix_type>(nonZeroElements,
122 rowIndices,
123 columnIndices,
124 numberOfNonzeroBlocks,
125 blockSize,
126 numberOfRows);
127 }
128
140 const GpuVector<int>& columnIndices,
141 size_t blockSize)
142 {
143 m_matrix = std::make_unique<matrix_type>(rowIndices, columnIndices, blockSize);
144 }
145
147 {
148 if (!other.m_matrix) {
149 throw std::runtime_error("Internal error, other.m_matrix is a nullptr.");
150 }
151 m_matrix = std::make_unique<matrix_type>(*other.m_matrix);
152 }
153
154 // We want to have this as non-mutable as possible, that is we do not want
155 // to deal with changing matrix sizes and sparsity patterns.
157
159
161
162 const matrix_type& get() const {
163 if (!m_matrix) {
164 throw std::runtime_error("GpuSparseMatrixWrapper: underlying matrix is nullptr.");
165 }
166 return *m_matrix;
167 }
168
178 template <class MatrixType>
179 static GpuSparseMatrixWrapper<T> fromMatrix(const MatrixType& matrix, bool copyNonZeroElementsDirectly = false)
180 {
181 GpuSparseMatrixWrapper<T> gpuSparseMatrixWrapper;
182 gpuSparseMatrixWrapper.m_matrix = std::make_unique<matrix_type>(
183 matrix_type::fromMatrix(matrix, copyNonZeroElementsDirectly));
184 return gpuSparseMatrixWrapper;
185 }
186
187 // Only participates in overload resolution when matrix_type == GpuSparseMatrix<T>
188 template <class M = matrix_type,
189 typename = std::enable_if_t<std::is_same_v<M, GpuSparseMatrix<T>>>>
191 {
192 m_matrix->setUpperTriangular();
193 }
194
198 template <class M = matrix_type,
199 typename = std::enable_if_t<std::is_same_v<M, GpuSparseMatrix<T>>>>
201 {
202 m_matrix->setLowerTriangular();
203 }
204
208 template <class M = matrix_type,
209 typename = std::enable_if_t<std::is_same_v<M, GpuSparseMatrix<T>>>>
211 {
212 m_matrix->setUnitDiagonal();
213 }
214
218 template <class M = matrix_type,
219 typename = std::enable_if_t<std::is_same_v<M, GpuSparseMatrix<T>>>>
221 {
222 m_matrix->setNonUnitDiagonal();
223 }
224
228 size_t N() const
229 {
230 return m_matrix->N();
231 }
232
237 size_t nonzeroes() const
238 {
239 return m_matrix->nonzeroes();
240 }
241
248 {
249 return m_matrix->getNonZeroValues();
250 }
251
258 {
259 return m_matrix->getNonZeroValues();
260 }
261
268 {
269 return m_matrix->getRowIndices();
270 }
271
278 {
279 return m_matrix->getRowIndices();
280 }
281
288 {
289 return m_matrix->getColumnIndices();
290 }
291
298 {
299 return m_matrix->getColumnIndices();
300 }
301
308 size_t dim() const
309 {
310 return m_matrix->dim();
311 }
312
316 size_t blockSize() const
317 {
318 return m_matrix->blockSize();
319 }
320
327 {
328 return m_matrix->getDescription();
329 }
330
336 virtual void mv(const GpuVector<T>& x, GpuVector<T>& y) const
337 {
338 m_matrix->mv(x, y);
339 }
340
346 virtual void umv(const GpuVector<T>& x, GpuVector<T>& y) const
347 {
348 m_matrix->umv(x, y);
349 }
350
351
357 virtual void usmv(T alpha, const GpuVector<T>& x, GpuVector<T>& y) const
358 {
359 m_matrix->usmv(alpha, x, y);
360 }
361
372 template <class MatrixType>
373 void updateNonzeroValues(const MatrixType& matrix, bool copyNonZeroElementsDirectly = false)
374 {
375 m_matrix->updateNonzeroValues(matrix, copyNonZeroElementsDirectly);
376 }
377
384 {
385 m_matrix->updateNonzeroValues(matrix.get());
386 }
387
388
409 template<class FunctionType>
410 auto dispatchOnBlocksize(FunctionType function) const
411 {
412 return m_matrix->dispatchOnBlocksize(function);
413 }
414
415private:
416
417 std::unique_ptr<matrix_type> m_matrix = nullptr;
418};
419} // namespace Opm::gpuistl
420#endif
The GpuSparseMatrixGeneric class uses cuSPARSE Generic API for sparse matrix operations.
Definition: GpuSparseMatrixGeneric.hpp:51
The GpuSparseMatrix class simple wrapper class for a CuSparse matrix.
Definition: GpuSparseMatrix.hpp:60
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
The GpuSparseMatrixWrapper Checks CUDA/HIP version and dispatches a version either using the old or t...
Definition: GpuSparseMatrixWrapper.hpp:61
GpuVector< int > & getRowIndices()
getRowIndices returns the row indices used to represent the BSR structure.
Definition: GpuSparseMatrixWrapper.hpp:267
virtual void usmv(T alpha, const GpuVector< T > &x, GpuVector< T > &y) const
umv computes y=alpha * Ax + y
Definition: GpuSparseMatrixWrapper.hpp:357
GpuSparseMatrixWrapper & operator=(const GpuSparseMatrixWrapper &)=delete
auto dispatchOnBlocksize(FunctionType function) const
Dispatches a function based on the block size of the matrix.
Definition: GpuSparseMatrixWrapper.hpp:410
matrix_type * operator->()
Definition: GpuSparseMatrixWrapper.hpp:78
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:383
detail::GpuSparseMatrixDescription & getDescription()
getDescription the cusparse matrix description.
Definition: GpuSparseMatrixWrapper.hpp:326
GpuSparseMatrix< T > matrix_type
Definition: GpuSparseMatrixWrapper.hpp:72
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
GpuVector< int > & getColumnIndices()
getColumnIndices returns the column indices used to represent the BSR structure.
Definition: GpuSparseMatrixWrapper.hpp:287
virtual void mv(const GpuVector< T > &x, GpuVector< T > &y) const
mv performs matrix vector multiply y = Ax
Definition: GpuSparseMatrixWrapper.hpp:336
void setUnitDiagonal()
setUnitDiagonal sets the CuSparse flag that this has unit diagional.
Definition: GpuSparseMatrixWrapper.hpp:210
GpuVector< T > & getNonZeroValues()
getNonZeroValues returns the GPU vector containing the non-zero values (ordered by block)
Definition: GpuSparseMatrixWrapper.hpp:247
void setLowerTriangular()
setLowerTriangular sets the CuSparse flag that this is an lower diagonal (with non-unit diagonal) mat...
Definition: GpuSparseMatrixWrapper.hpp:200
GpuSparseMatrixWrapper(const GpuSparseMatrixWrapper &other)
Definition: GpuSparseMatrixWrapper.hpp:146
size_t blockSize() const
blockSize size of the blocks
Definition: GpuSparseMatrixWrapper.hpp:316
const GpuVector< int > & getRowIndices() const
getRowIndices returns the row indices used to represent the BSR structure.
Definition: GpuSparseMatrixWrapper.hpp:277
void setUpperTriangular()
Definition: GpuSparseMatrixWrapper.hpp:190
size_t nonzeroes() const
nonzeroes behaves as the Dune::BCRSMatrix::nonzeros() function and returns the number of non zero blo...
Definition: GpuSparseMatrixWrapper.hpp:237
const GpuVector< T > & getNonZeroValues() const
getNonZeroValues returns the GPU vector containing the non-zero values (ordered by block)
Definition: GpuSparseMatrixWrapper.hpp:257
GpuSparseMatrixWrapper(const T *nonZeroElements, const int *rowIndices, const int *columnIndices, size_t numberOfNonzeroBlocks, size_t blockSize, size_t numberOfRows)
Definition: GpuSparseMatrixWrapper.hpp:114
const GpuVector< int > & getColumnIndices() const
getColumnIndices returns the column indices used to represent the BSR structure.
Definition: GpuSparseMatrixWrapper.hpp:297
virtual void umv(const GpuVector< T > &x, GpuVector< T > &y) const
umv computes y=Ax+y
Definition: GpuSparseMatrixWrapper.hpp:346
const matrix_type * operator->() const
Definition: GpuSparseMatrixWrapper.hpp:84
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:373
void setNonUnitDiagonal()
setNonUnitDiagonal sets the CuSparse flag that this has non-unit diagional.
Definition: GpuSparseMatrixWrapper.hpp:220
const matrix_type & get() const
Definition: GpuSparseMatrixWrapper.hpp:162
static constexpr int max_block_size
Maximum block size supported by this implementation.
Definition: GpuSparseMatrixWrapper.hpp:99
size_t N() const
N returns the number of rows (which is equal to the number of columns)
Definition: GpuSparseMatrixWrapper.hpp:228
T field_type
Definition: GpuSparseMatrixWrapper.hpp:63
GpuSparseMatrixWrapper(const GpuVector< int > &rowIndices, const GpuVector< int > &columnIndices, size_t blockSize)
Definition: GpuSparseMatrixWrapper.hpp:139
size_t dim() const
dim returns the dimension of the vector space on which this matrix acts
Definition: GpuSparseMatrixWrapper.hpp:308
The CuSparseResource class wraps a CuSparse resource in a proper RAII pattern.
Definition: CuSparseResource.hpp:55
Definition: AmgxInterface.hpp:38