wrong solution for multiple rhs

[bdortdivanlioglu]

Hello,

I wanted to report that I am getting a solution with large errors from Pardiso. I tried to simplify the code in case you want to replicate the result. I have two rhs. If I solve each system separate, I get the correct solution for both.

[Amatrix.txt](https://github.com/JuliaSparse/Pardiso.jl/files/3811826/Amatrix.txt)
[xx.txt](https://github.com/JuliaSparse/Pardiso.jl/files/3811827/xx.txt)
[yy.txt](https://github.com/JuliaSparse/Pardiso.jl/files/3811828/yy.txt)

using SparseArrays, LinearAlgebra, Pardiso, DelimitedFiles

# load the matrices of size [A is symmetric and positive definite -> 59x59] [xx -> 59x2] [yy -> 59x2] 
Amatrix = readdlm("Amatrix.txt");  xx = readdlm("xx.txt"); yy = readdlm("yy.txt");

# convert to sparse 
A = sparse(Amatrix)

# solve (basic usage way)
ps = PardisoSolver()
solve!(ps, xx, A, yy)

@show norm(A*xx-yy) # what I get: 0.6085402893836485 

I'm using Ubuntu 18.04/Pardiso 6.0/gcc version 7.4.0

[KristofferC]

Not a problem with the MKL solver:

julia> ps = MKLPardisoSolver()

julia> @show norm(A*xx-yy) # what I get: 0.6085402893836485
norm(A * xx - yy) = 1.9314667586919627e-14
1.9314667586919627e-14

And it doesn't seem to happen with random matrices.

The next step would be to try this with some other Pardiso wrapper (or in C) and see if it can be reproduced. Otherwise, we need to find if we are passing something wrong to pardiso from the Julia side, but that would surprise me since we get correct results for a bunch of other cases.