I think it would be useful to include a direct steady-state solver which computes the solution of the system A x = 0 in an exact way, just like QuTip steady-state method.
For instance, using Pardiso.jl, the steady-state can be found as
function steadystate_mkl(L::SparseSuperOpType)
n, m = length(L.basis_r[1]), length(L.basis_r[2])
ps = MKLPardisoSolver()
b = zeros(ComplexF64, n*m)
# unity trace condition for the linear solvers
w = mean(abs.(L.data))
b[1] = w
A = L.data + sparse(ones(n), 1:(n+1):n^2, fill(w, n), n*m, n*m)
x = solve(ps, A, b)
data = reshape(x, n, m)
DenseOperator(L.basis_r[1], L.basis_r[2], data)
end
which is just slightly slower than the iterative method but requires fewer allocations and is more exact.
I think it would be useful to include a direct steady-state solver which computes the solution of the system
A x = 0
in an exact way, just likeQuTip
steady-state method.For instance, using Pardiso.jl, the steady-state can be found as
which is just slightly slower than the iterative method but requires fewer allocations and is more exact.
Here you can find a working example for the steady-state of a "pumped cavity": https://gist.github.com/labay11/66934e15e648d06d8511fc557c259e59