bug - solving for multiple right-hand-sides

[hiemstar]

When I solve for multiple right-hand-sides I get an incorrect result after the first solve. See the following MWE:

using MUMPS
using MPI
using LinearAlgebra
using SparseArrays

MPI.Init()

icntl = get_icntl(verbose=false)

A = sprand(10, 10, .2) +  sparse(I, 10, 10)
mumps = Mumps{Float64}(mumps_unsymmetric, icntl, default_cntl64)
associate_matrix!(mumps, A)
factorize!(mumps)

for _ = 1:10
  rhs = rand(10)
  associate_rhs!(mumps, rhs)
  solve!(mumps)
  x = get_solution(mumps)
  println(norm(x - A \ rhs) / norm(x))
end

finalize(mumps)
MPI.Finalize()

output:

julia> 

7.447297301440223e-17
0.4924343334935974
0.7160770210505035
0.6776858374598822
0.6331079647874924
0.6580384625500728
0.4572616048461576
0.545888578220492
0.7025421228833388
0.5773467516678198

julia> 

My julia version is 1.8.3 on Ubuntu.