GESV and PARDISO give different solutions

[Alex2]

Hello!

I use

Windows 10

Intel(R) Core(TM) i5-3320M

Intel(R) Visual Fortran Compiler Professional Edition 11.1.072 Update 9 for Windows*

Intel(R) Math Kernel Library 10.2 Update 7 for Windows* OS

I was solving a functional equation, involving a repeated solution of a big system of linear equations. GESV does the job reliably well. The system is sparse with a very small proportion of non-zero elements, therefore I wanted to try the sparse solver, expecting it would work much more efficiently. This was my first experience using PARDISO, I cannot be sure that I set up everything correctly. It is indeed much more efficient, but gives an incorrect solution.

I tested the methods on a simpler problem, the coefficient matrix from the MKL manual, and both solvers give the same solution. It is a bit confusing. I did not dig into computational details of PARDISO, I chose default values for iparm.

I attach a zip file with Visual Studio project. To illustrate the case, I give an example of an incorrect solution from PARDISO in the simplest case that I could make (14 equations), the real problem that I was solving had 5,000 equations (but can have much more, depends on the discretization of the domain). In this example, there is a significant difference in the seventh element of solution vectors. The code also has a commented "Alternative simple example", based on the coefficient matrix from the MKL manual, that gives the same solutions from both methods.

Thank you in advance for help.

[Gennady_F_Intel]

What is the reason why you use 6 years old version of mkl?  You could try the current version 2018 u3 and check the results with this version.

I hope you don't expect to see exactly the same results from gesv and pardiso!

[mecej4]

If I recall correctly, with version 11 of the compiler package, you had to build by yourself the blas95 and lapack95 libraries, modules and DLLs from the sources and makefile that were supplied by Intel. Current releases come with those libraries, modules and DLLs already built by Intel.

I ran your code with 18.0.3, and obtained ||Ax - b|| < 1E-13 from both GESV and Pardiso. The two columns of results printed by your program agreed to the full 15 digits in most cases; where they differed, it was in the last digit.

 GESV
 ----
 ||Ax-b|| =   7.049916206369744E-014


 PARDISO
 -------
 Number of non-zero elements:  40 out of  196
 error =            0
 ||Ax-b|| =   5.084821452783217E-014


 Solution vectors from GESV and PARDISO
 --------------------------------------
  -462.805009448382       -462.805009448382
  -458.077216751934       -458.077216751934
  -457.848478087207       -457.848478087207
  -457.714674539041       -457.714674539042
  -457.619739410817       -457.619739410817
  -457.546102005193       -457.546102005193
  -457.485935864588       -457.485935864588
  -449.038739242094       -449.038739242094
  -444.310959586318       -444.310959586318
  -444.082220945388       -444.082220945388
  -443.948417408818       -443.948417408819
  -443.853482287929       -443.853482287929
  -443.779844887535       -443.779844887535
  -443.719678750926       -443.719678750926

[Alex2]

Gennady F. (Intel) wrote:

This is a fair point, and indeed I planned to try a newer version later. Other than this issue, the old version works great. Once I needed to compile on Windows a complicated C project made on a Unix system that uses MKL and GSL libraries. I could accomplish that using Cygwin and an older MKL version, while the updated version did not succeed.

I do not expect the exact same numbers of course, but do expect both methods to solve the system up to a reasonable error tolerance. What I found that GESV solves the system, but the "solution" vector from PARDISO is not a solution.

mecej4 wrote:

Thank you very much for running the code using a never version! I see that your output is exactly the same (all the digits) as I get, apart from the seventh element of the PARDISO solution. In my case it is -443.857328237060, and in your case it is -457.485935864588. So the problem is not present.

It is still interesting what is the reason for the discrepancy in the older version, but as Gennady suggested, it would be good indeed to transition to a newer version.