search

JuliaLinearAlgebra / IterativeSolvers.jl

Iterative algorithms for solving linear systems, eigensystems, and singular value problems
MIT License
403 stars 106 forks source link

GMRES throws error with left preconditioner #312

Closed alec-hoyland closed 2 years ago

alec-hoyland commented 2 years ago

This is a method error due to a view on A. I am working on a bugfix.

julia> using LinearAlgebra, IterativeSolvers

julia> A = rand(10,10); b = rand(10); M = diagm(diag(A));

julia> gmres(A, b);

julia> gmres(A, b; Pl=M);
ERROR: MethodError: no method matching ldiv!(::Matrix{Float64}, ::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})
Closest candidates are:
  ldiv!(::Any, ::Identity, ::Any) at ~/.julia/dev/IterativeSolvers/src/common.jl:32
  ldiv!(::StridedVecOrMat{T}, ::SuiteSparse.UMFPACK.UmfpackLU{T}, ::StridedVecOrMat{T}) where T<:Union{Float64, ComplexF64} at ~/Julia/share/julia/stdlib/v1.7/SuiteSparse/src/umfpack.jl:764
  ldiv!(::StridedVecOrMat{T}, ::Transpose{T, <:SuiteSparse.UMFPACK.UmfpackLU{T}}, ::StridedVecOrMat{T}) where T<:Union{Float64, ComplexF64} at ~/Julia/share/julia/stdlib/v1.7/SuiteSparse/src/umfpack.jl:766
  ...
Stacktrace:
 [1] init!(arnoldi::IterativeSolvers.ArnoldiDecomp{Float64, Matrix{Float64}}, x::Vector{Float64}, b::Vector{Float64}, Pl::Matrix{Float64}, Ax::Vector{Float64}; initially_zero::Bool)
   @ IterativeSolvers ~/.julia/dev/IterativeSolvers/src/gmres.jl:249
 [2] gmres_iterable!(x::Vector{Float64}, A::Matrix{Float64}, b::Vector{Float64}; Pl::Matrix{Float64}, Pr::Identity, abstol::Float64, reltol::Float64, restart::Int64, maxiter::Int64, initially_zero::Bool, orth_meth::ModifiedGramSchmidt)
   @ IterativeSolvers ~/.julia/dev/IterativeSolvers/src/gmres.jl:126
 [3] gmres!(x::Vector{Float64}, A::Matrix{Float64}, b::Vector{Float64}; Pl::Matrix{Float64}, Pr::Identity, abstol::Float64, reltol::Float64, restart::Int64, maxiter::Int64, log::Bool, initially_zero::Bool, verbose::Bool, orth_meth::ModifiedGramSchmidt)
   @ IterativeSolvers ~/.julia/dev/IterativeSolvers/src/gmres.jl:200
 [4] gmres(A::Matrix{Float64}, b::Vector{Float64}; kwargs::Base.Pairs{Symbol, Matrix{Float64}, Tuple{Symbol}, NamedTuple{(:Pl,), Tuple{Matrix{Float64}}}})
   @ IterativeSolvers ~/.julia/dev/IterativeSolvers/src/gmres.jl:143
 [5] top-level scope
   @ REPL[4]:1
alec-hoyland commented 2 years ago

Arbitrary preconditioners are hard to support. Better to add a preconditioner to Preconditioners.jl than to automatically factorize since you run into type piracy issues.