Request for Steady State Sparse GPU Solving

[MatthewBelzer]

Problem Description

I know there was some discussion on this in the past, but I was wondering if there are any plans to implement sparse GPU solving for steady state problems.

I ask because I am currently working on an optimization problem, where my objective function is significantly bottlenecked by solving the steady state of a Lindbladian. My optimization parameters are the parameters of my open quantum system's Hamiltonian and collapse operators, so the steady state density matrix has to be recalculated each iteration. Improving sparse GPU efficiency may help improve my optimization's convergence speed.

Proposed Solution

I am not sure if there is a way to fix this.

Alternate Solutions

No response

Additional Context

No response

[albertomercurio]

Hello,

in principle steadystate already support GPU objects if using the SteadyStateLinearSolver

using QuantumToolbox
using BenchmarkTools
using CUDA
CUDA.allowscalar(false)

# %%

F = 0.2
γ = 0.1
nth = 10

a = destroy(400)

H = a' * a + F * (a * a')
c_ops = [sqrt((1 + nth) * γ) * a, sqrt(nth * γ) * a']

L = liouvillian(H, c_ops)

ρss = steadystate(L)

# %%

L_gpu = cu(L)

ρss_gpu = steadystate(L_gpu, solver=SteadyStateLinearSolver())

It uses LinearSolve.jl with the KrylovJL_GMRES() algorithm as default. This calls Krylov.jl routines, which supports very well GPU, as they explain here.

However, the GPU support in steadystate can be still improved. Indeed, the CPU version introduces a preconditioner that significantly accelerates the convergence, while the GPU case does not. This means that the CPU case might still be faster than the GPU case.

As suggested by Krylov.jl, one may use some preconditioners of KrylovPreconditioners.jl, but I have never tested them.

Anyhow, the SteadyStateLinearSolver should already make the steady state calculation faster for large Liouvillians on CPU.

[albertomercurio]

Ok, with the latest version of QuantumToolbox.jl you can do something like

using QuantumToolbox
using KrylovPreconditioners
using CUDA
CUDA.allowscalar(false)

# %%

F = 0.5
γ = 0.1
nth = 15

a = destroy(500)

H = a' * a + F * (a * a')
c_ops = [sqrt((1 + nth) * γ) * a, sqrt(nth * γ) * a']

L = liouvillian(H, c_ops)

solver = SteadyStateLinearSolver()
ρss = steadystate(L, solver=solver)

# %%

L_gpu = CUSPARSE.CuSparseMatrixCSR(L)

solver = SteadyStateLinearSolver(Pr=KrylovPreconditioners.kp_ilu0)
ρss_gpu = steadystate(L_gpu, solver=solver)

Where I use KrylovPreconditioners.jl to also use some GPU-compatible preconditioner.

I hope this comment should fix you problem. If so, we can close this issue.

[ytdHuang]

@MatthewBelzer

[MatthewBelzer]

Hi, thank you for your responses. The first script you sent runs fine on my computer. However, when I run the the sparse gpu array code, it gives me this scalar indexing error:

ERROR: Scalar indexing is disallowed.
Invocation of getindex resulted in scalar indexing of a GPU array.
This is typically caused by calling an iterating implementation of a method.
Such implementations *do not* execute on the GPU, but very slowly on the CPU,
and therefore should be avoided.

If you want to allow scalar iteration, use `allowscalar` or `@allowscalar`
to enable scalar iteration globally or for the operations in question.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] errorscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
  [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
  [4] assertscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
  [5] getindex(A::CuArray{Int32, 1, CUDA.DeviceMemory}, I::Int64)
    @ GPUArrays ~/.julia/packages/GPUArrays/uiVyU/src/host/indexing.jl:50
  [6] getindex(A::CUDA.CUSPARSE.CuSparseMatrixCSR{ComplexF64, Int32}, i0::Int64, i1::Int64)
    @ CUDA.CUSPARSE ~/.julia/packages/CUDA/1kIOw/lib/cusparse/array.jl:398
  [7] getindex
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:280 [inlined]
  [8] getindex
    @ ./subarray.jl:320 [inlined]
  [9] iterate
    @ ./abstractarray.jl:1209 [inlined]
 [10] iterate
    @ ./abstractarray.jl:1207 [inlined]
 [11] _all
    @ ./reduce.jl:1309 [inlined]
 [12] all
    @ ./reducedim.jl:995 [inlined]
 [13] _istriu
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1317 [inlined]
 [14] istriu(A::UnitLowerTriangular{ComplexF64, CUDA.CUSPARSE.CuSparseMatrixCSR{ComplexF64, Int32}}, k::Int64)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1310
 [15] istriu
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1309 [inlined]
 [16] ldiv!
    @ ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:972 [inlined]
 [17] ldiv!
    @ ~/.julia/packages/KrylovPreconditioners/dc3YZ/ext/CUDA/ilu0.jl:78 [inlined]
 [18] macro expansion
    @ ~/.julia/packages/KrylovPreconditioners/dc3YZ/ext/CUDA/ilu0.jl:86 [inlined]
 [19] macro expansion
    @ ./timing.jl:421 [inlined]
 [20] ldiv!(y::CuArray{…}, ilu::KrylovPreconditionersCUDAExt.NVIDIA_ILU0{…}, x::CuArray{…})
    @ 

> KrylovPreconditionersCUDAExt ~/.julia/packages/KrylovPreconditioners/dc3YZ/ext/CUDA/ilu0.jl:85

 [21] mulorldiv!
    @ ~/.julia/packages/Krylov/wgV9p/src/krylov_utils.jl:307 [inlined]
 [22] gmres!(solver::Krylov.GmresSolver{…}, A::CUDA.CUSPARSE.CuSparseMatrixCSR{…}, b::CuArray{…}; M::UniformScaling{…}, N::KrylovPreconditionersCUDAExt.NVIDIA_ILU0{…}, ldiv::Bool, restart::Bool, reorthogonalization::Bool, atol::Float64, rtol::Float64, itmax::Int64, timemax::Float64, verbose::Int64, history::Bool, callback::Krylov.var"#987#995", iostream::Core.CoreSTDOUT)
    @ Krylov ~/.julia/packages/Krylov/wgV9p/src/gmres.jl:241
 [23] gmres!
    @ ~/.julia/packages/Krylov/wgV9p/src/gmres.jl:108 [inlined]
 [24] solve!
    @ ~/.julia/packages/Krylov/wgV9p/src/krylov_solve.jl:159 [inlined]
 [25] solve!(cache::LinearSolve.LinearCache{…}, alg::LinearSolve.KrylovJL{…}; kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/packages/LinearSolve/QuCP8/src/iterative_wrappers.jl:296
 [26] solve!
    @ ~/.julia/packages/LinearSolve/QuCP8/src/iterative_wrappers.jl:247 [inlined]
 [27] #solve!#8
    @ ~/.julia/packages/LinearSolve/QuCP8/src/common.jl:307 [inlined]
 [28] solve!
    @ ~/.julia/packages/LinearSolve/QuCP8/src/common.jl:306 [inlined]
 [29] #solve#7
    @ ~/.julia/packages/LinearSolve/QuCP8/src/common.jl:303 [inlined]
 [30] _steadystate(L::QuantumObject{…}, solver::SteadyStateLinearSolver{…}; kwargs::@Kwargs{})
    @ QuantumToolbox ~/.julia/packages/QuantumToolbox/YXgET/src/steadystate.jl:158
 [31] _steadystate
    @ ~/.julia/packages/QuantumToolbox/YXgET/src/steadystate.jl:130 [inlined]
 [32] #steadystate#366
    @ ~/.julia/packages/QuantumToolbox/YXgET/src/steadystate.jl:127 [inlined]
 [33] steadystate
    @ ~/.julia/packages/QuantumToolbox/YXgET/src/steadystate.jl:113 [inlined]
 [34] top-level scope
    @ Untitled-1:27
Some type information was truncated. Use `show(err)` to see complete types.

Line 27 is the line of pss_gpu with the sparse gpu steadystate() function.

[ytdHuang]

Could you check whether all the related packages in your environment are up-to-date ?

[albertomercurio]

It seems to be related to this issue. I fixed it installing the dev version of CUDA.jl with

(test) pkg> add CUDA#master

Until CUDA.jl releases a new version.

[ytdHuang]

@MatthewBelzer @albertomercurio

CUDA.jl has released new a version. Can the above problem be solved ?

[MatthewBelzer]

Hi, the GPU code above works now. Thanks for the help!