MKL 2025: Default iparm lead to bad results

[chnce]

I was chasing a bug in my software where the eigenvalues calculated for FEM matrices (using pypardiso as OPinv) were significantly off. With older MKL versions (2024) there were no issues.

Testing with mkl and pypardiso from conda-forge

mkl 2025.3.0 hac47afa_454 95.3 MiB conda https://conda.anaconda.org/conda-forge/
pypardiso 0.4.7 pyhd8ed1ab_0 15.4 KiB conda https://conda.anaconda.org/conda-forge/

Test script:

import numpy as np
import scipy.sparse as sp
import pypardiso


# FEM-like matrix, weakened diagonal (forces pivoting)
n = 100
e = np.ones(n)
em = np.ones(n - 1)

T = sp.diags([-em, 1e-3*e, -em], [-1, 0, 1], format="csc")
A = sp.kron(sp.eye(n, format="csc"), T) + sp.kron(T, sp.eye(n, format="csc"))

rng = np.random.default_rng(0)
x_true = rng.standard_normal(n * n)
b = A @ x_true

x = pypardiso.spsolve(A, b)

rel_res = np.linalg.norm(A @ x - b) / np.linalg.norm(b)
print("relative residual:", rel_res)

using the pypardiso default iparm leads to very large residuals.

setting iparm[0] = 1 seems to fix the issue.

The default iparm handling of MKL 2025 seems to have changed. Maybe pypardiso should initialise the iparm vector with the default values.
https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2025-2/pardiso-iparm-parameter.html