FIx eigenvalues

Author: albertomercurio

Project.toml

@@ -11,7 +11,6 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
@@ -39,7 +38,6 @@ FFTW = "1.5"
 Graphs = "1.7"
 IncompleteLU = "0.2"
 LinearAlgebra = "<0.0.1, 1"
-LinearMaps = "3"
 LinearSolve = "2"
 OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqTsit5 = "1"

docs/src/api.md

@@ -219,6 +219,12 @@ wigner
 negativity
 ```
 
+## [Linear Maps](@id doc-API:Linear-Maps)
+
+```@docs
+AbstractLinearMap
+```
+
 ## [Utility functions](@id doc-API:Utility-functions)
 
 ```@docs

src/QuantumToolbox.jl

@@ -40,7 +40,6 @@ import ArrayInterface: allowed_getindex, allowed_setindex!
 import FFTW: fft, fftshift
 import Graphs: connected_components, DiGraph
 import IncompleteLU: ilu
-import LinearMaps: LinearMap
 import Pkg
 import Random
 import SpecialFunctions: loggamma
@@ -54,6 +53,7 @@ BLAS.set_num_threads(1)
 include("utilities.jl")
 include("versioninfo.jl")
 include("progress_bar.jl")
+include("linear_maps.jl")
 
 # Quantum Object
 include("qobj/quantum_object.jl")

src/linear_maps.jl

@@ -0,0 +1,54 @@
+export AbstractLinearMap
+
+@doc raw"""
+    AbstractLinearMap{T, TS}
+
+Represents a general linear map with element type `T` and size `TS`.
+
+## Overview
+
+A **linear map** is a transformation `L` that satisfies:
+
+- **Additivity**: 
+    ```math
+    L(u + v) = L(u) + L(v)
+    ```
+- **Homogeneity**: 
+    ```math
+    L(cu) = cL(u)
+    ```
+
+It is typically represented as a matrix with dimensions given by `size`, and this abtract type helps to define this map when the matrix is not explicitly available.
+
+## Methods
+
+- `Base.eltype(A)`: Returns the element type `T`.
+- `Base.size(A)`: Returns the size `A.size`.
+- `Base.size(A, i)`: Returns the `i`-th dimension.
+
+## Example
+
+As an example, we now define the linear map used in the [`eigsolve_al`](@ref) function for Arnoldi-Lindblad eigenvalue solver:
+
+```julia-repl
+struct ArnoldiLindbladIntegratorMap{T,TS,TI} <: AbstractLinearMap{T,TS}
+    elty::Type{T}
+    size::TS
+    integrator::TI
+end
+
+function LinearAlgebra.mul!(y::AbstractVector, A::ArnoldiLindbladIntegratorMap, x::AbstractVector)
+    reinit!(A.integrator, x)
+    solve!(A.integrator)
+    return copyto!(y, A.integrator.u)
+end
+```
+
+where `integrator` is the ODE integrator for the time-evolution. In this way, we can diagonalize this linear map using the [`eigsolve`](@ref) function.
+"""
+abstract type AbstractLinearMap{T,TS} end
+
+Base.eltype(A::AbstractLinearMap{T}) where T = T
+
+Base.size(A::AbstractLinearMap) = A.size
+Base.size(A::AbstractLinearMap, i::Int) = A.size[i]

src/qobj/eigsolve.jl

@@ -84,9 +84,27 @@ function Base.show(io::IO, res::EigsolveResult)
     return show(io, MIME("text/plain"), res.vectors)
 end
 
