topic So far we haven't been able in Intel® oneAPI Math Kernel Library

MKL Pardiso fails to solve complex unsymmetric sparse matrix

areid2 — Wed, 09 Nov 2016 23:09:03 GMT

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?

Areid, this is also good

Gennady_F_Intel — Thu, 10 Nov 2016 02:16:15 GMT

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

Thanks for the help. I've

areid2 — Thu, 10 Nov 2016 05:32:00 GMT

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

Sorry, the attachment didn't

areid2 — Thu, 10 Nov 2016 05:34:03 GMT

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

Areid, Is it possible to send

Gennady_F_Intel — Thu, 10 Nov 2016 09:18:15 GMT

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

Hi,

areid2 — Tue, 22 Nov 2016 06:17:49 GMT

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

thanks Adam for the test. we

Gennady_F_Intel — Tue, 22 Nov 2016 14:39:18 GMT

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

i see with mkl 2017 u1

Gennady_F_Intel — Tue, 22 Nov 2016 15:31:41 GMT

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

areid2 — Tue, 22 Nov 2016 20:28:39 GMT

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.

Adam, could you please try to

Gennady_F_Intel — Wed, 23 Nov 2016 08:30:34 GMT

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

Thanks, I also get the

areid2 — Wed, 23 Nov 2016 08:45:36 GMT

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

yes, I agree with you, the

Gennady_F_Intel — Wed, 23 Nov 2016 08:54:10 GMT

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.

So far we haven't been able

areid2 — Thu, 24 Nov 2016 03:45:06 GMT

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.