py-mkl-pardiso

Python pybind11 wrapper for the Intel oneMKL PARDISO sparse direct solver.

pymklpardiso exposes double-precision real and complex support from Intel's PARDISO sparse direct solver to Python via pybind11. It works with SciPy sparse matrices in CSR format and NumPy arrays.

Supported platforms: Linux (x86_64), Windows (AMD64).

Installation

pip install py-mkl-pardiso

Or install from source (requires MKL):

git clone https://github.com/bodono/py-mkl-pardiso.git
cd py-mkl-pardiso
pip install -e ".[test]"

Building from source requires:

• Intel oneMKL (set MKLROOT if not auto-detected)
• A C++17 compiler
• Python >= 3.10
• pybind11 >= 2.12

Quick start

import numpy as np
import scipy.sparse as sp
from pymklpardiso import PardisoSolver, MTYPE_REAL_SYM_POSDEF

# Build a symmetric positive-definite matrix (upper triangle, CSR)
A_full = np.array([
    [4.0, 1.0],
    [1.0, 3.0],
])
A_upper = sp.csr_matrix(np.triu(A_full))
A_upper.sort_indices()

# Create solver — analyzes, factors, and is ready to solve
solver = PardisoSolver(A_upper, MTYPE_REAL_SYM_POSDEF)

b = np.array([1.0, 2.0])
x = solver.solve(b)
print(x)  # [0.09090909 0.63636364]

Complex matrices use the corresponding complex matrix type and NumPy complex128 values. Hermitian and complex-symmetric matrices are supplied as an upper-triangular CSR matrix, just like real symmetric matrices:

from pymklpardiso import MTYPE_COMPLEX_HERM_POSDEF

A_full = np.array([
    [4.0, 1.0 + 1.0j],
    [1.0 - 1.0j, 3.0],
])
A_upper = sp.csr_matrix(np.triu(A_full))

solver = PardisoSolver(A_upper, MTYPE_COMPLEX_HERM_POSDEF)
b = np.array([1.0 + 2.0j, 3.0 - 1.0j])
x = solver.solve(b)

Hermitian matrix diagonals must be real. The wrapper validates this whenever numeric values are loaded, including during construction, factor(), and refactor().

Refactoring workflow

When the sparsity pattern stays the same but values change (e.g., in an iterative algorithm), use refactor() to skip symbolic analysis:

solver = PardisoSolver(A_upper, MTYPE_REAL_SYM_POSDEF)
for new_values in value_generator:
    solver.refactor(new_values)
    x = solver.solve(b)

API reference

PardisoSolver(A, mtype, iparms=None, msglvl=0)

Create a PARDISO solver instance. The constructor extracts the CSR sparsity pattern from A, applies any iparms overrides, and runs symbolic analysis

• numeric factorization so the solver is ready to call solve().

Parameter	Type	Default	Description
A	sparse CSR	(required)	Square sparse matrix (any object with indptr, indices, data, shape). For symmetric and Hermitian types, pass only the upper triangle in CSR format.
mtype	int	(required)	Matrix type (see constants below).
iparms	dict	None	Optional {index: value} iparm overrides.
msglvl	int	0	Message level (0 = silent, 1 = print statistics).

Matrix type constants

Constant	Value	Description
MTYPE_REAL_STRUCT_SYM	1	Real structurally symmetric
MTYPE_REAL_SYM_POSDEF	2	Real symmetric positive definite
MTYPE_REAL_SYM_INDEF	-2	Real symmetric indefinite
MTYPE_REAL_NONSYM	11	Real nonsymmetric
MTYPE_COMPLEX_STRUCT_SYM	3	Complex structurally symmetric
MTYPE_COMPLEX_HERM_POSDEF	4	Complex Hermitian positive definite
MTYPE_COMPLEX_HERM_INDEF	-4	Complex Hermitian indefinite
MTYPE_COMPLEX_SYM	6	Complex symmetric
MTYPE_COMPLEX_NONSYM	13	Complex nonsymmetric

Core methods

solver.solve(b) Solve Ax = b. Accepts 1D (n,) or 2D (n, nrhs) arrays. Returns the solution as a new NumPy array (Fortran-contiguous for 2D). Real solvers return float64; complex solvers return complex128 and promote real right-hand sides automatically.

solver.solve_into(b, x) Solve Ax = b writing into pre-allocated x. For 2D arrays, both b and x must be Fortran-contiguous. The output must have dtype float64 for a real solver or complex128 for a complex solver.

solver.refactor(values) Re-factorize with new nonzero values (phase 22 only). Does not re-run symbolic analysis. Raises if symbolic analysis is invalid; use factor() to re-analyze from scratch. values must match the stored sparsity pattern exactly.

solver.factor(values) Re-analyze and re-factorize with new values (phases 11 + 22). This always runs fresh symbolic analysis. Use this for error recovery or when iparm changes require fresh symbolic analysis. values must match the stored sparsity pattern exactly.

Analysis invalidation: symbolic analysis (phase 11) becomes invalid when an iparm value changes via set_iparm() or set_iparm_all(), or after release() / reset(). Changing numeric values alone does not invalidate analysis. After a successful factor(), subsequent refactor() calls continue to work, including with value-dependent analysis settings such as iparm[10] = 1. When analysis is invalid, refactor() raises; call factor() to re-analyze and recover.

Other methods

Method	Description
solver.release()	Free PARDISO internal memory.
solver.n	Matrix dimension (property).
solver.nnz	Number of nonzeros (property).
solver.mtype	Matrix type (property).
solver.set_perm(perm)	Set fill-reducing permutation.
solver.clear_perm()	Clear permutation.
solver.has_perm()	Whether a permutation is set.
solver.set_iparm(idx, value)	Set a single iparm entry.
solver.get_iparm()	Get all 64 iparm values.
solver.get_iparm_value(idx)	Get a single iparm value.
solver.set_iparm_all(iparm)	Set all 64 iparm values.
solver.set_msglvl(msglvl)	Change message level.
solver.run_phase(phase)	Run an arbitrary PARDISO phase.
solver.run_phase_into(phase, b, x)	Run a phase with RHS/output arrays.

iparm notes

• iparm[0] is locked to 1 (user-supplied parameters).
• iparm[27] is locked to 0 (double precision). The wrapper supports float64 and complex128, not PARDISO's single-precision mode.
• iparm[34] is locked to 1 (zero-based indexing).
• See the MKL PARDISO iparm documentation for all parameters.

License

MIT