-function _map_ldiv(linsolve, y, x)
-    linsolve.b .= x
-    return y .= solve!(linsolve).u
+struct ArnoldiLindbladIntegratorMap{T,TS,TI} <: AbstractLinearMap{T,TS}
+    elty::Type{T}
+    size::TS
+    integrator::TI
+end
+
+function LinearAlgebra.mul!(y::AbstractVector, A::ArnoldiLindbladIntegratorMap, x::AbstractVector)
+    reinit!(A.integrator, x)
+    solve!(A.integrator)
+    return copyto!(y, A.integrator.u)
+end
+
+struct EigsolveInverseMap{T,TS,TI} <: AbstractLinearMap{T,TS}
+    elty::Type{T}
+    size::TS
+    linsolve::TI
+end
+
+function LinearAlgebra.mul!(y::AbstractVector, A::EigsolveInverseMap, x::AbstractVector)
+    A.linsolve.b .= x
+    return copyto!(y, solve!(A.linsolve).u)
 end
 
 function _permuteschur!(
@@ -293,7 +311,8 @@ function eigsolve(
 
         prob = LinearProblem(Aₛ, v0)
         linsolve = init(prob, solver; kwargs2...)
-        Amap = LinearMap{T}((y, x) -> _map_ldiv(linsolve, y, x), length(v0))
+
+        Amap = EigsolveInverseMap(T, size(A), linsolve)
 
         res = _eigsolve(Amap, v0, type, dims, k, krylovdim, tol = tol, maxiter = maxiter)
         vals = @. (1 + sigma * res.values) / res.values
@@ -370,13 +389,7 @@ function eigsolve_al(
 
     # prog = ProgressUnknown(desc="Applications:", showspeed = true, enabled=progress)
 
-    function arnoldi_lindblad_solve(ρ)
-        reinit!(integrator, ρ)
-        solve!(integrator)
-        return integrator.u
-    end
-
-    Lmap = LinearMap{eltype(MT1)}(arnoldi_lindblad_solve, size(L, 1), ismutating = false)
+    Lmap = ArnoldiLindbladIntegratorMap(eltype(MT1),size(L),integrator)
 
     res = _eigsolve(Lmap, mat2vec(ρ0), L.type, L.dims, k, krylovdim, maxiter = maxiter, tol = eigstol)
     # finish!(prog)

test/eigenvalues_and_operators.jl

@@ -59,7 +59,7 @@
     # eigen solve for general matrices
     vals, _, vecs = eigsolve(L.data, sigma = 0.01, k = 10, krylovdim = 50)
     vals2, vecs2 = eigen(sparse_to_dense(L.data))
-    vals3, state3, vecs3 = eigsolve_al(liouvillian(H, c_ops), 1 \ (40 * κ), k = 10, krylovdim = 50)
+    vals3, state3, vecs3 = eigsolve_al(L, 1 \ (40 * κ), k = 10, krylovdim = 50)
     idxs = sortperm(vals2, by = abs)
     vals2 = vals2[idxs][1:10]
     vecs2 = vecs2[:, idxs][:, 1:10]
@@ -91,4 +91,27 @@
     @test isapprox(abs2(vals2[1]), abs2(vals3[1]), atol = 1e-7)
     @test isapprox(vec2mat(vecs[1]).data * exp(-1im * angle(vecs[1][1])), vec2mat(vecs2[1]).data, atol = 1e-7)
     @test isapprox(vec2mat(vecs[1]).data * exp(-1im * angle(vecs[1][1])), vec2mat(state3[1]).data, atol = 1e-5)
+
+    @testset "Type Inference (eigen)" begin
+      N = 5
+      a = kron(destroy(N), qeye(N))
+      a_d = a'
+      b = kron(qeye(N), destroy(N))
+      b_d = b'
+  
+      ωc = 1
+      ωb = 1
+      g = 0.01
+      κ = 0.1
+      n_thermal = 0.01
+  
+      H = ωc * a_d * a + ωb * b_d * b + g * (a + a_d) * (b + b_d)
+      c_ops = [√((1 + n_thermal) * κ) * a, √κ * b, √(n_thermal * κ) * a_d]
+      L = liouvillian(H, c_ops)
+
+      @inferred eigenstates(H, sparse = false)
+      @inferred eigenstates(H, sparse = true)
+      @inferred eigenstates(L, sparse = true)
+      @inferred eigsolve_al(L, 1 \ (40 * κ), k = 10)
+    end
 end

test/generalized_master_equation.jl

@@ -26,7 +26,7 @@
     c_ops = [sqrt(0.01) * a2, sqrt(0.01) * sm2]
     L2 = liouvillian(H_d, c_ops)
 
-    @test (expect(Xp' * Xp, steadystate(L1)) < 1e-10 && expect(Xp' * Xp, steadystate(L2)) > 1e-3)
+    @test (expect(Hermitian(Xp' * Xp), steadystate(L1)) < 1e-10 && expect(Hermitian(Xp' * Xp), steadystate(L2)) > 1e-3)
 
     H = 1 * a' * a + 1 * sz / 2 + 1e-5 * (a * sp + a' * sm)