Changed TimeEvolution to use KrylovKit instead of ExponentialUtilities

lib/YaoBlocks/Project.toml

@@ -1,14 +1,14 @@
 name = "YaoBlocks"
 uuid = "418bc28f-b43b-5e0b-a6e7-61bbc1a2c1df"
-version = "0.13.9"
+version = "0.13.10"
 
 [deps]
 BitBasis = "50ba71b6-fa0f-514d-ae9a-0916efc90dcf"
 CacheServers = "a921213e-d44a-5460-ac04-5d720a99ba71"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LegibleLambdas = "f1f30506-32fe-5131-bd72-7c197988f9e5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LuxurySparse = "d05aeea4-b7d4-55ac-b691-9e7fabb07ba2"
@@ -26,7 +26,7 @@ BitBasis = "0.8"
 CacheServers = "0.2"
 ChainRulesCore = "1.11"
 DocStringExtensions = "0.8, 0.9"
-ExponentialUtilities = "1.5"
+KrylovKit = "0.5, 0.6"
 LegibleLambdas = "0.2, 0.3"
 LuxurySparse = "0.7"
 MLStyle = "0.3, 0.4"

lib/YaoBlocks/src/YaoBlocks.jl

@@ -14,7 +14,8 @@ using StatsBase, TupleTools, InteractiveUtils
 using MLStyle: @match
 using LinearAlgebra: eigen!
 using SparseArrays, LuxurySparse
-using ExponentialUtilities, Random, CacheServers
+using Random, CacheServers
+import KrylovKit: exponentiate
 using DocStringExtensions
 import StaticArrays: SMatrix
 

lib/YaoBlocks/src/primitive/time_evolution.jl

@@ -87,21 +87,27 @@ function LinearAlgebra.mul!(y::AbstractVector, A::BlockMap{T,GT}, x::AbstractVec
     return y
 end
 
+compatible_multiplicative_operand(::AbstractArray, source::AbstractArray) = source
+
 function YaoAPI.unsafe_apply!(reg::AbstractArrayReg{D,T}, te::TimeEvolution) where {D,T}
     if isdiagonal(te.H)
         # `compatible_multiplicative_operand` is used to ensure GPU compatibility.
-        reg.state .*= exp.((-im * te.dt) .* ExponentialUtilities.compatible_multiplicative_operand(reg.state, diag(mat(T, te.H))))
+        reg.state .*= exp.((-im * te.dt) .* compatible_multiplicative_operand(reg.state, diag(mat(T, te.H))))
         return reg
     end
     st = state(reg)
     dt = real(te.dt) == 0 ? imag(te.dt) : -im * te.dt
     A = BlockMap(T, te.H)
     @inbounds for j = 1:size(st, 2)
         v = view(st, :, j)
-        # TODO: opnorm could be estimated by inspecting the block operation.
-        # We should fix it in the future in improve the accuracy.
-        Ks = arnoldi(A, v; tol = te.tol, ishermitian = true, opnorm = 1.0)
-        expv!(v, dt, Ks)
+        v_out, info = exponentiate(A, dt, v,
+                                   tol=te.tol,
+                                   krylovdim=min(1000, size(A,1)),
+                                   ishermitian=true,
+                                   eager=true,)
+        # info.converged is 1 if it converged and 0 otherwise
+        Bool(info.converged) || @warn "Application of $(typeof(te)) did not converge. error = $(info.normres) but tol=$(te.tol)"
+        v .= v_out
     end
     return reg
 end
@@ -124,11 +130,8 @@ function Base.adjoint(te::TimeEvolution)
 end
 Base.copy(te::TimeEvolution) = TimeEvolution(te.H, te.dt, tol = te.tol)
 
-YaoAPI.isdiagonal(r::TimeEvolution) = isdiagonal(r.H)
-function YaoAPI.isunitary(te::TimeEvolution)
-    iszero(imag(te.dt)) || return false
-    return true
-end
+YaoAPI.isdiagonal(te::TimeEvolution) = isdiagonal(te.H)
+YaoAPI.isunitary(te::TimeEvolution) = iszero(imag(te.dt))
 
 iparams_range(::TimeEvolution{D,T}) where {D,T} = ((typemin(T), typemax(T)),)