GpuSparseMatrix.hpp
Go to the documentation of this file.
1/*
2 Copyright 2022-2023 SINTEF AS
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_GPUSPARSEMATRIX_HPP
20#define OPM_GPUSPARSEMATRIX_HPP
21#include <cusparse.h>
22#include <iostream>
23#include <memory>
24#include <stdexcept>
25#include <type_traits>
26#include <opm/common/ErrorMacros.hpp>
32#include <vector>
33
34namespace Opm::gpuistl
35{
36
57template <typename T>
59
60{
61public:
62 using field_type = T;
63
72 static constexpr int max_block_size = 6;
73
87 GpuSparseMatrix(const T* nonZeroElements,
88 const int* rowIndices,
89 const int* columnIndices,
90 size_t numberOfNonzeroBlocks,
91 size_t blockSize,
92 size_t numberOfRows);
93
105 const GpuVector<int>& columnIndices,
106 size_t blockSize);
107
109
110 // We want to have this as non-mutable as possible, that is we do not want
111 // to deal with changing matrix sizes and sparsity patterns.
113
115
125 template <class MatrixType>
126 static GpuSparseMatrix<T> fromMatrix(const MatrixType& matrix, bool copyNonZeroElementsDirectly = false);
127
132
137
142
143
148
152 size_t N() const
153 {
154 // Technically this safe conversion is not needed since we enforce these to be
155 // non-negative in the constructor, but keeping them for added sanity for now.
156 //
157 // We don't believe this will yield any performance penality (it's used too far away from the inner loop),
158 // but should that be false, they can be removed.
159 return detail::to_size_t(m_numberOfRows);
160 }
161
166 size_t nonzeroes() const
167 {
168 // Technically this safe conversion is not needed since we enforce these to be
169 // non-negative in the constructor, but keeping them for added sanity for now.
170 //
171 // We don't believe this will yield any performance penality (it's used too far away from the inner loop),
172 // but should that be false, they can be removed.
173 return detail::to_size_t(m_numberOfNonzeroBlocks);
174 }
175
182 {
183 if (m_genericMatrixForBlockSize1) {
184 return m_genericMatrixForBlockSize1->getNonZeroValues();
185 }
186 return m_nonZeroElements;
187 }
188
195 {
196 if (m_genericMatrixForBlockSize1) {
197 return m_genericMatrixForBlockSize1->getNonZeroValues();
198 }
199 return m_nonZeroElements;
200 }
201
208 {
209 if (m_genericMatrixForBlockSize1) {
210 return m_genericMatrixForBlockSize1->getRowIndices();
211 }
212 return m_rowIndices;
213 }
214
221 {
222 if (m_genericMatrixForBlockSize1) {
223 return m_genericMatrixForBlockSize1->getRowIndices();
224 }
225 return m_rowIndices;
226 }
227
234 {
235 if (m_genericMatrixForBlockSize1) {
236 return m_genericMatrixForBlockSize1->getColumnIndices();
237 }
238 return m_columnIndices;
239 }
240
247 {
248 if (m_genericMatrixForBlockSize1) {
249 return m_genericMatrixForBlockSize1->getColumnIndices();
250 }
251 return m_columnIndices;
252 }
253
260 size_t dim() const
261 {
262 // Technically this safe conversion is not needed since we enforce these to be
263 // non-negative in the constructor, but keeping them for added sanity for now.
264 //
265 // We don't believe this will yield any performance penality (it's used too far away from the inner loop),
266 // but should that be false, they can be removed.
267 return detail::to_size_t(m_blockSize) * detail::to_size_t(m_numberOfRows);
268 }
269
273 size_t blockSize() const
274 {
275 // Technically this safe conversion is not needed since we enforce these to be
276 // non-negative in the constructor, but keeping them for added sanity for now.
277 //
278 // We don't believe this will yield any performance penality (it's used too far away from the inner loop),
279 // but should that be false, they can be removed.
280 return detail::to_size_t(m_blockSize);
281 }
282
289 {
290 return *m_matrixDescription;
291 }
292
298 virtual void mv(const GpuVector<T>& x, GpuVector<T>& y) const;
299
305 virtual void umv(const GpuVector<T>& x, GpuVector<T>& y) const;
306
307
314 virtual void usmv(T alpha, const GpuVector<T>& x, GpuVector<T>& y) const;
315
326 template <class MatrixType>
327 void updateNonzeroValues(const MatrixType& matrix, bool copyNonZeroElementsDirectly = false);
328
335
336
357 template<class FunctionType>
358 auto dispatchOnBlocksize(FunctionType function) const
359 {
360 return dispatchOnBlocksizeImpl<max_block_size>(function);
361 }
362
363private:
364 GpuVector<T> m_nonZeroElements;
365 GpuVector<int> m_columnIndices;
366 GpuVector<int> m_rowIndices;
367
368 // Notice that we store these three as int to make sure we are cusparse compatible.
369 //
370 // This gives the added benefit of checking the size constraints at construction of the matrix
371 // rather than in some call to cusparse.
372 const int m_numberOfNonzeroBlocks;
373 const int m_numberOfRows;
374 const int m_blockSize;
375
376 detail::GpuSparseMatrixDescriptionPtr m_matrixDescription;
377 detail::CuSparseHandle& m_cusparseHandle;
378
379 // For blockSize == 1, we use the generic API
380 std::unique_ptr<GpuSparseMatrixGeneric<T>> m_genericMatrixForBlockSize1;
381
382 template <class VectorType>
383 void assertSameSize(const VectorType& vector) const;
384
385 template<int blockSizeCompileTime, class FunctionType>
386 auto dispatchOnBlocksizeImpl(FunctionType function) const
387 {
388 if (blockSizeCompileTime == m_blockSize) {
389 return function(std::integral_constant<int, blockSizeCompileTime>());
390 }
391
392 if constexpr (blockSizeCompileTime > 1) {
393 return dispatchOnBlocksizeImpl<blockSizeCompileTime - 1>(function);
394 } else {
395 OPM_THROW(std::runtime_error, fmt::format("Unsupported block size: {}", m_blockSize));
396 }
397 }
398};
399} // namespace Opm::gpuistl
400#endif
The GpuSparseMatrix class simple wrapper class for a CuSparse matrix.
Definition: GpuSparseMatrix.hpp:60
GpuSparseMatrix & operator=(const GpuSparseMatrix &)=delete
virtual void mv(const GpuVector< T > &x, GpuVector< T > &y) const
mv performs matrix vector multiply y = Ax
const GpuVector< int > & getRowIndices() const
getRowIndices returns the row indices used to represent the BSR structure.
Definition: GpuSparseMatrix.hpp:220
GpuSparseMatrix(const GpuVector< int > &rowIndices, const GpuVector< int > &columnIndices, size_t blockSize)
const GpuVector< int > & getColumnIndices() const
getColumnIndices returns the column indices used to represent the BSR structure.
Definition: GpuSparseMatrix.hpp:246
void setNonUnitDiagonal()
setNonUnitDiagonal sets the CuSparse flag that this has non-unit diagional.
T field_type
Definition: GpuSparseMatrix.hpp:62
void updateNonzeroValues(const MatrixType &matrix, bool copyNonZeroElementsDirectly=false)
updateNonzeroValues updates the non-zero values by using the non-zero values of the supplied matrix
void updateNonzeroValues(const GpuSparseMatrix< T > &matrix)
updateNonzeroValues updates the non-zero values by using the non-zero values of the supplied matrix
detail::GpuSparseMatrixDescription & getDescription()
getDescription the cusparse matrix description.
Definition: GpuSparseMatrix.hpp:288
GpuSparseMatrix(const GpuSparseMatrix &)
GpuVector< int > & getRowIndices()
getRowIndices returns the row indices used to represent the BSR structure.
Definition: GpuSparseMatrix.hpp:207
const GpuVector< T > & getNonZeroValues() const
getNonZeroValues returns the GPU vector containing the non-zero values (ordered by block)
Definition: GpuSparseMatrix.hpp:194
size_t nonzeroes() const
nonzeroes behaves as the Dune::BCRSMatrix::nonzeros() function and returns the number of non zero blo...
Definition: GpuSparseMatrix.hpp:166
size_t blockSize() const
blockSize size of the blocks
Definition: GpuSparseMatrix.hpp:273
void setUpperTriangular()
setUpperTriangular sets the CuSparse flag that this is an upper diagonal (with unit diagonal) matrix.
size_t N() const
N returns the number of rows (which is equal to the number of columns)
Definition: GpuSparseMatrix.hpp:152
virtual void umv(const GpuVector< T > &x, GpuVector< T > &y) const
umv computes y=Ax+y
void setLowerTriangular()
setLowerTriangular sets the CuSparse flag that this is an lower diagonal (with non-unit diagonal) mat...
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
static constexpr int max_block_size
Maximum block size supported by this implementation.
Definition: GpuSparseMatrix.hpp:72
auto dispatchOnBlocksize(FunctionType function) const
Dispatches a function based on the block size of the matrix.
Definition: GpuSparseMatrix.hpp:358
void setUnitDiagonal()
setUnitDiagonal sets the CuSparse flag that this has unit diagional.
virtual void usmv(T alpha, const GpuVector< T > &x, GpuVector< T > &y) const
umv computes y=alpha * Ax + y
GpuVector< int > & getColumnIndices()
getColumnIndices returns the column indices used to represent the BSR structure.
Definition: GpuSparseMatrix.hpp:233
GpuSparseMatrix(const T *nonZeroElements, const int *rowIndices, const int *columnIndices, size_t numberOfNonzeroBlocks, size_t blockSize, size_t numberOfRows)
GpuVector< T > & getNonZeroValues()
getNonZeroValues returns the GPU vector containing the non-zero values (ordered by block)
Definition: GpuSparseMatrix.hpp:181
size_t dim() const
dim returns the dimension of the vector space on which this matrix acts
Definition: GpuSparseMatrix.hpp:260
The CuSparseHandle class provides a singleton for the simulator universal cuSparseHandle.
Definition: CuSparseHandle.hpp:41
The CuSparseResource class wraps a CuSparse resource in a proper RAII pattern.
Definition: CuSparseResource.hpp:55
std::shared_ptr< CuSparseResource< cusparseMatDescr_t > > GpuSparseMatrixDescriptionPtr
Definition: CuMatrixDescription.hpp:35
__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
Definition: AmgxInterface.hpp:38