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().

Matrix type constants

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

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