MPI/Distributed linear system solver?

New to Julia
,
21

Hi @Ribeiro!

I suggest the following code for your problem:

using MKL
using Krylov

n = 1000     # Dimension of your problem
T = Float64  # Precision of your problem
A = rand(T, n, n)
b = rand(T, n)

Id = Diagonal(ones(T,n))
F = lu(A)

# GMRES
solver_gmres = GmresSolver(n, n, 50, Vector{T})
solver_gmres.x .= 0.0
for i = 1 : 1000
  println(i)
  A .= A .+ Id  # replace by your update of A
  b .= 2 .* b   # replace by your update of b
  gmres!(solver_gmres, A, b, solver_gmres.x, N=F, atol=1e-12, rtol=0.0, ldiv=true)
  .... # Use solution solver_gmres.x, the solution of Ax = b
  if solver_gmres.stats.niter >= 50  # Tune this parameter if you recompute a preconditioner too often
    println("Recompute the preconditioner")
    F = lu(A)
  end
end

A few comments:

  • GMRES will reuse the same memory, so you don’t need to allocate storage for every call of gmres.
  • The previous solution is reused as an initial guess.
  • GMRES stops when ||r_k|| <= 1e-12.
  • The preconditioner is updated if it starts to be less efficient.
7 Likes