Merge branch 'main' into liouvillian2

Author: ytdHuang

CHANGELOG.md

@@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+- 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]) 
 - Generalize the definition of `liouvillian`. It no longer expects the Hamiltonian to be Hermitian. ([#581])
 
 ## [v0.38.1]
@@ -359,4 +360,5 @@ Release date: 2024-11-13
 [#575]: https://github.com/qutip/QuantumToolbox.jl/issues/575
 [#576]: https://github.com/qutip/QuantumToolbox.jl/issues/576
 [#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

src/QuantumToolbox.jl

@@ -15,6 +15,7 @@ import SciMLBase:
     remake,
     u_modified!,
     NullParameters,
+    LinearProblem,
     ODEFunction,
     SDEFunction,
     ODEProblem,
@@ -46,7 +47,8 @@ import SciMLOperators:
     IdentityOperator,
     update_coefficients!,
     concretize
-import LinearSolve: LinearProblem, SciMLLinearSolveAlgorithm, KrylovJL_MINRES, KrylovJL_GMRES
+import LinearSolve:
+    SciMLLinearSolveAlgorithm, KrylovJL_MINRES, KrylovJL_GMRES, UMFPACKFactorization, OperatorAssumptions
 import DiffEqBase: get_tstops
 import DiffEqCallbacks: PeriodicCallback, FunctionCallingCallback, FunctionCallingAffect, TerminateSteadyState
 import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm

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,17 +370,24 @@ function eigsolve(
         vals = res.values
     else
         Aₛ = A - sigma * I
-        solver === nothing && (solver = isH ? KrylovJL_MINRES() : 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)
-        linsolve = init(prob, solver; kwargs2...)
+        prob = LinearProblem{true}(Aₛ, v0)
+
+        linsolve = if isnothing(solver)
+            if typeof(A) <: SparseMatrixCSC && isprimitivetype(T)
+                init(prob, UMFPACKFactorization(); kwargs2...)
+            else
+                init(prob; kwargs2...)
+            end
+        else
+            init(prob, solver; kwargs2...)
+        end
 
         Amap = EigsolveInverseMap(T, size(A), linsolve)
 

test/ext-test/cpu/arbitrary_precision/Project.toml

@@ -1,3 +1,4 @@
 [deps]
 GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
 QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
+Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"

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()
@@ -88,7 +89,9 @@ if (GROUP == "Arbitrary-Precision")
     Pkg.update()
 
     using QuantumToolbox
+    using LinearAlgebra
     using SparseArrays
+    using Sparspak
     using GenericSchur
 
     QuantumToolbox.about()