Use LinearSolve's internal methods for preconditioners

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Change default solver detection in `eigensolve` when using `sigma` keyword argument (shift-inverse algorithm). If the operator is a `SparseMatrixCSC`, the default solver is `UMFPACKFactorization`, otherwise it is automatically chosen by LinearSolve.jl, depending on the type of the operator. ([#580]) 
 - Add keyword argument `assume_hermitian` to `liouvillian`. This allows users to disable the assumption that the Hamiltonian is Hermitian. ([#581])
+- Use LinearSolve's internal methods for preconditioners in `SteadyStateLinearSolver`. ([#588])
 
 ## [v0.38.1]
 Release date: 2025-10-27
@@ -362,3 +363,4 @@ Release date: 2024-11-13
 [#579]: https://github.com/qutip/QuantumToolbox.jl/issues/579
 [#580]: https://github.com/qutip/QuantumToolbox.jl/issues/580
 [#581]: https://github.com/qutip/QuantumToolbox.jl/issues/581
+[#588]: https://github.com/qutip/QuantumToolbox.jl/issues/588

docs/src/users_guide/steadystate.md

@@ -41,7 +41,7 @@ Although it is not obvious, the [`SteadyStateDirectSolver()`](@ref SteadyStateDi
 
 To overcome this, one can provide preconditioner that solves for an approximate inverse for the (modified) Liouvillian, thus better conditioning the problem, leading to faster convergence. The left and right preconditioner can be specified as the keyword argument `Pl` and `Pr`, respectively:
 ```julia
-solver = SteadyStateLinearSolver(alg = KrylovJL_GMRES(), Pl = Pl, Pr = Pr)
+solver = SteadyStateLinearSolver(alg = KrylovJL_GMRES())
 ```
 The use of a preconditioner can actually make these iterative methods faster than the other solution methods. The problem with precondioning is that it is only well defined for Hermitian matrices. Since the Liouvillian is non-Hermitian, the ability to find a good preconditioner is not guaranteed. Moreover, if a preconditioner is found, it is not guaranteed to have a good condition number. 
 

src/steadystate.jl

@@ -23,23 +23,20 @@ A solver which solves [`steadystate`](@ref) by finding the zero (or lowest) eige
 struct SteadyStateEigenSolver <: SteadyStateSolver end
 
 @doc raw"""
-    SteadyStateLinearSolver(alg = KrylovJL_GMRES(), Pl = nothing, Pr = nothing)
+    SteadyStateLinearSolver(
+        alg = KrylovJL_GMRES(; precs = (A, p) -> (ilu(A, τ = 0.01), I))
+    )
 
 A solver which solves [`steadystate`](@ref) by finding the inverse of Liouvillian [`SuperOperator`](@ref) using the `alg`orithms given in [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/).
 
 # Arguments
 - `alg::SciMLLinearSolveAlgorithm=KrylovJL_GMRES()`: algorithms given in [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/)
-- `Pl::Union{Function,Nothing}=nothing`: left preconditioner, see documentation [Solving for Steady-State Solutions](@ref doc:Solving-for-Steady-State-Solutions) for more details.
-- `Pr::Union{Function,Nothing}=nothing`: right preconditioner, see documentation [Solving for Steady-State Solutions](@ref doc:Solving-for-Steady-State-Solutions) for more details.
+
+# Note
+Refer to [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/) for more details about the available algorithms. For example, the preconditioners can be defined directly in the solver like: `SteadyStateLinearSolver(alg = KrylovJL_GMRES(; precs = (A, p) -> (I, Diagonal(A))))`.
 """
-Base.@kwdef struct SteadyStateLinearSolver{
-    MT<:Union{SciMLLinearSolveAlgorithm,Nothing},
-    PlT<:Union{Function,Nothing},
-    PrT<:Union{Function,Nothing},
-} <: SteadyStateSolver
-    alg::MT = KrylovJL_GMRES()
-    Pl::PlT = nothing
-    Pr::PrT = nothing
+Base.@kwdef struct SteadyStateLinearSolver{MT<:Union{SciMLLinearSolveAlgorithm,Nothing}} <: SteadyStateSolver
+    alg::MT = KrylovJL_GMRES(; precs = (A, p) -> (ilu(A, τ = 0.01), I))
 end
 
 @doc raw"""
@@ -159,14 +156,8 @@ function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateLinear
     L_tmp = L_tmp + Tn
 
     (haskey(kwargs, :Pl) || haskey(kwargs, :Pr)) && error("The use of preconditioners must be defined in the solver.")
-    if !isnothing(solver.Pl)
-        kwargs = merge((; kwargs...), (Pl = solver.Pl(L_tmp),))
-    elseif isa(L_tmp, SparseMatrixCSC)
-        kwargs = merge((; kwargs...), (Pl = ilu(L_tmp, τ = 0.01),))
-    end
-    !isnothing(solver.Pr) && (kwargs = merge((; kwargs...), (Pr = solver.Pr(L_tmp),)))
 
-    prob = LinearProblem(L_tmp, v0)
+    prob = LinearProblem{true}(L_tmp, v0)
     ρss_vec = solve(prob, solver.alg; kwargs...).u
 
     ρss = reshape(ρss_vec, N, N)
@@ -407,16 +398,10 @@ function _steadystate_fourier(
     v0[n_max*N+1] = weight
 
     (haskey(kwargs, :Pl) || haskey(kwargs, :Pr)) && error("The use of preconditioners must be defined in the solver.")
-    if !isnothing(solver.Pl)
-        kwargs = merge((; kwargs...), (Pl = solver.Pl(M),))
-    elseif isa(M, SparseMatrixCSC)
-        kwargs = merge((; kwargs...), (Pl = ilu(M, τ = 0.01),))
-    end
-    !isnothing(solver.Pr) && (kwargs = merge((; kwargs...), (Pr = solver.Pr(M),)))
     !haskey(kwargs, :abstol) && (kwargs = merge((; kwargs...), (abstol = tol,)))
     !haskey(kwargs, :reltol) && (kwargs = merge((; kwargs...), (reltol = tol,)))
 
-    prob = LinearProblem(M, v0)
+    prob = LinearProblem{true}(M, v0)
     ρtot = solve(prob, solver.alg; kwargs...).u
 
     offset1 = n_max * N