Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

bug - solving for multiple right-hand-sides #112

Open
hiemstar opened this issue Dec 20, 2022 · 3 comments
Open

bug - solving for multiple right-hand-sides #112

hiemstar opened this issue Dec 20, 2022 · 3 comments

Comments

@hiemstar
Copy link

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.

@dpo
Copy link
Member

dpo commented Dec 20, 2022

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

@hiemstar
Copy link
Author

hiemstar commented Dec 20, 2022 via email

@ovanvincq
Copy link

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants