My Project
Loading...
Searching...
No Matches
coloringAndReorderingUtils.hpp
1/*
2 Copyright 2024 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_COLORING_AND_REORDERING_UTILS_HPP
20#define OPM_COLORING_AND_REORDERING_UTILS_HPP
21
22#include <fmt/core.h>
23#include <memory>
24#include <opm/common/ErrorMacros.hpp>
25#include <opm/grid/utility/SparseTable.hpp>
26#include <opm/simulators/linalg/gpuistl/detail/safe_conversion.hpp>
27#include <tuple>
28#include <vector>
29/*
30This file contains a collection of utility functions used in the GPU implementation of ILU and DILU
31The functions deal with creating the mappings between reordered and natural indices, as well as
32extracting sparsity structures from dune matrices and creating gpusparsematrix indices
33*/
35{
36inline std::vector<int>
37createReorderedToNatural(const Opm::SparseTable<size_t>& levelSets)
38{
39 auto res = std::vector<int>(Opm::gpuistl::detail::to_size_t(levelSets.dataSize()));
40 int globCnt = 0;
41 for (auto row : levelSets) {
42 for (auto col : row) {
44 fmt::format("Internal error. globCnt = {}, res.size() = {}", globCnt, res.size()));
45 res[globCnt++] = static_cast<int>(col);
46 }
47 }
48 return res;
49}
50
51inline std::vector<int>
52createNaturalToReordered(const Opm::SparseTable<size_t>& levelSets)
53{
54 auto res = std::vector<int>(Opm::gpuistl::detail::to_size_t(levelSets.dataSize()));
55 int globCnt = 0;
56 for (auto row : levelSets) {
57 for (auto col : row) {
59 fmt::format("Internal error. globCnt = {}, res.size() = {}", globCnt, res.size()));
60 res[col] = globCnt++;
61 }
62 }
63 return res;
64}
65
66template <class M, class field_type, class GPUM>
67inline std::unique_ptr<GPUM>
68createReorderedMatrix(const M& naturalMatrix, const std::vector<int>& reorderedToNatural)
69{
70 M reorderedMatrix(naturalMatrix.N(), naturalMatrix.N(), naturalMatrix.nonzeroes(), M::row_wise);
71 for (auto dstRowIt = reorderedMatrix.createbegin(); dstRowIt != reorderedMatrix.createend(); ++dstRowIt) {
72 auto srcRow = naturalMatrix.begin() + reorderedToNatural[dstRowIt.index()];
73 // For elements in A
74 for (auto elem = srcRow->begin(); elem != srcRow->end(); elem++) {
75 dstRowIt.insert(elem.index());
76 }
77 }
78 return std::unique_ptr<GPUM>(new auto(GPUM::fromMatrix(reorderedMatrix, true)));
79}
80
81template <class M, class field_type, class GPUM>
82inline std::tuple<std::unique_ptr<GPUM>, std::unique_ptr<GPUM>>
83extractLowerAndUpperMatrices(const M& naturalMatrix, const std::vector<int>& reorderedToNatural)
84{
85 const size_t new_nnz = (naturalMatrix.nonzeroes() - naturalMatrix.N()) / 2;
86
87 M reorderedLower(naturalMatrix.N(), naturalMatrix.N(), new_nnz, M::row_wise);
88 M reorderedUpper(naturalMatrix.N(), naturalMatrix.N(), new_nnz, M::row_wise);
89
90 for (auto lowerIt = reorderedLower.createbegin(), upperIt = reorderedUpper.createbegin();
91 lowerIt != reorderedLower.createend();
92 ++lowerIt, ++upperIt) {
93
94 auto srcRow = naturalMatrix.begin() + reorderedToNatural[lowerIt.index()];
95
96 for (auto elem = srcRow->begin(); elem != srcRow->end(); ++elem) {
97 if (elem.index() < srcRow.index()) { // add index to lower matrix if under the diagonal
98 lowerIt.insert(elem.index());
99 } else if (elem.index() > srcRow.index()) { // add element to upper matrix if above the diagonal
100 upperIt.insert(elem.index());
101 }
102 }
103 }
104
105 return {std::unique_ptr<GPUM>(new auto(GPUM::fromMatrix(reorderedLower, true))),
106 std::unique_ptr<GPUM>(new auto(GPUM::fromMatrix(reorderedUpper, true)))};
107}
108} // namespace Opm::gpuistl::detail
109
110#endif
Contains wrappers to make the CuBLAS library behave as a modern C++ library with function overlading.
Definition autotuner.hpp:29
__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
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242