I have written a FORTRAN program that uses Pardiso to solve some ill-conditioned system of equations. 95% of the time, the code works just fine. But for some parameters, it calculates completely wrong solutions. I could trace this back to a specific matrix that isn't solved correctly by pardiso. The matrix has a condition number of about 10^8, so in double precision it should definitely be possible to calculate some approximate solution.
I have attached a very short example program that simply loads this specific matrix and the right-hand side from a file (the file is "unformatted", it's only called "txt" because this forum doesn't like "dat" extensions) and solves the system with Dgesv from the MKL version of LAPACK and with Pardiso. It then calculates the norm of the residual for both solutions. For the solution computed by Dgesv, you get about 2*10^-7, which is perfectly ok for my needs. For Pardiso, the norm is about 4*10^+6, so the "solution" is complete nonsense.
I'd be very grateful if you could check if I'm doing something wrong here or if there really is a bug in Pardiso.
We're still using MKL 10.3 update 9 in my company, so if this is a bug, but it has already been fixed, I'd be happy to know it so I can ask our administrator to install the latest version as soon as possible.
Am I supposed to ask this question on Intel Premier Support? I don't really understand why I got no answer, even though it should be trivial to at least download the sample code, compile it with "ifort -mkl" and reproduce the problem (or tell me that the problem is not reproducible).
(I do not speak for Intel) I am able to duplicate the problem with MKL 10.3.6 and 11.0.2 on Windows. MKL 10.3.10 gives acceptable results, and the results agree with (i) the results given by the IMSL routine LSLXG for your matrix, and (ii) the results from Pardiso 4.1.2 (from the Pardiso Web site). The first ten elements of the solution vector are
[bash] 1 = 0.0000e+000 2 = 0.0000e+000 3 = -4.2686e+000 4 = -8.5265e+000 5 = -1.2763e+001 6 = -1.6968e+001 7 = -2.1131e+001 8 = -2.5241e+001 9 = -2.9288e+001 10 = -3.3263e+001 [/bash] MKL 10.3.10 gives the 2-norm of the residual as 1.107979e-019. I used the attached C program to read and process your data file. With the same program, MKL 11.0.2 gives the 2-norm of the residual as 2.642830e+005, and the first ten components of the solution as [bash] 1 = 0.0000e+000 2 = 0.0000e+000 3 = -4.2647e+000 4 = -8.5187e+000 5 = -1.2751e+001 6 = -1.6952e+001 7 = -2.1111e+001 8 = -2.5217e+001 9 = -2.9261e+001 10 = -3.3233e+001 [/bash] I do not think that there are any undefined variables at play here, since the IMSL results agree with those given by MKL 10.3.10, although the norm of the residuals is not quite the same (2.967E-10).
As to why nobody downloaded and ran your program, my guess is that doing so would help little, because the call to Pardiso has many arguments. In particular, one would have to check the values that you passed in IPARM against the Pardiso manual, because changing one or more of those values can significantly change the outcome of the calculations inside Pardiso.
I have reproduced your problem. Dgesv use global permutation for handling the zero pivots so it returns good result for matrix with many zero pivots. Otherwise pardiso doesn't use global permutation so it's strongly recommended to use additional iterative refinement algorithm during solving step. To achieve required result you may set iparm(8)=15 (iterative refinement).
Thanks for replying. Since a week had passed with no response from Intel, it appeared that the topic failed to arouse any interest!
I am afraid that the problem is not simply resolved by increasing the number of iterative refinement steps. Here are the results from various versions of MKL, obtained using unsym.c and lars.dat (attached to my previous post in this thread), with unsym.c modified to have iparm = 15 (zero-based indices used in C). For all the runs, I used the /Qopenmp compiler option and set OMP_NUM_THREADS=1.
MKL version 2-norm of residual
10.3.1 /IA32 8.120103e+000
11.0.2 /IA32 4.928226e+002
11.0.2 /x64 2.303653e-002
For the same matrix, Pardiso 4.1.2 gives the 2-norm of the residuals as 4.667119e-007 after 5 steps of refinement and 2.933695e-010 after 15 steps of refinement. Thus, a large number of refinement steps is not necessary in order to obtain a reasonably small residual norm.
The test matrix has a big conditional number (about 10^8) so the permutation on reordering step has a big influence for solution of the system. So you may try to use minimum degree algorithm on reordering step.
Which OS do you use linux or windows? It helps me to reproduse your results.
As I wrote in my first reply in this thread, the results that I posted were obtained on Windows 7/X64. Here are similar results obtained on Suse Linux 12.2, using unsym.c and lars.dat, compiled with -mkl -openmp, and run with OMP_NUM_THREADS=1:
MKL 11.0.1, IA32: 5.766764e-02 MKL 11.0.1, intel64: 3.102445e-02 Pardiso 4.1.2, intel64: 4.2221E-10
Although I accept your comments on pivoting and iterative refinement, from a user's point of view the results given by MKL 11.0.1 on Linux would be seen as inferior compared to those given by Pardiso 4.1.2 from Basel/Lugano.
As I wrote before it is strongly recommended to use minimum degree algorithm on reordering step and iterative refinement for this matrix.
I have tested this matrix on linux 64, windows 32/x64 and used MKL11.0.3, MKL10.3.12, MKL10.3.9. When I set iparm(2)=0 (minimum degree algorithm) and iparm(8)=15 I obtain the 2-norm of relative residual is equal to 2.33117335E-008 with all configurations.