I've run into a matrix where Pardiso is producing an incorrect result. I've checked that exact same matrix can be solved in other matrix solvers like Matlab and with both dense and sparse solvers. Pardiso doesn't return an error code, but the solution that it provides is clearly quite bad.
I have a Premier support account where I've tried to report this issue, but for some reason when I create a new case it doesn't show up under my account. I do get an issue reference number. Does anyone know if there is a phone number for Premier support? Any other ideas how I can report this issue?
Thanks for the help. I've attached a Matlab data file that contains the matrices A, b and Pardiso_x. Pardiso_x is the solution that Pardiso produces to the problem Ax=b. There's also a bit of Matlab code that computes b from A and Pardiso_x and compares that with the original b which shows it isn't a good solution.
We have tried to solve this problem with mtype = 13 and a sparse matrix input. We have tried the default initialization of iparam (using pardisoinit) as well as with iparam=0 (zero-based indexing on iparam) to workaround an issue in 11.3 update 3, as well as with 11.3 update 4.
You can also see that Matlab can easily invert this matrix as a dense matrix and find a good solution, so it doesn't seem to be an issue with the conditioning of the matrix. Pardiso also didn't produce an error when it computed this result.
If you aren't a Matlab user then let me know what formats you could support to check out this test case.
I've put together a self contained C++ example that reproduces the problem. We've tried on Linux and Windows with different versions of MKL and most configurations fail but in different ways. Some configurations fail with an error, some just return a bad result. There is a second copy of the matrix commented out in the file that is nearly identical with some differences in the 14th decimal place. Strange thing is that it often gives a different result even though it is only slightly different. In some configurations one of the 2 matrices will solve and the other won't!
I'd appreciate if anyone has some insight into this bug. We've tried lots of different cases with matrices similar to these and many of them give good solutions. It just seems that certain matrices fail to solve but it's not at all clear what is different about these matrices compared to the ones that work. It would appear pardiso is very sensitive to some very small differences in these cases. The correct answer that we're comparing with was computed by converting to a dense matrix and solving with Matlab.
i see with mkl 2017 u1 version of mkl on win64 system, the ****************** max relative difference with correct solution: 734.89125654341649
Do you see the similar difference?
Yes, I see almost the same value on mkl 2016 u2 under Linux too. When I found a configuration that gave the correct result, the number is more like 1e-5 which seems more reasonable.
The way that the relative difference is calculated means that the difference between the pardiso solution and the correct solution at some specific element in x is about 700x larger than the largest value in the correct solution. Since these are complex values we are using absolute values. So max relative error is max_element(abs(pardiso_x - correct_x))/max_element(abs(correct_x))
Thanks for looking into this.
Thanks, I also get the correct result with that setting. We'll try it on a few more of our test cases to see whether they all work with that setting. I'll let you know how it goes.
The weighted matching option appears to be the default for a nonsymmetric matrix and supposed to increase reliability. Do you think this is a bug, or does that technique not work well on some matrices?
So far we haven't been able to find any test cases that fail with iparm=0, so it would appear that works around the problem and the bug must be in the weighted matching.