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.