Merge branch 'cite' of https://github.com/qutip/QuantumToolbox.jl into cite

Author: ytdHuang

.github/workflows/Code-Quality.yml

@@ -40,6 +40,7 @@ jobs:
       matrix:
         version:
           - 'lts'
+          - '1'
         os:
           - 'ubuntu-latest'
         arch:

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v5
       - name: Check spelling
-        uses: crate-ci/[email protected].2
+        uses: crate-ci/[email protected].3

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
 - Add `QuantumToolbox.cite()` for bibtex generator of `QuantumToolbox.jl`. ([#544])
+- Add `sortby` and `rev` keyword arguments to eigensolvers. ([#546])
 
 ## [v0.35.0]
 Release date: 2025-09-03
@@ -317,3 +318,4 @@ Release date: 2024-11-13
 [#537]: https://github.com/qutip/QuantumToolbox.jl/issues/537
 [#539]: https://github.com/qutip/QuantumToolbox.jl/issues/539
 [#544]: https://github.com/qutip/QuantumToolbox.jl/issues/544
+[#546]: https://github.com/qutip/QuantumToolbox.jl/issues/546

Makefile

@@ -12,10 +12,10 @@ changelog:
 	${JULIA} -e 'using Changelog; Changelog.generate(Changelog.CommonMark(), "CHANGELOG.md"; repo = "qutip/QuantumToolbox.jl")'
 
 test:
-	${JULIA} --project -e 'using Pkg; Pkg.resolve(); Pkg.test()'
+	${JULIA} --project -e 'using Pkg; Pkg.update(); Pkg.test()'
 
 docs:
-	${JULIA} --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+	${JULIA} --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.update()'
 	${JULIA} --project=docs docs/make.jl
 
 vitepress:

ext/QuantumToolboxMakieExt.jl

@@ -24,10 +24,8 @@ import Makie:
     scatter!,
     arrows3d!,
     text!,
-    Point3f,
     mesh!,
     RGBf,
-    Vec3f,
     Point3f,
     NoShading,
     cameracontrols,
@@ -660,7 +658,7 @@ function _plot_arcs!(b::Bloch, lscene)
             dot(cross(v1, vm), n) < 0 && (θ -= 2π)
         end
         t_range = range(0, θ, length = 100)
-        arc_points = [Point3f((v1 * cos(t) + cross(n, v1) * sin(t))) for t in t_range]
+        arc_points = [Point3f(v1 * cos(t) + cross(n, v1) * sin(t)) for t in t_range]
         lines!(lscene, arc_points; color = "blue", linestyle = :solid)
     end
     return nothing
@@ -685,13 +683,13 @@ function _plot_vectors!(b::Bloch, lscene)
 
         arrows3d!(
             lscene,
-            [Point3f(0, 0, 0)],
-            [v],
+            Point3f(0),
+            Point3f(v),
             color = color,
             shaftradius = b.vector_width,
             tiplength = b.vector_tiplength,
             tipradius = b.vector_tipradius,
-            rasterize = 3,
+            # rasterize = 3, #TODO: maybe uncomment this after https://github.com/MakieOrg/Makie.jl/issues/5259 is fixed
         )
     end
     return nothing

src/qobj/boolean_functions.jl

@@ -10,53 +10,47 @@ export isunitary
 
 Checks if the [`QuantumObject`](@ref) `A` is a [`Bra`](@ref). Default case returns `false` for any other inputs.
 """
-isbra(A::QuantumObject) = isbra(typeof(A))
-isbra(::Type{<:QuantumObject{Bra,N}}) where {N} = true
+isbra(A::QuantumObject{Bra}) = true
 isbra(A) = false # default case
 
 @doc raw"""
     isket(A)
 
 Checks if the [`QuantumObject`](@ref) `A` is a [`Ket`](@ref). Default case returns `false` for any other inputs.
 """
-isket(A::QuantumObject) = isket(typeof(A))
-isket(::Type{<:QuantumObject{Ket,N}}) where {N} = true
+isket(A::QuantumObject{Ket}) = true
 isket(A) = false # default case
 
 @doc raw"""
     isoper(A)
 
 Checks if the [`AbstractQuantumObject`](@ref) `A` is a [`Operator`](@ref). Default case returns `false` for any other inputs.
 """
-isoper(A::AbstractQuantumObject) = isoper(typeof(A))
-isoper(::Type{<:AbstractQuantumObject{Operator,N}}) where {N} = true
+isoper(A::AbstractQuantumObject{Operator}) = true
 isoper(A) = false # default case
 
 @doc raw"""
     isoperbra(A)
 
 Checks if the [`QuantumObject`](@ref) `A` is a [`OperatorBra`](@ref). Default case returns `false` for any other inputs.
 """
-isoperbra(A::QuantumObject) = isoperbra(typeof(A))
-isoperbra(::Type{<:QuantumObject{OperatorBra,N}}) where {N} = true
+isoperbra(A::QuantumObject{OperatorBra}) = true
 isoperbra(A) = false # default case
 
 @doc raw"""
     isoperket(A)
 
 Checks if the [`QuantumObject`](@ref) `A` is a [`OperatorKet`](@ref). Default case returns `false` for any other inputs.
 """
-isoperket(A::QuantumObject) = isoperket(typeof(A))
-isoperket(::Type{<:QuantumObject{OperatorKet,N}}) where {N} = true
+isoperket(A::QuantumObject{OperatorKet}) = true
 isoperket(A) = false # default case
 
 @doc raw"""
     issuper(A)
 
 Checks if the [`AbstractQuantumObject`](@ref) `A` is a [`SuperOperator`](@ref). Default case returns `false` for any other inputs.
 """
-issuper(A::AbstractQuantumObject) = issuper(typeof(A))
-issuper(::Type{<:AbstractQuantumObject{<:SuperOperatorType,N}}) where {N} = true
+issuper(A::AbstractQuantumObject{<:SuperOperatorType}) = true
 issuper(A) = false # default case
 
 @doc raw"""

src/qobj/eigsolve.jl

@@ -149,12 +149,12 @@ function _permuteschur!(
     return T, Q
 end
 
-function _update_schur_eigs!(Hₘ, Uₘ, Uₘᵥ, f, m, β, sorted_vals)
+function _update_schur_eigs!(Hₘ, Uₘ, Uₘᵥ, f, m, β, sorted_vals, sortby, rev)
     F = hessenberg!(Hₘ)
     copyto!(Uₘ, Hₘ)
     LAPACK.orghr!(1, m, Uₘ, F.τ)
     Tₘ, Uₘ, values = hseqr!(Hₘ, Uₘ)
-    sortperm!(sorted_vals, values, by = abs, rev = true)
+    sortperm!(sorted_vals, values, by = sortby, rev = rev)
     _permuteschur!(Tₘ, Uₘ, sorted_vals)
     mul!(f, Uₘᵥ, β)
 
@@ -170,6 +170,8 @@ function _eigsolve(
     m::Int = max(20, 2 * k + 1);
     tol::Real = 1e-8,
     maxiter::Int = 200,
+    sortby::Function = abs2,
+    rev = true,
 ) where {T<:BlasFloat,ObjType<:Union{Nothing,Operator,SuperOperator}}
     n = size(A, 2)
     V = similar(b, n, m + 1)
@@ -209,7 +211,7 @@ function _eigsolve(
 
     M = typeof(cache0)
 
-    Tₘ, Uₘ = _update_schur_eigs!(Hₘ, Uₘ, Uₘᵥ, f, m, β, sorted_vals)
+    Tₘ, Uₘ = _update_schur_eigs!(Hₘ, Uₘ, Uₘᵥ, f, m, β, sorted_vals, sortby, rev)
 
     numops = m
     iter = 0
@@ -235,7 +237,7 @@ function _eigsolve(
 
         # println( A * Vₘ ≈ Vₘ * M(Hₘ) + qₘ * M(transpose(βeₘ)) )     # SHOULD BE TRUE
 
-        Tₘ, Uₘ = _update_schur_eigs!(Hₘ, Uₘ, Uₘᵥ, f, m, β, sorted_vals)
+        Tₘ, Uₘ = _update_schur_eigs!(Hₘ, Uₘ, Uₘᵥ, f, m, β, sorted_vals, sortby, rev)
 
         numops += m - k - 1
         iter += 1
@@ -263,6 +265,8 @@ end
         tol::Real = 1e-8,
         maxiter::Int = 200,
         solver::Union{Nothing, SciMLLinearSolveAlgorithm} = nothing,
+        sortby::Function = abs2,
+        rev::Bool = true,
         kwargs...)
 
 Solve for the eigenvalues and eigenvectors of a matrix `A` using the Arnoldi method.
@@ -276,6 +280,8 @@ Solve for the eigenvalues and eigenvectors of a matrix `A` using the Arnoldi met
 - `tol::Real`: the tolerance for the Arnoldi method. Default is `1e-8`.
 - `maxiter::Int`: the maximum number of iterations for the Arnoldi method. Default is `200`.
 - `solver::Union{Nothing, SciMLLinearSolveAlgorithm}`: the linear solver algorithm. Default is `nothing`.
+- `sortby::Function`: the function to sort eigenvalues. Default is `abs2`.
+- `rev::Bool`: whether to sort in descending order. Default is `true`.
 - `kwargs`: Additional keyword arguments passed to the solver.
 
 # Notes
@@ -293,6 +299,8 @@ function eigsolve(
     tol::Real = 1e-8,
     maxiter::Int = 200,
     solver::Union{Nothing,SciMLLinearSolveAlgorithm} = nothing,
+    sortby::Function = abs2,
+    rev::Bool = true,
     kwargs...,
 )
     return eigsolve(
@@ -306,6 +314,8 @@ function eigsolve(
         tol = tol,
         maxiter = maxiter,
         solver = solver,
+        sortby = sortby,
+        rev = rev,
         kwargs...,
     )
 end
@@ -321,14 +331,27 @@ function eigsolve(
     tol::Real = 1e-8,
     maxiter::Int = 200,
     solver::Union{Nothing,SciMLLinearSolveAlgorithm} = nothing,
+    sortby::Function = abs2,
+    rev::Bool = true,
     kwargs...,
 )
     T = eltype(A)
     isH = ishermitian(A)
     v0 === nothing && (v0 = normalize!(rand(T, size(A, 1))))
 
     if sigma === nothing
-        res = _eigsolve(A, v0, type, dimensions, eigvals, krylovdim, tol = tol, maxiter = maxiter)
+        res = _eigsolve(
+            A,
+            v0,
+            type,
+            dimensions,
+            eigvals,
+            krylovdim,
+            tol = tol,
+            maxiter = maxiter,
+            sortby = sortby,
+            rev = rev,
+        )
         vals = res.values
     else
         Aₛ = A - sigma * I
@@ -346,7 +369,18 @@ function eigsolve(
 
         Amap = EigsolveInverseMap(T, size(A), linsolve)
 
-        res = _eigsolve(Amap, v0, type, dimensions, eigvals, krylovdim, tol = tol, maxiter = maxiter)
+        res = _eigsolve(
+            Amap,
+            v0,
+            type,
+            dimensions,
+            eigvals,
+            krylovdim,
+            tol = tol,
+            maxiter = maxiter,
+            sortby = sortby,
+            rev = rev,
+        )
         vals = @. (1 + sigma * res.values) / res.values
     end
 
@@ -368,6 +402,8 @@ end
         krylovdim::Int = min(10, size(H, 1)),
         maxiter::Int = 200,
         eigstol::Real = 1e-6,
+        sortby::Function = abs2,
+        rev::Bool = true,
         kwargs...,
     )
 
@@ -384,6 +420,8 @@ Solve the eigenvalue problem for a Liouvillian superoperator `L` using the Arnol
 - `krylovdim`: The dimension of the Krylov subspace.
 - `maxiter`: The maximum number of iterations for the eigsolver.
 - `eigstol`: The tolerance for the eigsolver.
+- `sortby::Function`: the function to sort eigenvalues. Default is `abs2`.
+- `rev::Bool`: whether to sort in descending order. Default is `true`.
 - `kwargs`: Additional keyword arguments passed to the differential equation solver.
 
 # Notes
@@ -407,6 +445,8 @@ function eigsolve_al(
     krylovdim::Int = min(10, size(H, 1)),
     maxiter::Int = 200,
     eigstol::Real = 1e-6,
+    sortby::Function = abs2,
+    rev::Bool = true,
     kwargs...,
 ) where {HOpType<:Union{Operator,SuperOperator}}
     L_evo = _mesolve_make_L_QobjEvo(H, c_ops)
@@ -422,8 +462,18 @@ function eigsolve_al(
 
     Lmap = ArnoldiLindbladIntegratorMap(eltype(H), size(L_evo), integrator)
 
-    res =
-        _eigsolve(Lmap, mat2vec(ρ0), L_evo.type, L_evo.dimensions, eigvals, krylovdim, maxiter = maxiter, tol = eigstol)
+    res = _eigsolve(
+        Lmap,
+        mat2vec(ρ0),
+        L_evo.type,
+        L_evo.dimensions,
+        eigvals,
+        krylovdim,
+        maxiter = maxiter,
+        tol = eigstol,
+        sortby = sortby,
+        rev = rev,
+    )
 
     vals = similar(res.values)
     vecs = similar(res.vectors)