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!
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
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.
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.