topic Re:Pardiso fails with high residual in Intel® oneAPI Math Kernel Library

Pardiso fails with high residual

Dan78 — Sat, 16 Dec 2023 12:21:30 GMT

Hi all,

I have a EM simulation tool that solved Maxwell equations in integral form. The matrices are ill-conditioned but sparse under specific circumstances. I have used MUMPS sparse direct solver and it worked very well but since we use Intel MKL for all other operations in our simulator, I decided to move to Pardiso. The problem is that when I use Pardiso the residual is very high and the solution is not correct at all. I have tried to re-produce the error using a small 6x6 sparse matrix which is attached. I have not been able to find the root cause but is there anyone that can help me with this issues? My original matrix is in COO format and I use MKL auxiliary routines to convert COO to CSR which is understandable by Pardiso. The coefficient matrix in my equation is unsymmetric and complex. I compile the code as below:

$ icpx pardiso_fail.cc -lmkl_core -lmkl_intel_ilp64 -lmkl_intel_thread -liomp5 -DMKL_ILP64

I have compiled the example that is with MKL and it works correctly with very low residual.

Any help is highly appreciated.

Regards,

Dan

Re: Pardiso fails with high residual

mecej4 — Sat, 16 Dec 2023 15:40:38 GMT

Something is wrong with your example matrix. You wrote that the matrix is 6 X 6, but you set n = 8 in the code.

Re: Pardiso fails with high residual

Dan78 — Sat, 16 Dec 2023 19:03:44 GMT

@mecej4Thanks for pointing this out. I had missed to update this from the example that is included with MKL. I updated the code to have correct values and the result is still unreliable. I have attached the corrected code and now the result is:

Reordering completed ...
Number of nonzeros in factors = 28
Number of factorization MFLOPS = 0
Factorization completed ...
Solving system with iparm[11] = 0
The solution of the system is:
6.01275e+06,6.01276e+06
6.01275e+06,6.01275e+06
6.01275e+06,6.01276e+06
6.01275e+06,6.01275e+06
1,1
1,1

Relative residual = 1.1547
Error: residual is too high!

Please find the correct code attached.

Regards,

Dan

Re: Pardiso fails with high residual

Dan78 — Sat, 16 Dec 2023 19:09:28 GMT

Sorry I could not attached the corrected source code. Here it is attached. I got error when I tried to attach the file so I changed the extension from ".cpp" to ".c". If you want to compile it please rename it to "pardiso_fail.cpp".

Thanks,

Dan

Re: Pardiso fails with high residual

Dan78 — Mon, 18 Dec 2023 17:32:27 GMT

Hi all,

Some updates here. I reduced the sparsity of my matrix which ended up in better condition number and then Pardiso started to give reliable results with low residual. I am suspecting that Pardiso is not able to solve equations when coefficient matrices are kind of ill-posed, highly indefinite and the condition number is high. Like I mentioned before, MUMPS has been able so solve such equations. I think Pardiso internally uses METIS for reordering like MUMPS does. MUMPS also supports Scotch which might perform better in some cases. I remember that I experienced similar problems in accuracy when I used Pardiso long time ago (Not MKL's version but the one that was implemented by Olaf Schenk). If any adjustment in the parameters would make the solution reliable, I would be more than happy to know about it.

Regards,

Dan

Re: Pardiso fails with high residual

mecej4 — Tue, 19 Dec 2023 14:19:10 GMT

I don't know what to make of your equations.

The first eq. is x_5 = 1 + i, the second is - x_5 = 1 + i. Similarly, the third and fourth equations are inconsistent. In other words, the matrix is rank-deficient, and the MKL Pardiso solver is not an appropriate tool.

Perhaps the small test problem matrix is not representative of the matrices that your large application generates. Otherwise, you have to rethink what "obtaining a solution" means, such as "least squares", i.e., norm(A.x - b), in the case of over-determined equation sets, or "least norm solution", i.e., norm(x) for under-determined equation sets. Pardiso (Intel's or Panua's) does not cover such problems.

Re: Pardiso fails with high residual

Dan78 — Tue, 19 Dec 2023 20:40:23 GMT

Hi @mecej4

Yes I get your point. This small matrix is probably ill-posed. I have some regularization techniques for such matrices but I didn't applied because I wanted to know whether Pardiso can be an option for us or not. I did some more tests and this time with larger matrices i.e., ~41000 x 41000 which are ~99% sparse. This time Pardiso failed like before with very large residual values but MUMPS solved exactly same equation with very good accuracy. I didn't calculate residual for MUMPS but when I plotted the results in MATLAB the results from Pardiso are very off while MUMPS are quite close to the solution for dense matrices using LU dense solver.

I use latest version of MUMPS i.e., 5.6.2 and compiled it with SCOTCH as reordering package and using ifx and icx compiler with MKL. I think Pardiso uses METIS for reordering.

Conclusion is that while Pardiso is a bit faster than MUMPS, but the results from MUMPS are very accurate. Again, I remember that I had done similar tests with Pardiso long time ago and got inaccurate results. I know the matrices that I use are ill-conditioned since they are coming from integral equations so this might be challenging for some solvers.

Regards,

Dan

Re:Pardiso fails with high residual

ShanmukhS_Intel — Wed, 27 Dec 2023 16:32:43 GMT

Hi Danesh,

Some updates here. I reduced the sparsity of my matrix which ended up in better condition number and then Pardiso started to give reliable results with low residual

>> Thanks for letting us know. We are looking into your issue. We will get back to you soon with an update.

Best Regards,

Shanmukh.SS

Re: Re:Pardiso fails with high residual

Dan78 — Wed, 27 Dec 2023 16:44:20 GMT

Hi @ShanmukhS_Intel,

Thanks for the update. Well, the challenge if to solve ill-conditioned coefficient matrices with a sparse direct solver with low residual. Otherwise, for a well-conditioned, or relatively ill-conditioned where you can find an efficient pre-conditioner, an iterative solver with lower cost could be used.

I have completely moved to MUMPS and now the solutions I get with exactly same matrices is very accurate. MUMPS calculates residual and backward error as a part of calculations, if the flag is set (Pardiso doesn't support this features I guess) and my residual and backward error is less than 1e-9. Also, when I plot the results, the difference is very low comparing to exact solution. MUMPS never failed even when I sparsified my matrix up to 99% and the error was still very low. I used SCOTCH for reordering and 3-steps if iterative refinement.

It would be good if MKL could also support MUMPS since this would be appropriate for the cases with large condition number of the cases that the equation is ill-posed.

Regards,

Dan