From fa4304f945118ba9f52321576f8df48a38dbe1f0 Mon Sep 17 00:00:00 2001 From: capric98 <42015599+capric98@users.noreply.github.com> Date: Wed, 23 Jul 2025 20:01:09 -0500 Subject: [PATCH 1/3] manually set iparm --- pypardiso/pardiso_wrapper.py | 28 +++++++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/pypardiso/pardiso_wrapper.py b/pypardiso/pardiso_wrapper.py index be9174b..4ee9190 100644 --- a/pypardiso/pardiso_wrapper.py +++ b/pypardiso/pardiso_wrapper.py @@ -142,6 +142,24 @@ def __init__(self, mtype=11, phase=13, size_limit_storage=5e7): self.size_limit_storage = size_limit_storage self._solve_transposed = False + self._init_iparm() # manually setup iparm + + def _init_iparm(self): + # https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2025-2/pardiso-iparm-parameter.html + self.iparm[0 ] = 1 # manually setting + self.iparm[1 ] = 103 # OMP version of fill-in reduce, maximum level of optimization (L=10) + self.iparm[5 ] = 0 # set to 1 if overwrite_b + self.iparm[7 ] = 2 # default 2-steps refinement + self.iparm[9 ] = 13 # set to 8 if symmetric + self.iparm[10] = 1 # set to 0 if symmetric + self.iparm[11] = 0 # Ax=b, no transpose + self.iparm[12] = 1 # enable matching for nonsymmetric + self.iparm[17] = 0 # disable nnz report + self.iparm[20] = 1 # default pivoting + self.iparm[23] = 0 # cannot use 10, idk why + self.iparm[34] = 1 # zero-based indexing + self.iparm[36] = 0 # CSR format + def factorize(self, A): """ Factorize the matrix A, the factorization will automatically be used if the same matrix A is passed to the @@ -270,9 +288,13 @@ def _call_pardiso(self, A, b): c_int32_p = ctypes.POINTER(ctypes.c_int32) c_float64_p = ctypes.POINTER(ctypes.c_double) - # 1-based indexing - ia = A.indptr + 1 - ja = A.indices + 1 + ia = A.indptr + ja = A.indices + + if self.iparm[34] == 0: + # 1-based indexing + ia = ia + 1 + ja = ja + 1 self._mkl_pardiso(self.pt.ctypes.data_as(ctypes.POINTER(self._pt_type[0])), # pt ctypes.byref(ctypes.c_int32(1)), # maxfct From 87995d5ad4dca509a6bd4381bd8de596252d7d9c Mon Sep 17 00:00:00 2001 From: capric98 <42015599+capric98@users.noreply.github.com> Date: Wed, 23 Jul 2025 20:15:00 -0500 Subject: [PATCH 2/3] consist with fortran style --- pypardiso/pardiso_wrapper.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/pypardiso/pardiso_wrapper.py b/pypardiso/pardiso_wrapper.py index 4ee9190..0a782c3 100644 --- a/pypardiso/pardiso_wrapper.py +++ b/pypardiso/pardiso_wrapper.py @@ -145,20 +145,20 @@ def __init__(self, mtype=11, phase=13, size_limit_storage=5e7): self._init_iparm() # manually setup iparm def _init_iparm(self): - # https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2025-2/pardiso-iparm-parameter.html - self.iparm[0 ] = 1 # manually setting - self.iparm[1 ] = 103 # OMP version of fill-in reduce, maximum level of optimization (L=10) - self.iparm[5 ] = 0 # set to 1 if overwrite_b - self.iparm[7 ] = 2 # default 2-steps refinement - self.iparm[9 ] = 13 # set to 8 if symmetric - self.iparm[10] = 1 # set to 0 if symmetric - self.iparm[11] = 0 # Ax=b, no transpose - self.iparm[12] = 1 # enable matching for nonsymmetric - self.iparm[17] = 0 # disable nnz report - self.iparm[20] = 1 # default pivoting - self.iparm[23] = 0 # cannot use 10, idk why - self.iparm[34] = 1 # zero-based indexing - self.iparm[36] = 0 # CSR format + # https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2025-2/pardiso.html + self.set_iparm(1 , 1 ) # manually setting + self.set_iparm(2 , 103) # OMP version of fill-in reduce, maximum level of optimization (L=10) + self.set_iparm(6 , 0 ) # set to 1 if overwrite_b + self.set_iparm(8 , 2 ) # default 2-steps refinement + self.set_iparm(10, 13 ) # set to 8 if symmetric + self.set_iparm(11, 1 ) # set to 0 if symmetric + self.set_iparm(12, 0 ) # Ax=b, no transpose + self.set_iparm(13, 1 ) # enable matching for nonsymmetric + self.set_iparm(18, 0 ) # disable nnz report + self.set_iparm(21, 1 ) # default pivoting + self.set_iparm(24, 0 ) # cannot use 10, idk why + self.set_iparm(35, 1 ) # zero-based indexing + self.set_iparm(37, 0 ) # CSR format def factorize(self, A): """ @@ -291,7 +291,7 @@ def _call_pardiso(self, A, b): ia = A.indptr ja = A.indices - if self.iparm[34] == 0: + if self.get_iparm(35) == 0: # 1-based indexing ia = ia + 1 ja = ja + 1 From 0ca605a35f1fe5bb4faa671d7d6ca454c708b800 Mon Sep 17 00:00:00 2001 From: capric98 <42015599+capric98@users.noreply.github.com> Date: Wed, 23 Jul 2025 20:44:49 -0500 Subject: [PATCH 3/3] support overwrite_b --- pypardiso/pardiso_wrapper.py | 26 +++++++++++++++++++------- pypardiso/scipy_aliases.py | 4 ++-- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/pypardiso/pardiso_wrapper.py b/pypardiso/pardiso_wrapper.py index 0a782c3..a9529c1 100644 --- a/pypardiso/pardiso_wrapper.py +++ b/pypardiso/pardiso_wrapper.py @@ -145,7 +145,7 @@ def __init__(self, mtype=11, phase=13, size_limit_storage=5e7): self._init_iparm() # manually setup iparm def _init_iparm(self): - # https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2025-2/pardiso.html + # https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2025-2/pardiso-iparm-parameter.html self.set_iparm(1 , 1 ) # manually setting self.set_iparm(2 , 103) # OMP version of fill-in reduce, maximum level of optimization (L=10) self.set_iparm(6 , 0 ) # set to 1 if overwrite_b @@ -181,7 +181,7 @@ def factorize(self, A): b = np.zeros((A.shape[0], 1)) self._call_pardiso(A, b) - def solve(self, A, b): + def solve(self, A, b, overwrite_b: bool=False): """ solve Ax=b for x @@ -196,14 +196,24 @@ def solve(self, A, b): """ self._check_A(A) - b = self._check_b(A, b) + b_fcontig = self._check_b(A, b) + + iparm6 = self.get_iparm(6) + if not(b_fcontig is b): + if overwrite_b: + warnings.warn('overwrite_b takes no effect: input b is not f_contiguous or the dtype is not compatible', + SparseEfficiencyWarning) + overwrite_b = True # already made a copy of b, safe to overwrite + + if overwrite_b: + self.set_iparm(6, 1) # overwrite if self._is_already_factorized(A): self.set_phase(33) else: self.set_phase(13) - x = self._call_pardiso(A, b) + x = self._call_pardiso(A, b_fcontig, overwrite_b=overwrite_b) # it is possible to call the solver with empty columns, but computationally expensive to check this # beforehand, therefore only the result is checked for infinite elements. @@ -211,6 +221,7 @@ def solve(self, A, b): # warnings.warn('The result contains infinite elements. Make sure that matrix A contains no empty columns.', # PyPardisoWarning) # --> this check doesn't work consistently, maybe add an advanced input check method for A + self.set_iparm(6, iparm6) # revert return x @@ -281,9 +292,10 @@ def _check_b(self, A, b): return b - def _call_pardiso(self, A, b): + def _call_pardiso(self, A, b, overwrite_b: bool=False): - x = np.zeros_like(b) + # PARDISO will fill x with 0 if iparm(6)=1 + x = np.empty_like(b) if overwrite_b else np.zeros_like(b) pardiso_error = ctypes.c_int32(0) c_int32_p = ctypes.POINTER(ctypes.c_int32) c_float64_p = ctypes.POINTER(ctypes.c_double) @@ -316,7 +328,7 @@ def _call_pardiso(self, A, b): if pardiso_error.value != 0: raise PyPardisoError(pardiso_error.value) else: - return np.ascontiguousarray(x) # change memory-layout back from fortran to c order + return b if overwrite_b else x # no need to force c_contiguous return def get_iparms(self): """Returns a dictionary of iparms""" diff --git a/pypardiso/scipy_aliases.py b/pypardiso/scipy_aliases.py index 7271c39..d48bb8e 100644 --- a/pypardiso/scipy_aliases.py +++ b/pypardiso/scipy_aliases.py @@ -9,7 +9,7 @@ pypardiso_solver = PyPardisoSolver() -def spsolve(A, b, factorize=True, squeeze=True, solver=pypardiso_solver, *args, **kwargs): +def spsolve(A, b, factorize=True, squeeze=True, solver=pypardiso_solver, overwrite_b: bool=False, *args, **kwargs): """ This function mimics scipy.sparse.linalg.spsolve, but uses the Pardiso solver instead of SuperLU/UMFPACK @@ -45,7 +45,7 @@ def spsolve(A, b, factorize=True, squeeze=True, solver=pypardiso_solver, *args, if factorize and not solver._is_already_factorized(A): solver.factorize(A) - x = solver.solve(A, b) + x = solver.solve(A, b, overwrite_b=overwrite_b) if squeeze: return x.squeeze() # scipy spsolve always returns vectors with shape (n,) indstead of (n,1)