ezp
lightweight C++ wrapper for selected distributed solvers for linear systems
Loading...
Searching...
No Matches
pardiso.hpp
Go to the documentation of this file.
1/*******************************************************************************
2 * Copyright (C) 2025-2026 Theodore Chang
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 ******************************************************************************/
38#ifndef PARDISO_HPP
39#define PARDISO_HPP
40
41#ifdef EZP_MKL
42
43#include "abstract/sparse_solver.hpp"
44
45#include <mpl/mpl.hpp>
46
47namespace ezp {
48 enum matrix_type : std::int8_t {
49 real_and_symmetric_positive_definite = 2,
50 real_and_symmetric_indefinite = -2,
51 real_and_structurally_symmetric = 1,
52 real_and_nonsymmetric = 11,
53 complex_and_hermitian_positive_definite = 4,
54 complex_and_hermitian_indefinite = -4,
55 complex_and_structurally_symmetric = 3,
56 complex_and_symmetric = 6,
57 complex_and_nonsymmetric = 13
58 };
59
60 enum message_level : std::int8_t {
61 no_output = 0,
62 print_statistical_information = 1
63 };
64
65 template<data_t DT, index_t IT> class pardiso final {
66 static constexpr IT one{1}, negone{-1}, PARDISO_ANA_FACT{12}, PARDISO_SOLVE{33}, PARDISO_RELEASE{-1};
67
68 std::int64_t pt[64]{};
69
70 IT iparm[64]{};
71
72 sparse_csr_mat<DT, IT> a_mat{};
73
74 const mpl::communicator& comm_world{mpl::environment::comm_world()};
75
76 const int comm = MPI_Comm_c2f(comm_world.native_handle());
77
78 const IT mtype, msglvl;
79
80 bool is_allocated{false};
81
82 auto alloc(sparse_csr_mat<DT, IT>&& A) {
83 dealloc();
84 is_allocated = true;
85 a_mat = std::move(A);
86 IT error{-1};
87 if constexpr(sizeof(IT) == 4) cluster_sparse_solver(pt, &one, &one, &mtype, &PARDISO_ANA_FACT, &a_mat.n, a_mat.data, a_mat.row_ptr, a_mat.col_idx, nullptr, &negone, iparm, &msglvl, nullptr, nullptr, &comm, &error);
88 else if constexpr(sizeof(IT) == 8) cluster_sparse_solver_64(pt, &one, &one, &mtype, &PARDISO_ANA_FACT, &a_mat.n, a_mat.data, a_mat.row_ptr, a_mat.col_idx, nullptr, &negone, iparm, &msglvl, nullptr, nullptr, &comm, &error);
89 return sync_error(error);
90 }
91
92 auto dealloc() {
93 if(!is_allocated) return;
94 is_allocated = false;
95 IT error;
96 if constexpr(sizeof(IT) == 4) cluster_sparse_solver(pt, &one, &one, &mtype, &PARDISO_RELEASE, &negone, nullptr, nullptr, nullptr, nullptr, &negone, iparm, &msglvl, nullptr, nullptr, &comm, &error);
97 else if constexpr(sizeof(IT) == 8) cluster_sparse_solver_64(pt, &one, &one, &mtype, &PARDISO_RELEASE, &negone, nullptr, nullptr, nullptr, nullptr, &negone, iparm, &msglvl, nullptr, nullptr, &comm, &error);
98 for(auto& i : pt) i = 0;
99 }
100
101 auto sync_error(IT error) {
102 comm_world.allreduce(mpl::min<IT>(), error);
103 return error;
104 }
105
106 public:
107 explicit pardiso(const IT mtype, const IT msglvl = 0)
108 : pardiso(matrix_type{static_cast<std::int8_t>(mtype)}, message_level{static_cast<std::int8_t>(msglvl)}) {}
109
110 explicit pardiso(const matrix_type mtype, const message_level msglvl = no_output)
111 : mtype(mtype)
112 , msglvl(msglvl) {
113 if constexpr(floating_t<DT>) {
114 if(mtype != real_and_symmetric_positive_definite && mtype != real_and_symmetric_indefinite && mtype != real_and_structurally_symmetric && mtype != real_and_nonsymmetric) throw std::invalid_argument("matrix type does not match the data type");
115 }
116 else {
117 if(mtype != complex_and_hermitian_positive_definite && mtype != complex_and_hermitian_indefinite && mtype != complex_and_structurally_symmetric && mtype != complex_and_symmetric && mtype != complex_and_nonsymmetric) throw std::invalid_argument("matrix type does not match the data type");
118 }
119 }
120
121 pardiso(const pardiso& other)
122 : a_mat(other.a_mat)
123 , mtype(other.mtype)
124 , msglvl(other.msglvl) { std::copy_n(other.iparm, 64, iparm); }
125 pardiso(pardiso&&) = delete;
126 pardiso& operator=(const pardiso&) = delete;
127 pardiso& operator=(pardiso&&) = delete;
128
129 ~pardiso() { dealloc(); }
130
131 auto& operator()(const IT index) { return iparm[index]; }
132
133 auto& iparm_default_value(const auto config) {
134 iparm[0] = config;
135 return *this;
136 };
137 auto& iparm_reducing_ordering(const auto config) {
138 iparm[1] = config;
139 return *this;
140 };
141 auto& iparm_user_permutation(const auto config) {
142 iparm[4] = config;
143 return *this;
144 };
145 auto& iparm_iterative_refinement(const auto config) {
146 iparm[7] = config;
147 return *this;
148 };
149 auto& iparm_pivoting_perturbation(const auto config) {
150 iparm[9] = config;
151 return *this;
152 };
153 auto& iparm_scaling(const auto config) {
154 iparm[10] = config;
155 return *this;
156 };
157 auto& iparm_transpose_matrix(const auto config) {
158 iparm[11] = config;
159 return *this;
160 };
161 auto& iparm_weighted_matching(const auto config) {
162 iparm[12] = config;
163 return *this;
164 };
165 auto& iparm_nnz_factor(const auto config) {
166 iparm[17] = config;
167 return *this;
168 };
169 auto& iparm_pivoting_type(const auto config) {
170 iparm[20] = config;
171 return *this;
172 };
173 auto& iparm_matrix_checker(const auto config) {
174 iparm[26] = config;
175 return *this;
176 };
177 auto& iparm_precision(const auto config) {
178 iparm[27] = config;
179 return *this;
180 };
181 auto& iparm_partial_solve(const auto config) {
182 iparm[30] = config;
183 return *this;
184 };
185 auto& iparm_zero_based_indexing(const auto config) {
186 iparm[34] = config;
187 return *this;
188 };
189 auto& iparm_schur_complement(const auto config) {
190 iparm[35] = config;
191 return *this;
192 };
193 auto& iparm_out_of_core(const auto config) {
194 iparm[59] = config;
195 return *this;
196 };
197
198 IT solve(sparse_csr_mat<DT, IT>&& A, full_mat<DT, IT>&& B) {
199 if(A.n != B.n_rows) return -1;
200
201 iparm[39] = 0; // force centralised input/output
202
203 if constexpr(std::is_same_v<DT, double> || std::is_same_v<DT, complex16>) iparm[27] = 0;
204 else if constexpr(std::is_same_v<DT, float> || std::is_same_v<DT, complex8>) iparm[27] = 1;
205
206 auto error = sync_error(alloc(std::move(A)));
207 if(error < 0) return error;
208
209 return solve(std::move(B));
210 }
211
212 IT solve(full_mat<DT, IT>&& B) {
213 if(a_mat.n != B.n_rows) return -1;
214
215 std::vector<DT> b_ref;
216 if(0 == comm_world.rank()) b_ref.resize(B.n_rows * B.n_cols);
217
218 void *b_ptr, *x_ptr;
219
220 if(0 == iparm[0] || 0 == iparm[5]) {
221 // b unchanged, x has solution
222 if(0 == comm_world.rank()) std::copy_n(B.data, b_ref.size(), b_ref.data());
223 b_ptr = b_ref.data();
224 x_ptr = B.data;
225 }
226 else {
227 // b has solution, x is still used
228 b_ptr = B.data;
229 x_ptr = b_ref.data();
230 }
231
232 IT error{-1};
233 if constexpr(sizeof(IT) == 4) cluster_sparse_solver(pt, &one, &one, &mtype, &PARDISO_SOLVE, &a_mat.n, a_mat.data, a_mat.row_ptr, a_mat.col_idx, nullptr, &B.n_cols, iparm, &msglvl, b_ptr, x_ptr, &comm, &error);
234 else if constexpr(sizeof(IT) == 8) cluster_sparse_solver_64(pt, &one, &one, &mtype, &PARDISO_SOLVE, &a_mat.n, a_mat.data, a_mat.row_ptr, a_mat.col_idx, nullptr, &B.n_cols, iparm, &msglvl, b_ptr, x_ptr, &comm, &error);
235
236 return sync_error(error);
237 }
238 };
239} // namespace ezp
240
241#endif
242
243#endif
244
pardiso
Solver for general sparse matrices.