MUMPS.jl icon indicating copy to clipboard operation
MUMPS.jl copied to clipboard

bug - solving for multiple right-hand-sides

Open hiemstar opened this issue 2 years ago • 3 comments

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.

hiemstar avatar Dec 20 '22 13:12 hiemstar

Thank you for the report. If I understand well, the issue is not with multiple right-hand sides, but rather with successive solves, where the factorization is reused?!

dpo avatar Dec 20 '22 15:12 dpo

Yes, that is correct.

On Tue, Dec 20, 2022 at 4:23 PM Dominique @.***> wrote:

Thank you for the report. If I understand well, the issue is not with multiple right-hand sides, but rather with successive solves, where the factorization is reused?!

— Reply to this email directly, view it on GitHub https://github.com/JuliaSmoothOptimizers/MUMPS.jl/issues/112#issuecomment-1359559147, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACKARLUPH5TWFUL672XGWCLWOHFQ3ANCNFSM6AAAAAATEQS4JE . You are receiving this because you authored the thread.Message ID: @.***>

hiemstar avatar Dec 20 '22 15:12 hiemstar

This can be solved by setting the phase to 4 at the end of each step with the function set_job!:

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

Output with Julia 1.9.1 on Ubuntu:

2.6134993330539642e-16
1.7513297027755214e-16
2.1165999061806354e-16
2.060185629889662e-16
1.2115267855126429e-16
1.843508341922951e-16
2.1968948894896043e-16
2.0330531384006366e-16
2.2474742915123943e-16
1.6937101875768352e-16

ovanvincq avatar Jul 03 '23 13:07 ovanvincq