From f80036acc0e08f110a4105105266e1c365541060 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 20 Jul 2026 15:05:11 -0400 Subject: [PATCH 1/3] use backends and allocators in DMRG --- .../derivatives/hamiltonian_derivatives.jl | 10 +++-- src/algorithms/derivatives/mpo_derivatives.jl | 21 +++++++--- .../derivatives/projection_derivatives.jl | 13 +++++-- src/algorithms/groundstate/dmrg.jl | 38 +++++++++++++++---- 4 files changed, 60 insertions(+), 22 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index 438fe87ec..e9be54228 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -37,14 +37,15 @@ end function AC_hamiltonian( site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian, above::_HAM_MPS_TYPES, envs; - prepare::Bool = true + prepare::Bool = true, + backend::AbstractBackend = DefaultBackend(), allocator = BufferAllocator() ) @assert below === above "JordanMPO assumptions break" GL = leftenv(envs, site, below) GR = rightenv(envs, site, below) W = operator[site] H_AC = JordanMPO_AC_Hamiltonian(GL, W, GR) - return prepare ? prepare_operator!!(H_AC) : H_AC + return prepare ? prepare_operator!!(H_AC, backend, allocator) : H_AC end function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTensor) @@ -189,14 +190,15 @@ end function AC2_hamiltonian( site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian, above::_HAM_MPS_TYPES, envs; - prepare::Bool = true + prepare::Bool = true, + backend::AbstractBackend = DefaultBackend(), allocator = BufferAllocator() ) @assert below === above "JordanMPO assumptions break" GL = leftenv(envs, site, below) GR = rightenv(envs, site + 1, below) W1, W2 = operator[site], operator[site + 1] H_AC2 = JordanMPO_AC2_Hamiltonian(GL, W1, W2, GR) - return prepare ? prepare_operator!!(H_AC2) : H_AC2 + return prepare ? prepare_operator!!(H_AC2, backend, allocator) : H_AC2 end for f in (:AC_hamiltonian, :AC2_hamiltonian) diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl index db78fe768..dceab226d 100644 --- a/src/algorithms/derivatives/mpo_derivatives.jl +++ b/src/algorithms/derivatives/mpo_derivatives.jl @@ -22,21 +22,30 @@ MPO_AC2_Hamiltonian(GL, O1, O2, GR) = MPODerivativeOperator(GL, (O1, O2), GR) # Constructors # ------------ -function C_hamiltonian(site::Int, below, operator, above, envs; prepare::Bool = true) +function C_hamiltonian( + site::Int, below, operator, above, envs; prepare::Bool = true, + backend::AbstractBackend = DefaultBackend(), allocator = BufferAllocator() + ) H_C = MPO_C_Hamiltonian(leftenv(envs, site + 1, below), rightenv(envs, site, below)) - return prepare ? prepare_operator!!(H_C) : H_C + return prepare ? prepare_operator!!(H_C, backend, allocator) : H_C end -function AC_hamiltonian(site::Int, below, operator, above, envs; prepare::Bool = true) +function AC_hamiltonian( + site::Int, below, operator, above, envs; prepare::Bool = true, + backend::AbstractBackend = DefaultBackend(), allocator = BufferAllocator() + ) O = isnothing(operator) ? nothing : operator[site] H_AC = MPO_AC_Hamiltonian(leftenv(envs, site, below), O, rightenv(envs, site, below)) - return prepare ? prepare_operator!!(H_AC) : H_AC + return prepare ? prepare_operator!!(H_AC, backend, allocator) : H_AC end -function AC2_hamiltonian(site::Int, below, operator, above, envs; prepare::Bool = true) +function AC2_hamiltonian( + site::Int, below, operator, above, envs; prepare::Bool = true, + backend::AbstractBackend = DefaultBackend(), allocator = BufferAllocator() + ) O1, O2 = isnothing(operator) ? (nothing, nothing) : (operator[site], operator[site + 1]) H_AC2 = MPO_AC2_Hamiltonian( leftenv(envs, site, below), O1, O2, rightenv(envs, site + 1, below) ) - return prepare ? prepare_operator!!(H_AC2) : H_AC2 + return prepare ? prepare_operator!!(H_AC2, backend, allocator) : H_AC2 end # Properties diff --git a/src/algorithms/derivatives/projection_derivatives.jl b/src/algorithms/derivatives/projection_derivatives.jl index b3a21fb39..52d70bdf3 100644 --- a/src/algorithms/derivatives/projection_derivatives.jl +++ b/src/algorithms/derivatives/projection_derivatives.jl @@ -14,18 +14,23 @@ Projection_AC2_Hamiltonian(GL, A1, A2, GR) = ProjectionDerivativeOperator(GL, (A # ------------ function AC_hamiltonian( site::Int, below, operator::ProjectionOperator, above, envs; - prepare::Bool = true + prepare::Bool = true, + backend::AbstractBackend = DefaultBackend(), allocator = BufferAllocator() ) GL = leftenv(envs, site, below) GR = rightenv(envs, site, below) H_AC = Projection_AC_Hamiltonian(GL, operator.ket.AC[site], GR) - return prepare ? prepare_operator!!(H_AC) : H_AC + return prepare ? prepare_operator!!(H_AC, backend, allocator) : H_AC end -function AC2_hamiltonian(site::Int, below, operator::ProjectionOperator, above, envs; prepare::Bool = true) +function AC2_hamiltonian( + site::Int, below, operator::ProjectionOperator, above, envs; + prepare::Bool = true, + backend::AbstractBackend = DefaultBackend(), allocator = BufferAllocator() + ) GL = leftenv(envs, site, below) GR = rightenv(envs, site + 1, below) H_AC2 = Projection_AC2_Hamiltonian(GL, operator.ket.AC[site], operator.ket.AR[site + 1], GR) - return prepare ? prepare_operator!!(H_AC2) : H_AC2 + return prepare ? prepare_operator!!(H_AC2, backend, allocator) : H_AC2 end # Actions diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index de8ac4bd8..d4a45c46f 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -24,8 +24,12 @@ given `trscheme`). ## Fields $(TYPEDFIELDS) + +!!! warning + The default allocator uses a buffer that is created at algorithm construction time. + As a result, it is **not thread-safe** to reuse the algorithm across multiple concurrent solves. """ -struct DMRG{A, F, E, G} <: Algorithm +struct DMRG{A, F, E, G, B, AL} <: Algorithm "tolerance for convergence criterium" tol::Float64 @@ -46,12 +50,19 @@ struct DMRG{A, F, E, G} <: Algorithm "factorization used for the post-update gauge: a QR algorithm (no truncation) or a truncated SVD" alg_gauge::G + + "backend for tensor contractions and index manipulations" + backend::B + + "allocator for handling (intermediate) tensor objects" + allocator::AL end function DMRG(; tol = Defaults.tol, maxiter = Defaults.maxiter, alg_eigsolve = (;), verbosity = Defaults.verbosity, finalize = Defaults._finalize, alg_expand = nothing, trscheme = notrunc(), - alg_svd = Defaults.alg_svd(), alg_orth = Defaults.alg_orth() + alg_svd = Defaults.alg_svd(), alg_orth = Defaults.alg_orth(), + backend = DefaultBackend(), allocator = BufferAllocator() ) # single-site DMRG defaults to the per-bond adaptive controller (`AdaptiveKrylov`); pass # `alg_eigsolve = (; adaptive = false, ...)` to opt out (the splat overrides the default). @@ -63,7 +74,7 @@ function DMRG(; if !isnothing(alg_expand) && !_truncates(alg_gauge) @warn "DMRG with `alg_expand` but no truncation (`trscheme = notrunc()`): the bond dimension will grow unboundedly each sweep." end - return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize, alg_expand, alg_gauge) + return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize, alg_expand, alg_gauge, backend, allocator) end @@ -83,7 +94,7 @@ function local_update!( alg_eigsolve = adapt_solver(alg.alg_eigsolve; decay_rate, g_local = ϵ_local, g_global = ϵ_global, eps_trunc = ϵ_trunc) ac_old = ψ.AC[site] λ, AC′, info = @timeit timeroutput "AC_eigsolve" begin - H_effective = AC_hamiltonian(site, ψ, O, ψ, envs) + H_effective = AC_hamiltonian(site, ψ, O, ψ, envs; alg.backend, alg.allocator) fixedpoint(H_effective, ac_old, :SR, alg_eigsolve) end @@ -106,8 +117,12 @@ Two-site DMRG algorithm for finding the dominant eigenvector. ## Fields $(TYPEDFIELDS) + +!!! warning + The default allocator uses a buffer that is created at algorithm construction time. + As a result, it is **not thread-safe** to reuse the algorithm across multiple concurrent solves. """ -struct DMRG2{A, G, F} <: Algorithm +struct DMRG2{A, G, F, B, AL} <: Algorithm "tolerance for convergence criterium" tol::Float64 @@ -125,12 +140,19 @@ struct DMRG2{A, G, F} <: Algorithm "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`" finalize::F + + "backend for tensor contractions and index manipulations" + backend::B + + "allocator for handling (intermediate) tensor objects" + allocator::AL end # TODO: find better default truncation function DMRG2(; tol = Defaults.tol, maxiter = Defaults.maxiter, verbosity = Defaults.verbosity, alg_eigsolve = (;), alg_svd = Defaults.alg_svd(), trscheme, - finalize = Defaults._finalize + finalize = Defaults._finalize, + backend = DefaultBackend(), allocator = BufferAllocator() ) # two-site DMRG defaults to the per-bond adaptive controller (`AdaptiveKrylov`); pass # `alg_eigsolve = (; adaptive = false, ...)` to opt out (the splat overrides the default). @@ -138,7 +160,7 @@ function DMRG2(; Defaults.alg_eigsolve(; adaptive = true, alg_eigsolve...) : alg_eigsolve # two-site DMRG always truncates the enlarged bond back down, so the gauge is a truncated SVD alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme) - return DMRG2(tol, maxiter, verbosity, alg_eigsolve′, alg_gauge, finalize) + return DMRG2(tol, maxiter, verbosity, alg_eigsolve′, alg_gauge, finalize, backend, allocator) end function local_update!( @@ -147,7 +169,7 @@ function local_update!( ϵ_global, ϵ_trunc, decay_rate, timeroutput ) - Heff = @timeit timeroutput "AC2_hamiltonian" AC2_hamiltonian(pos, ψ, O, ψ, envs) + Heff = @timeit timeroutput "AC2_hamiltonian" AC2_hamiltonian(pos, ψ, O, ψ, envs; alg.backend, alg.allocator) kind = direction === Val(:right) ? :ACAR : :ALAC ac2 = AC2(ψ, pos; kind) From 6ac438dd1b7bb45f16af54c1376f957056d9c39b Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 21 Jul 2026 09:02:02 -0400 Subject: [PATCH 2/3] use backends and allocators in more parts of JordanMPO --- .../derivatives/hamiltonian_derivatives.jl | 142 ++++++++++-------- 1 file changed, 81 insertions(+), 61 deletions(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index e9be54228..b35b71392 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -12,26 +12,32 @@ const _HAM_MPS_TYPES = Union{ Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. """ -struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator +struct JordanMPO_AC_Hamiltonian{O1, O2, O3, Bk, Al} <: DerivativeOperator D::Union{O1, Missing} # onsite I::Union{O1, Missing} # not started E::Union{O1, Missing} # finished C::Union{O2, Missing} # starting B::Union{O2, Missing} # ending A::Union{O3, Missing} # continuing + backend::Bk # contraction backend used by the matvec + allocator::Al # scratch-buffer allocator used by the matvec - function JordanMPO_AC_Hamiltonian{O1, O2, O3}( + function JordanMPO_AC_Hamiltonian{O1, O2, O3, Bk, Al}( D::Union{O1, Missing}, I::Union{O1, Missing}, E::Union{O1, Missing}, - C::Union{O2, Missing}, B::Union{O2, Missing}, A::Union{O3, Missing} - ) where {O1, O2, O3} - return new{O1, O2, O3}(D, I, E, C, B, A) + C::Union{O2, Missing}, B::Union{O2, Missing}, A::Union{O3, Missing}, + backend::Bk, allocator::Al + ) where {O1, O2, O3, Bk, Al} + return new{O1, O2, O3, Bk, Al}(D, I, E, C, B, A, backend, allocator) end end -function JordanMPO_AC_Hamiltonian{O1, O2, O3}(D, I, E, C, B, A) where {O1, O2, O3} - return JordanMPO_AC_Hamiltonian{O1, O2, O3}( +function JordanMPO_AC_Hamiltonian{O1, O2, O3}( + D, I, E, C, B, A, backend = DefaultBackend(), allocator = BufferAllocator() + ) where {O1, O2, O3} + return JordanMPO_AC_Hamiltonian{O1, O2, O3, typeof(backend), typeof(allocator)}( ismissing(D) ? D : convert(O1, D), ismissing(I) ? I : convert(O1, I), ismissing(E) ? E : convert(O1, E), ismissing(C) ? C : convert(O2, C), - ismissing(B) ? E : convert(O2, B), ismissing(A) ? A : convert(O3, A) + ismissing(B) ? E : convert(O2, B), ismissing(A) ? A : convert(O3, A), + backend, allocator ) end @@ -44,11 +50,14 @@ function AC_hamiltonian( GL = leftenv(envs, site, below) GR = rightenv(envs, site, below) W = operator[site] - H_AC = JordanMPO_AC_Hamiltonian(GL, W, GR) + H_AC = JordanMPO_AC_Hamiltonian(GL, W, GR; backend, allocator) return prepare ? prepare_operator!!(H_AC, backend, allocator) : H_AC end -function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTensor) +function JordanMPO_AC_Hamiltonian( + GL::MPSTensor, W::JordanMPOTensor, GR::MPSTensor; + backend::AbstractBackend = DefaultBackend(), allocator = BufferAllocator() + ) # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind # them once and reuse the locals throughout WA, WB, WC, WD = W.A, W.B, W.C, W.D @@ -65,7 +74,7 @@ function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTens # starting C = if nonzero_length(WC) > 0 GR_2 = GR[2:(end - 1)] - @plansor starting[-1 -2; -3 -4] ≔ WC[-1; -3 1] * GR_2[-4 1; -2] + @plansor backend = backend allocator = allocator starting[-1 -2; -3 -4] ≔ WC[-1; -3 1] * GR_2[-4 1; -2] only(starting) else missing @@ -74,7 +83,7 @@ function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTens # ending B = if nonzero_length(WB) > 0 GL_2 = GL[2:(end - 1)] - @plansor ending[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * WB[1 -2; -4] + @plansor backend = backend allocator = allocator ending[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * WB[1 -2; -4] only(ending) else missing @@ -94,7 +103,7 @@ function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTens # specialization for nearest neighbours nonzero_length(WA) == 0 && (A = missing) - return JordanMPO_AC_Hamiltonian{O1, O2, O3}(D, I, E, C, B, A) + return JordanMPO_AC_Hamiltonian{O1, O2, O3}(D, I, E, C, B, A, backend, allocator) end function prepare_operator!!( @@ -108,11 +117,11 @@ function prepare_operator!!( missing elseif !ismissing(C) Id = TensorKit.id(storagetype(C), space(C, 2)) - @plansor C[-1 -2; -3 -4] += H.D[-1; -3] * Id[-2; -4] + @plansor backend = backend allocator = allocator C[-1 -2; -3 -4] += H.D[-1; -3] * Id[-2; -4] missing elseif !ismissing(B) Id = TensorKit.id(storagetype(B), space(B, 1)) - @plansor B[-1 -2; -3 -4] += Id[-1; -3] * H.D[-2; -4] + @plansor backend = backend allocator = allocator B[-1 -2; -3 -4] += Id[-1; -3] * H.D[-2; -4] missing else H.D @@ -123,7 +132,7 @@ function prepare_operator!!( missing elseif !ismissing(C) Id = id(storagetype(C), space(C, 1)) - @plansor C[-1 -2; -3 -4] += Id[-1; -3] * H.I[-4; -2] + @plansor backend = backend allocator = allocator C[-1 -2; -3 -4] += Id[-1; -3] * H.I[-4; -2] missing else H.I @@ -134,7 +143,7 @@ function prepare_operator!!( missing elseif !ismissing(B) Id = id(storagetype(B), space(B, 2)) - @plansor B[-1 -2; -3 -4] += H.E[-1; -3] * Id[-2; -4] + @plansor backend = backend allocator = allocator B[-1 -2; -3 -4] += H.E[-1; -3] * Id[-2; -4] missing else H.E @@ -143,7 +152,7 @@ function prepare_operator!!( O3′ = prepared_operator_type(O3, typeof(backend), typeof(allocator)) A = ismissing(H.A) ? H.A : prepare_operator!!(H.A, backend, allocator) - return JordanMPO_AC_Hamiltonian{O1, O2, O3′}(D, I, E, C, B, A)::JordanMPO_AC_Hamiltonian{O1, O2, O3′} + return JordanMPO_AC_Hamiltonian{O1, O2, O3′}(D, I, E, C, B, A, backend, allocator)::JordanMPO_AC_Hamiltonian{O1, O2, O3′} end @@ -155,7 +164,7 @@ end Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs. In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor. """ -struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator +struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4, Bk, Al} <: DerivativeOperator II::Union{O1, Missing} # not_started IC::Union{O2, Missing} # starting right ID::Union{O1, Missing} # onsite right @@ -166,25 +175,30 @@ struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator BE::Union{O2, Missing} # ending left DE::Union{O1, Missing} # onsite left EE::Union{O1, Missing} # finished + backend::Bk # contraction backend used by the matvec + allocator::Al # scratch-buffer allocator used by the matvec - function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( + function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4, Bk, Al}( II::Union{O1, Missing}, IC::Union{O2, Missing}, ID::Union{O1, Missing}, CB::Union{O2, Missing}, CA::Union{O3, Missing}, AB::Union{O3, Missing}, AA::Union{O4, Missing}, - BE::Union{O2, Missing}, DE::Union{O1, Missing}, EE::Union{O1, Missing} - ) where {O1, O2, O3, O4} - return new{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + BE::Union{O2, Missing}, DE::Union{O1, Missing}, EE::Union{O1, Missing}, + backend::Bk, allocator::Al + ) where {O1, O2, O3, O4, Bk, Al} + return new{O1, O2, O3, O4, Bk, Al}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE, backend, allocator) end end function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( - II, IC, ID, CB, CA, AB, AA, BE, DE, EE + II, IC, ID, CB, CA, AB, AA, BE, DE, EE, + backend = DefaultBackend(), allocator = BufferAllocator() ) where {O1, O2, O3, O4} - return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( + return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4, typeof(backend), typeof(allocator)}( ismissing(II) ? II : convert(O1, II), ismissing(IC) ? IC : convert(O2, IC), ismissing(ID) ? ID : convert(O1, ID), ismissing(CB) ? CB : convert(O2, CB), ismissing(CA) ? CA : convert(O3, CA), ismissing(AB) ? AB : convert(O3, AB), ismissing(AA) ? AA : convert(O4, AA), ismissing(BE) ? BE : convert(O2, BE), - ismissing(DE) ? DE : convert(O1, DE), ismissing(EE) ? EE : convert(O1, EE) + ismissing(DE) ? DE : convert(O1, DE), ismissing(EE) ? EE : convert(O1, EE), + backend, allocator ) end @@ -197,7 +211,7 @@ function AC2_hamiltonian( GL = leftenv(envs, site, below) GR = rightenv(envs, site + 1, below) W1, W2 = operator[site], operator[site + 1] - H_AC2 = JordanMPO_AC2_Hamiltonian(GL, W1, W2, GR) + H_AC2 = JordanMPO_AC2_Hamiltonian(GL, W1, W2, GR; backend, allocator) return prepare ? prepare_operator!!(H_AC2, backend, allocator) : H_AC2 end @@ -210,7 +224,10 @@ for f in (:AC_hamiltonian, :AC2_hamiltonian) end end -function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::JordanMPOTensor, GR::MPSTensor) +function JordanMPO_AC2_Hamiltonian( + GL::MPSTensor, W1::JordanMPOTensor, W2::JordanMPOTensor, GR::MPSTensor; + backend::AbstractBackend = DefaultBackend(), allocator = BufferAllocator() + ) # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind # them once and reuse the locals throughout A1, B1, C1, D1 = W1.A, W1.B, W1.C, W1.D @@ -224,7 +241,7 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda # starting right IC = if nonzero_length(C2) > 0 - @plansor IC_[-1 -2; -3 -4] ≔ C2[-1; -3 1] * GR[2:(end - 1)][-4 1; -2] + @plansor backend = backend allocator = allocator IC_[-1 -2; -3 -4] ≔ C2[-1; -3 1] * GR[2:(end - 1)][-4 1; -2] only(IC_) else missing @@ -238,7 +255,7 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda # starting left - ending right CB = if nonzero_length(C1) > 0 && nonzero_length(B2) > 0 - @plansor CB_[-1 -2; -3 -4] ≔ C1[-1; -3 1] * B2[1 -2; -4] + @plansor backend = backend allocator = allocator CB_[-1 -2; -3 -4] ≔ C1[-1; -3 1] * B2[1 -2; -4] # have to convert to complex if hamiltonian is real but states are complex scalartype(GL) <: Complex ? complex(only(CB_)) : only(CB_) else @@ -247,7 +264,7 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda # starting left - continuing right CA = if nonzero_length(C1) > 0 && nonzero_length(A2) > 0 - @plansor CA_[-1 -2 -3; -4 -5 -6] ≔ C1[-1; -4 2] * A2[2 -2; -5 1] * + @plansor backend = backend allocator = allocator CA_[-1 -2 -3; -4 -5 -6] ≔ C1[-1; -4 2] * A2[2 -2; -5 1] * GR[2:(end - 1)][-6 1; -3] only(CA_) else @@ -256,7 +273,7 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda # continuing left - ending right AB = if nonzero_length(A1) > 0 && nonzero_length(B2) > 0 - @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL[2:(end - 1)][-1 2; -4] * A1[2 -2; -5 1] * + @plansor backend = backend allocator = allocator AB_[-1 -2 -3; -4 -5 -6] ≔ GL[2:(end - 1)][-1 2; -4] * A1[2 -2; -5 1] * B2[1 -3; -6] only(AB_) else @@ -265,7 +282,7 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda # ending left BE = if nonzero_length(B1) > 0 - @plansor BE_[-1 -2; -3 -4] ≔ GL[2:(end - 1)][-1 2; -3] * B1[2 -2; -4] + @plansor backend = backend allocator = allocator BE_[-1 -2; -3 -4] ≔ GL[2:(end - 1)][-1 2; -3] * B1[2 -2; -4] only(BE_) else missing @@ -303,7 +320,8 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda II, IC, ID, CB, CA, AB, AA, - BE, DE, EE + BE, DE, EE, + backend, allocator ) end @@ -317,11 +335,11 @@ function prepare_operator!!( CB::Union{Missing, O2} = if !ismissing(CA) && !ismissing(H.CB) Id = TensorKit.id(storagetype(H.CB), space(CA, 3)) - @plansor CA[-1 -2 -3; -4 -5 -6] += H.CB[-1 -2; -4 -5] * Id[-3; -6] + @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += H.CB[-1 -2; -4 -5] * Id[-3; -6] missing elseif !ismissing(AB) && !ismissing(H.CB) Id = TensorKit.id(storagetype(H.CB), space(AB, 1)) - @plansor AB[-1 -2 -3; -4 -5 -6] += H.CB[-2 -3; -5 -6] * Id[-1; -4] + @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += H.CB[-2 -3; -5 -6] * Id[-1; -4] missing else H.CB @@ -330,7 +348,7 @@ function prepare_operator!!( # starting right IC::Union{Missing, O2} = if !ismissing(CA) && !ismissing(H.IC) Id = TensorKit.id(storagetype(H.IC), space(CA, 1)) - @plansor CA[-1 -2 -3; -4 -5 -6] += Id[-1; -4] * H.IC[ -2 -3; -5 -6] + @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += Id[-1; -4] * H.IC[ -2 -3; -5 -6] missing else H.IC @@ -339,7 +357,7 @@ function prepare_operator!!( # ending left BE::Union{Missing, O2} = if !ismissing(AB) && !ismissing(H.BE) Id = TensorKit.id(storagetype(H.BE), space(AB, 3)) - @plansor AB[-1 -2 -3; -4 -5 -6] += H.BE[-1 -2; -4 -5] * Id[-3; -6] + @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += H.BE[-1 -2; -4 -5] * Id[-3; -6] missing else H.BE @@ -348,12 +366,12 @@ function prepare_operator!!( # onsite left DE::Union{Missing, O1} = if !ismissing(BE) && !ismissing(H.DE) Id = TensorKit.id(storagetype(H.DE), space(BE, 1)) - @plansor BE[-1 -2; -3 -4] += Id[-1; -3] * H.DE[-2; -4] + @plansor backend = backend allocator = allocator BE[-1 -2; -3 -4] += Id[-1; -3] * H.DE[-2; -4] missing elseif !ismissing(AB) && !ismissing(H.DE) Id1 = id(storagetype(H.DE), space(AB, 1)) Id2 = id(storagetype(H.DE), space(AB, 3)) - @plansor AB[-1 -2 -3; -4 -5 -6] += Id1[-1; -4] * H.DE[-2; -5] * Id2[-3; -6] + @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += Id1[-1; -4] * H.DE[-2; -5] * Id2[-3; -6] missing # TODO: could also try in CA? else @@ -363,12 +381,12 @@ function prepare_operator!!( # onsite right ID::Union{Missing, O1} = if !ismissing(IC) && !ismissing(H.ID) Id = TensorKit.id(storagetype(H.ID), space(IC, 2)) - @plansor IC[-1 -2; -3 -4] += H.ID[-1; -3] * Id[-2; -4] + @plansor backend = backend allocator = allocator IC[-1 -2; -3 -4] += H.ID[-1; -3] * Id[-2; -4] missing elseif !ismissing(CA) && !ismissing(H.ID) Id1 = TensorKit.id(storagetype(H.ID), space(CA, 1)) Id2 = TensorKit.id(storagetype(H.ID), space(CA, 3)) - @plansor CA[-1 -2 -3; -4 -5 -6] += Id1[-1; -4] * H.ID[-2; -5] * Id2[-3; -6] + @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += Id1[-1; -4] * H.ID[-2; -5] * Id2[-3; -6] missing else H.ID @@ -377,11 +395,11 @@ function prepare_operator!!( # finished II::Union{Missing, O1} = if !ismissing(IC) && !ismissing(H.II) I = id(storagetype(H.II), space(IC, 1)) - @plansor IC[-1 -2; -3 -4] += I[-1; -3] * H.II[-2; -4] + @plansor backend = backend allocator = allocator IC[-1 -2; -3 -4] += I[-1; -3] * H.II[-2; -4] II = missing elseif !ismissing(CA) && !ismissing(H.II) I = id(storagetype(H.II), space(CA, 1) ⊗ space(CA, 2)) - @plansor CA[-1 -2 -3; -4 -5 -6] += I[-1 -2; -4 -5] * H.II[-3; -6] + @plansor backend = backend allocator = allocator CA[-1 -2 -3; -4 -5 -6] += I[-1 -2; -4 -5] * H.II[-3; -6] II = missing else H.II @@ -390,11 +408,11 @@ function prepare_operator!!( # unstarted EE::Union{Missing, O1} = if !ismissing(BE) && !ismissing(H.EE) I = id(storagetype(H.EE), space(BE, 2)) - @plansor BE[-1 -2; -3 -4] += H.EE[-1; -3] * I[-2; -4] + @plansor backend = backend allocator = allocator BE[-1 -2; -3 -4] += H.EE[-1; -3] * I[-2; -4] EE = missing elseif !ismissing(AB) && !ismissing(H.EE) I = id(storagetype(H.EE), space(AB, 2) ⊗ space(AB, 3)) - @plansor AB[-1 -2 -3; -4 -5 -6] += H.EE[-1; -4] * I[-2 -3; -5 -6] + @plansor backend = backend allocator = allocator AB[-1 -2 -3; -4 -5 -6] += H.EE[-1; -4] * I[-2 -3; -5 -6] EE = missing else H.EE @@ -403,35 +421,37 @@ function prepare_operator!!( O4′ = prepared_operator_type(O4, typeof(backend), typeof(allocator)) AA = prepare_operator!!(H.AA, backend, allocator) - return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4′}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4′}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE, backend, allocator) end # Actions # ------- function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) + backend, allocator = H.backend, H.allocator y = ismissing(H.A) ? zerovector(x) : H.A(x) - ismissing(H.D) || @plansor y[-1 -2; -3] += x[-1 1; -3] * H.D[-2; 1] - ismissing(H.E) || @plansor y[-1 -2; -3] += H.E[-1; 1] * x[1 -2; -3] - ismissing(H.I) || @plansor y[-1 -2; -3] += x[-1 -2; 1] * H.I[1; -3] - ismissing(H.C) || @plansor y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1] - ismissing(H.B) || @plansor y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3] + ismissing(H.D) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += x[-1 1; -3] * H.D[-2; 1] + ismissing(H.E) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += H.E[-1; 1] * x[1 -2; -3] + ismissing(H.I) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += x[-1 -2; 1] * H.I[1; -3] + ismissing(H.C) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1] + ismissing(H.B) || @plansor backend = backend allocator = allocator y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3] return y end function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) + backend, allocator = H.backend, H.allocator y = ismissing(H.AA) ? zerovector(x) : H.AA(x) - ismissing(H.II) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 -4] * H.II[-3; 1] - ismissing(H.IC) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 2] * H.IC[-4 -3; 2 1] - ismissing(H.ID) || @plansor y[-1 -2; -3 -4] += x[-1 -2; -3 1] * H.ID[-4; 1] - ismissing(H.CB) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 2] * H.CB[-2 -4; 1 2] - ismissing(H.CA) || @plansor y[-1 -2; -3 -4] += x[-1 1; 3 2] * H.CA[-2 -4 -3; 1 2 3] - ismissing(H.AB) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 3] * H.AB[-1 -2 -4; 1 2 3] - ismissing(H.BE) || @plansor y[-1 -2; -3 -4] += x[1 2; -3 -4] * H.BE[-1 -2; 1 2] - ismissing(H.DE) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1] - ismissing(H.EE) || @plansor y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1] + ismissing(H.II) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 -2; 1 -4] * H.II[-3; 1] + ismissing(H.IC) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 -2; 1 2] * H.IC[-4 -3; 2 1] + ismissing(H.ID) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 -2; -3 1] * H.ID[-4; 1] + ismissing(H.CB) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 1; -3 2] * H.CB[-2 -4; 1 2] + ismissing(H.CA) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 1; 3 2] * H.CA[-2 -4 -3; 1 2 3] + ismissing(H.AB) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[1 2; -3 3] * H.AB[-1 -2 -4; 1 2 3] + ismissing(H.BE) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[1 2; -3 -4] * H.BE[-1 -2; 1 2] + ismissing(H.DE) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1] + ismissing(H.EE) || @plansor backend = backend allocator = allocator y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1] return y end From d1138d02079c0072557363889862fbf8180864d8 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 21 Jul 2026 10:00:13 -0400 Subject: [PATCH 3/3] fix typo --- src/algorithms/derivatives/hamiltonian_derivatives.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index b35b71392..0aa520aa2 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -36,7 +36,7 @@ function JordanMPO_AC_Hamiltonian{O1, O2, O3}( return JordanMPO_AC_Hamiltonian{O1, O2, O3, typeof(backend), typeof(allocator)}( ismissing(D) ? D : convert(O1, D), ismissing(I) ? I : convert(O1, I), ismissing(E) ? E : convert(O1, E), ismissing(C) ? C : convert(O2, C), - ismissing(B) ? E : convert(O2, B), ismissing(A) ? A : convert(O3, A), + ismissing(B) ? B : convert(O2, B), ismissing(A) ? A : convert(O3, A), backend, allocator ) end