Pardiso regression in 2025.1.0

[urob]

After upgrading MKL from 2024.2.2 to 2025.1.0 the Pardiso solver produces incorrect results for a variety of problems.

I am attaching a reproducer matrix A.txt and vector b.txt. Solving the system with default iparms in 2024.2.2 produces a solution x such that Ax - b is indeed numerically zero.

Solving the same system in 2025.1.0 produces solutions for which Ax - b is clearly non-zero (with values ranging into the thousands). This is true for the default iparms and any alternative I tried.

[c_sim]

Hi,

Thank you for reporting this bug. We have successfully reproduced it and are currently working to identify the root cause. In the meantime, we have two workarounds:

1) Use the new factorization algorithm, i.e.,  set iparm[23]=1. Additionally, ensure that matching and scaling are disabled by setting iparm[10]=0 and iparm[12]=0.

2) Use default factorization algorithm, i.e., set iparm[23]=0, and execute phase 1 in a serial manner using only one thread.

Attached, you will find a reproducer and the two workarounds mentioned above. Please let us know if these solutions are effective for you.

We will provide updates as soon as we identify the issue and have a confirmed date for the bug fix.

Thank you once again for reporting.

Kind Regards,

Chris

[urob]

Thanks Chris for reproducing and providing the workarounds! I can confirm that both workarounds are effective.  

One minor observation about workaround 1: While the solutions are numerically correct, they are a few orders of magnitude less accurate than the default options under 2024.2.2 (residuals of the order 1e-8 as opposed to 1e-19). Setting iparms[12] = 1 as suggested in the code comment does not help, because that induces Pardiso to use the classic factorization algorithm regardless of the value chosen for iparms[23].

Thanks again for looking into this

[c_sim]

Hi,

We have identified the issue, and a fix will be included in the next release of Intel oneMKL.

To answer your previous question: As you noted, enabling scaling or matching (iparm[10]=1 and/or iparm[12]=1) will make PARDISO to select the classic factorization path, where this bug is present. Therefore, if you need to enable these features, currently only workaround 2 would be the option.

Thank you for bringing this up.

Kind Regards,

Chris

[urob]

Thanks for the update, that's great to hear! Looking forward to the next release.