GpuSeqILU0.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_GPUSEQILU0_HPP
20#define OPM_GPUSEQILU0_HPP
21
22#include <dune/istl/preconditioner.hh>
30
31
32
33namespace Opm::gpuistl
34{
50template <class M, class X, class Y, int l = 1>
52{
53public:
55 using matrix_type = typename std::remove_const<M>::type ;
57 using domain_type = X;
59 using range_type = Y;
61 using field_type = typename X::field_type;
62
69 template <typename T = M, typename std::enable_if_t<is_gpu_matrix_v<T>, int> = 0>
70 GpuSeqILU0(const M& A, field_type w);
71
72
79 template <typename T = M, typename std::enable_if_t<!is_gpu_matrix_v<T>, int> = 0>
80 GpuSeqILU0(const M& A, field_type w);
81
82
85 virtual void pre(X& x, Y& b) override;
86
88 virtual void apply(X& v, const Y& d) override;
89
92 virtual void post(X& x) override;
93
95 virtual Dune::SolverCategory::Category category() const override;
96
98 virtual void update() override;
99
100
102 static constexpr bool shouldCallPre()
103 {
104 return false;
105 }
106
108 static constexpr bool shouldCallPost()
109 {
110 return false;
111 }
112
113 virtual bool hasPerfectUpdate() const override {
114 return true;
115 }
116
117
118private:
120 const M& m_underlyingMatrix;
122 field_type m_w;
123
128
129 GpuVector<field_type> m_temporaryStorage;
130
131
137
138 std::unique_ptr<GpuVector<field_type>> m_buffer;
139 detail::CuSparseHandle& m_cuSparseHandle;
140
141 bool m_analysisDone = false;
142
143 void analyzeMatrix();
144 size_t findBufferSize();
145
146 void createILU();
147
148 void updateILUConfiguration();
149};
150
151
152// Case where matrix is a CPU matrix as input
153template <class M, class X, class Y, int l>
154template <typename T, typename std::enable_if_t<!is_gpu_matrix_v<T>, int>>
156 : m_underlyingMatrix(A)
157 , m_w(w)
158 , m_LU(GpuSparseMatrix<field_type>::fromMatrix(detail::makeMatrixWithNonzeroDiagonal(A)))
159 , m_temporaryStorage(m_LU.N() * m_LU.blockSize())
160 , m_descriptionL(detail::createLowerDiagonalDescription())
161 , m_descriptionU(detail::createUpperDiagonalDescription())
162 , m_cuSparseHandle(detail::CuSparseHandle::getInstance())
163{
164 // Some sanity check
165 OPM_ERROR_IF(A.N() != m_LU.N(),
166 fmt::format("CuSparse matrix not same size as DUNE matrix. {} vs {}.", m_LU.N(), A.N()));
167 OPM_ERROR_IF(
168 A[0][0].N() != m_LU.blockSize(),
169 fmt::format("CuSparse matrix not same blocksize as DUNE matrix. {} vs {}.", m_LU.blockSize(), A[0][0].N()));
170 OPM_ERROR_IF(
171 A.N() * A[0][0].N() != m_LU.dim(),
172 fmt::format("CuSparse matrix not same dimension as DUNE matrix. {} vs {}.", m_LU.dim(), A.N() * A[0][0].N()));
173 OPM_ERROR_IF(A.nonzeroes() != m_LU.nonzeroes(),
174 fmt::format("CuSparse matrix not same number of non zeroes as DUNE matrix. {} vs {}. ",
175 m_LU.nonzeroes(),
176 A.nonzeroes()));
177
178
179 updateILUConfiguration();
180}
181
182// Case where matrix is a GPU matrix as input
183template <class M, class X, class Y, int l>
184template <typename T, typename std::enable_if_t<is_gpu_matrix_v<T>, int>>
185GpuSeqILU0<M, X, Y, l>::GpuSeqILU0(const M& A, field_type w)
186 : m_underlyingMatrix(A)
187 , m_w(w)
188 , m_LU(A)
189 , m_temporaryStorage(m_LU.N() * m_LU.blockSize())
190 , m_descriptionL(detail::createLowerDiagonalDescription())
191 , m_descriptionU(detail::createUpperDiagonalDescription())
192 , m_cuSparseHandle(detail::CuSparseHandle::getInstance())
193{
194 // Some sanity check
195 OPM_ERROR_IF(A.N() != m_LU.N(),
196 fmt::format("CuSparse matrix not same size as DUNE matrix. {} vs {}.", m_LU.N(), A.N()));
197 OPM_ERROR_IF(A.nonzeroes() != m_LU.nonzeroes(),
198 fmt::format("CuSparse matrix not same number of non zeroes as DUNE matrix. {} vs {}. ",
199 m_LU.nonzeroes(),
200 A.nonzeroes()));
201
202
203 updateILUConfiguration();
204}
205
206} // end namespace Opm::gpuistl
207
208#endif
Interface class adding the update() method to the preconditioner interface.
Definition: PreconditionerWithUpdate.hpp:32
Sequential ILU0 preconditioner on the GPU through the CuSparse library.
Definition: GpuSeqILU0.hpp:52
virtual void pre(X &x, Y &b) override
Prepare the preconditioner.
typename std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: GpuSeqILU0.hpp:55
virtual void update() override
Updates the matrix data.
Y range_type
The range type of the preconditioner.
Definition: GpuSeqILU0.hpp:59
typename X::field_type field_type
The field type of the preconditioner.
Definition: GpuSeqILU0.hpp:61
virtual void post(X &x) override
Post processing.
virtual void apply(X &v, const Y &d) override
Apply the preconditoner.
X domain_type
The domain type of the preconditioner.
Definition: GpuSeqILU0.hpp:57
static constexpr bool shouldCallPost()
Definition: GpuSeqILU0.hpp:108
virtual Dune::SolverCategory::Category category() const override
Category of the preconditioner (see SolverCategory::Category)
virtual bool hasPerfectUpdate() const override
Definition: GpuSeqILU0.hpp:113
static constexpr bool shouldCallPre()
Definition: GpuSeqILU0.hpp:102
GpuSeqILU0(const M &A, field_type w)
Constructor.
Definition: GpuSeqILU0.hpp:155
size_t nonzeroes() const
nonzeroes behaves as the Dune::BCRSMatrix::nonzeros() function and returns the number of non zero blo...
Definition: GpuSparseMatrix.hpp:144
size_t blockSize() const
blockSize size of the blocks
Definition: GpuSparseMatrix.hpp:233
size_t N() const
N returns the number of rows (which is equal to the number of columns)
Definition: GpuSparseMatrix.hpp:130
size_t dim() const
dim returns the dimension of the vector space on which this matrix acts
Definition: GpuSparseMatrix.hpp:220
The CuSparseHandle class provides a singleton for the simulator universal cuSparseHandle.
Definition: CuSparseHandle.hpp:41
GpuSparseMatrixDescriptionPtr createLowerDiagonalDescription()
createLowerDiagonalDescription creates a lower diagonal matrix description
Definition: CuMatrixDescription.hpp:60
GpuSparseMatrixDescriptionPtr createUpperDiagonalDescription()
createUpperDiagonalDescription creates an upper diagonal matrix description
Definition: CuMatrixDescription.hpp:75
std::shared_ptr< CuSparseResource< cusparseMatDescr_t > > GpuSparseMatrixDescriptionPtr
Definition: CuMatrixDescription.hpp:35
const Matrix makeMatrixWithNonzeroDiagonal(const Matrix &matrix, const typename Matrix::field_type replacementValue=std::numeric_limits< typename Matrix::field_type >::epsilon())
makeMatrixWithNonzeroDiagonal creates a new matrix with the zero diagonal elements (when viewed as a ...
Definition: fix_zero_diagonal.hpp:40
Definition: autotuner.hpp:30