Direct steadystate solver

[labay11]

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.

Here you can find a working example for the steady-state of a "pumped cavity": https://gist.github.com/labay11/66934e15e648d06d8511fc557c259e59

[PhilipVinc]

https://docs.qojulia.org/api/#QuantumOptics.steadystate.iterative

probably you can pass a function as an iterative method even if that solver is not iterative