Change default solver for eigsolve and add GPU tests

Author: albertomercurio

Project.toml

@@ -11,7 +11,6 @@ DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
-IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
@@ -52,7 +51,6 @@ Distributed = "1"
 FFTW = "1.5"
 GPUArrays = "10, 11"
 Graphs = "1.7"
-IncompleteLU = "0.2"
 KernelAbstractions = "0.9"
 LaTeXStrings = "1.2"
 LinearAlgebra = "1"

src/QuantumToolbox.jl

@@ -46,7 +46,8 @@ import SciMLOperators:
     IdentityOperator,
     update_coefficients!,
     concretize
-import LinearSolve: LinearProblem, SciMLLinearSolveAlgorithm, KrylovJL_MINRES, KrylovJL_GMRES
+import LinearSolve:
+    LinearProblem, SciMLLinearSolveAlgorithm, KrylovJL_MINRES, KrylovJL_GMRES, UMFPACKFactorization, OperatorAssumptions
 import DiffEqBase: get_tstops
 import DiffEqCallbacks: PeriodicCallback, FunctionCallingCallback, FunctionCallingAffect, TerminateSteadyState
 import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm
@@ -58,7 +59,6 @@ import ArrayInterface: allowed_getindex, allowed_setindex!
 import Distributed: RemoteChannel
 import FFTW: fft, ifft, fftfreq, fftshift
 import Graphs: connected_components, DiGraph
-import IncompleteLU: ilu
 import LaTeXStrings: @L_str
 import Pkg
 import ProgressMeter: Progress, next!

src/qobj/eigsolve.jl

@@ -352,10 +352,9 @@ function eigsolve(
     kwargs...,
 )
     T = eltype(A)
-    isH = ishermitian(A)
-    v0 === nothing && (v0 = normalize!(rand(T, size(A, 1))))
+    isnothing(v0) && (v0 = normalize!(rand(T, size(A, 1))))
 
-    if sigma === nothing
+    if isnothing(sigma)
         res = _eigsolve(
             A,
             v0,
@@ -371,16 +370,15 @@ function eigsolve(
         vals = res.values
     else
         Aₛ = A - sigma * I
-        solver === nothing && (solver = isH ? KrylovJL_MINRES() : KrylovJL_GMRES())
+        isnothing(solver) && (solver = typeof(A) <: SparseMatrixCSC ? UMFPACKFactorization() : KrylovJL_GMRES())
 
         kwargs2 = (; kwargs...)
-        condition = !haskey(kwargs2, :Pl) && typeof(A) <: SparseMatrixCSC
-        condition && (kwargs2 = merge(kwargs2, (Pl = ilu(Aₛ, τ = 0.01),)))
 
         !haskey(kwargs2, :abstol) && (kwargs2 = merge(kwargs2, (abstol = tol * 1e-6,)))
         !haskey(kwargs2, :reltol) && (kwargs2 = merge(kwargs2, (reltol = tol * 1e-6,)))
+        !haskey(kwargs2, :assumptions) && (kwargs2 = merge(kwargs2, (assumptions = OperatorAssumptions(true),)))
 
-        prob = LinearProblem(Aₛ, v0)
+        prob = LinearProblem{true}(Aₛ, v0)
         linsolve = init(prob, solver; kwargs2...)
 
         Amap = EigsolveInverseMap(T, size(A), linsolve)

test/ext-test/gpu/Project.toml

@@ -1,6 +1,8 @@
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
 
 [compat]
-CUDA = "5"
+CUDA = "5"

test/ext-test/gpu/cuda_ext.jl

@@ -264,3 +264,35 @@ end
         @inferred ptrace(ψ, 2)
     end
 end
+
+@testset "CUDA eigsolve" begin
+    N = 30
+    Δ = 0.5
+    U = 0.1
+    κ = 0.1
+    F = 0.5
+
+    a = destroy(N)
+    H = Δ * a' * a + U / 2 * a' * a' * a * a + F * (a + a')
+
+    c_ops = [sqrt(κ) * a]
+
+    L = liouvillian(H, c_ops)
+    L_gpu = CuSparseMatrixCSR(L)
+
+    vals_cpu, vecs_cpu = eigenstates(L; sparse = true, sigma = 0.01, eigvals = 4, krylovdim = 30)
+    vals_gpu, vecs_gpu = eigenstates(
+        L_gpu;
+        sparse = true,
+        sigma = 0.01,
+        eigvals = 4,
+        krylovdim = 30,
+        solver = LUFactorization(),
+        v0 = CUDA.rand(ComplexF64, size(L_gpu, 1)),
+    )
+
+    @test vals_cpu ≈ vals_gpu atol = 1e-8
+    @test all(zip(vecs_cpu, vecs_gpu)) do (v_cpu, v_gpu)
+        return isapprox(abs(dot(v_cpu.data, Array(v_gpu.data))), 1; atol = 1e-8)
+    end
+end

test/runtests.jl

@@ -74,7 +74,8 @@ if (GROUP == "CUDA_Ext")
     import StaticArraysCore: SVector
     using CUDA
     using CUDA.CUSPARSE
-    # CUDA.allowscalar(false) # This is already set in the extension script
+    using CUDSS
+    using LinearSolve
 
     QuantumToolbox.about()
     CUDA.versioninfo()