Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

MKL Pardiso fails to solve complex unsymmetric sparse matrix

areid2
New Contributor I
769 Views

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?

0 Kudos
12 Replies
Gennady_F_Intel
Moderator
769 Views

Areid, this is also good place to report the issue!  Could you give us some more details to check the problem on our side? 

0 Kudos
areid2
New Contributor I
769 Views

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[1]=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.

Thanks,

Adam

0 Kudos
areid2
New Contributor I
769 Views

Sorry, the attachment didn't work with the last comment. It should be attached to this comment.

0 Kudos
Gennady_F_Intel
Moderator
769 Views

Areid, Is it possible to send us matrix without matlab call? We see matlab run but doesn’t see the way how you call pardiso to reproduce issue on our side

0 Kudos
areid2
New Contributor I
769 Views

Hi,

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.

Thanks,

Adam

 

 

0 Kudos
Gennady_F_Intel
Moderator
769 Views

thanks Adam for the test. we will check asap with the latest versions.

0 Kudos
Gennady_F_Intel
Moderator
769 Views

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?

0 Kudos
areid2
New Contributor I
769 Views

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.

0 Kudos
Gennady_F_Intel
Moderator
769 Views

Adam, could you please try to disable matching (ipatm[12]==0; With this setting, I see the max relative difference with correct solution: 9.9900102025820376e-010

0 Kudos
areid2
New Contributor I
769 Views

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?

Thanks,

Adam

0 Kudos
Gennady_F_Intel
Moderator
769 Views

yes, I agree with you, the matching is the default setting and it looks like this case is the real issue but we need to investigate this case. 

0 Kudos
areid2
New Contributor I
769 Views

So far we haven't been able to find any test cases that fail with iparm[12]=0, so it would appear that works around the problem and the bug must be in the weighted matching.

0 Kudos
Reply