ezp
lightweight C++ wrapper for selected distributed solvers for linear systems
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 
47 namespace 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.