diff --git a/ext/MPSKitAdaptExt.jl b/ext/MPSKitAdaptExt.jl index 78f4b41af..912b46569 100644 --- a/ext/MPSKitAdaptExt.jl +++ b/ext/MPSKitAdaptExt.jl @@ -1,6 +1,6 @@ module MPSKitAdaptExt -using TensorKit: space, spacetype +using TensorKit: space, spacetype, scalartype using MPSKit using BlockTensorKit: nonzero_pairs using Adapt @@ -31,8 +31,12 @@ end # inline to improve type stability with closures @inline Adapt.adapt_structure(to, mpo::MPO) = MPO(map(adapt(to), mpo.O)) -@inline Adapt.adapt_structure(to, W::MPSKit.JordanMPOTensor) = - MPSKit.JordanMPOTensor(space(W), adapt(to, W.A), adapt(to, W.B), adapt(to, W.C), adapt(to, W.D)) +@inline function Adapt.adapt_structure(to, W::MPSKit.JordanMPOTensor) + tensors = adapt(to, W.tensors) + T = scalartype(tensors) + scalars = Dict{CartesianIndex{4}, T}(K => convert(T, c) for (K, c) in W.scalars) + return MPSKit.JordanMPOTensor(tensors, scalars) +end @inline Adapt.adapt_structure(to, mpo::MPOHamiltonian) = MPOHamiltonian(map(x -> adapt(to, x), mpo.W)) @inline Adapt.adapt_structure(to, ml::MPSKit.Multiline) = MPSKit.Multiline(map(adapt(to), ml.data)) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index c39977731..a4e1fe605 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -18,7 +18,7 @@ export QP, LeftGaugedQP, RightGaugedQP # operators: export AbstractMPO export MPO, FiniteMPO, InfiniteMPO -export JordanMPOTensor, JordanMPOTensorMap +export JordanMPOTensor export MPOHamiltonian, FiniteMPOHamiltonian, InfiniteMPOHamiltonian export MultilineMPO export UntimedOperator, TimedOperator, MultipliedOperator, LazySum diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index fb93d67b6..47ee2b25a 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -48,8 +48,12 @@ function AC_hamiltonian( end function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTensor) + # 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 + # onsite - D = nonzero_length(W.D) > 0 ? only(W.D) : missing + D = nonzero_length(WD) > 0 ? only(WD) : missing # not started I = size(W, 4) == 1 ? missing : removeunit(GR[1], 2) @@ -58,25 +62,25 @@ function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTens E = size(W, 1) == 1 ? missing : removeunit(GL[end], 2) # starting - C = if nonzero_length(W.C) > 0 + C = if nonzero_length(WC) > 0 GR_2 = GR[2:(end - 1)] - @plansor starting[-1 -2; -3 -4] ≔ W.C[-1; -3 1] * GR_2[-4 1; -2] + @plansor starting[-1 -2; -3 -4] ≔ WC[-1; -3 1] * GR_2[-4 1; -2] only(starting) else missing end # ending - B = if nonzero_length(W.B) > 0 + B = if nonzero_length(WB) > 0 GL_2 = GL[2:(end - 1)] - @plansor ending[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * W.B[1 -2; -4] + @plansor ending[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * WB[1 -2; -4] only(ending) else missing end # continuing - A = MPO_AC_Hamiltonian(GL[2:(end - 1)], W.A, GR[2:(end - 1)]) + A = MPO_AC_Hamiltonian(GL[2:(end - 1)], WA, GR[2:(end - 1)]) # obtaining storagetype of environments since these should have already mixed # the types of the operator and state @@ -87,7 +91,7 @@ function JordanMPO_AC_Hamiltonian(GL::MPSTensor, W::JordanMPOTensor, GR::MPSTens O3 = typeof(A) # specialization for nearest neighbours - nonzero_length(W.A) == 0 && (A = missing) + nonzero_length(WA) == 0 && (A = missing) return JordanMPO_AC_Hamiltonian{O1, O2, O3}(D, I, E, C, B, A) end @@ -196,6 +200,11 @@ function AC2_hamiltonian( end function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::JordanMPOTensor, GR::MPSTensor) + # 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 + A2, B2, C2, D2 = W2.A, W2.B, W2.C, W2.D + # not started II = size(W2, 4) == 1 ? missing : transpose(removeunit(GR[1], 2)) @@ -203,22 +212,22 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda EE = size(W1, 1) == 1 ? missing : removeunit(GL[end], 2) # starting right - IC = if nonzero_length(W2.C) > 0 - @plansor IC_[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR[2:(end - 1)][-4 1; -2] + IC = if nonzero_length(C2) > 0 + @plansor IC_[-1 -2; -3 -4] ≔ C2[-1; -3 1] * GR[2:(end - 1)][-4 1; -2] only(IC_) else missing end # onsite left - DE = nonzero_length(W1.D) > 0 ? only(W1.D) : missing + DE = nonzero_length(D1) > 0 ? only(D1) : missing # onsite right - ID = nonzero_length(W2.D) > 0 ? only(W2.D) : missing + ID = nonzero_length(D2) > 0 ? only(D2) : missing # starting left - ending right - CB = if nonzero_length(W1.C) > 0 && nonzero_length(W2.B) > 0 - @plansor CB_[-1 -2; -3 -4] ≔ W1.C[-1; -3 1] * W2.B[1 -2; -4] + CB = if nonzero_length(C1) > 0 && nonzero_length(B2) > 0 + @plansor 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 @@ -226,8 +235,8 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda end # starting left - continuing right - CA = if nonzero_length(W1.C) > 0 && nonzero_length(W2.A) > 0 - @plansor CA_[-1 -2 -3; -4 -5 -6] ≔ W1.C[-1; -4 2] * W2.A[2 -2; -5 1] * + 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] * GR[2:(end - 1)][-6 1; -3] only(CA_) else @@ -235,24 +244,24 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda end # continuing left - ending right - AB = if nonzero_length(W1.A) > 0 && nonzero_length(W2.B) > 0 - @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL[2:(end - 1)][-1 2; -4] * W1.A[2 -2; -5 1] * - W2.B[1 -3; -6] + 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] * + B2[1 -3; -6] only(AB_) else missing end # ending left - BE = if nonzero_length(W1.B) > 0 - @plansor BE_[-1 -2; -3 -4] ≔ GL[2:(end - 1)][-1 2; -3] * W1.B[2 -2; -4] + BE = if nonzero_length(B1) > 0 + @plansor BE_[-1 -2; -3 -4] ≔ GL[2:(end - 1)][-1 2; -3] * B1[2 -2; -4] only(BE_) else missing end # continuing - continuing - AA = MPO_AC2_Hamiltonian(GL[2:(end - 1)], W1.A, W2.A, GR[2:(end - 1)]) + AA = MPO_AC2_Hamiltonian(GL[2:(end - 1)], A1, A2, GR[2:(end - 1)]) S = spacetype(GL) M = storagetype(GL) @@ -261,16 +270,16 @@ function JordanMPO_AC2_Hamiltonian(GL::MPSTensor, W1::JordanMPOTensor, W2::Jorda O3 = tensormaptype(S, 3, 3, M) O4 = typeof(AA) - if nonzero_length(W1.A) == 0 && nonzero_length(W2.A) == 0 + if nonzero_length(A1) == 0 && nonzero_length(A2) == 0 AA = missing else - mask1 = falses(size(W1.A, 1), size(W1.A, 4)) - for I in nonzero_keys(W1.A) + mask1 = falses(size(A1, 1), size(A1, 4)) + for I in nonzero_keys(A1) mask1[I[1], I[4]] = true end - mask2 = falses(size(W2.A, 1), size(W2.A, 4)) - for I in nonzero_keys(W2.A) + mask2 = falses(size(A2, 1), size(A2, 4)) + for I in nonzero_keys(A2) mask2[I[1], I[4]] = true end diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl index af95282cc..8dd20f6fb 100644 --- a/src/operators/jordanmpotensor.jl +++ b/src/operators/jordanmpotensor.jl @@ -1,70 +1,77 @@ """ - JordanMPOTensor{E,S,TA,TB,TC,TD} <: AbstractBlockTensorMap{E,S,2,2} + JordanMPOTensor{T,S,A} <: AbstractBlockTensorMap{T,S,2,2} -A tensor map that represents the upper triangular block form of a matrix product operator (MPO). +A single tensor of a matrix product operator (MPO) in upper triangular (Jordan) block form, +as used to represent the local tensors of an [`MPOHamiltonian`](@ref). +The virtual (row, column) structure is ```math \\begin{pmatrix} 1 & C & D \\\\ -0 & A & B \\\\ -0 & 0 & 1 +⋅ & A & B \\\\ +⋅ & ⋅ & 1 \\end{pmatrix} ``` -""" -struct JordanMPOTensor{ - E, S, - TA <: AbstractTensorMap{E, S, 2, 2}, - TB <: AbstractTensorMap{E, S, 2, 1}, - TC <: AbstractTensorMap{E, S, 1, 2}, - TD <: AbstractTensorMap{E, S, 1, 1}, - } <: AbstractBlockTensorMap{E, S, 2, 2} - V::TensorMapSumSpace{S, 2, 2} - A::SparseBlockTensorMap{TA, E, S, 2, 2, 4} - B::SparseBlockTensorMap{TB, E, S, 2, 1, 3} - C::SparseBlockTensorMap{TC, E, S, 1, 2, 3} - D::SparseBlockTensorMap{TD, E, S, 1, 1, 2} - # uninitialized constructor - function JordanMPOTensor{E, S, TA, TB, TC, TD}( - ::UndefInitializer, V::TensorMapSumSpace{S, 2, 2} - ) where {E, S, TA, TB, TC, TD} - allVs = eachspace(V) - # Note that this is a bit of a hack using end to get the last index: - # it should be 1 or end depending on this being an "edge" tensor or a "bulk" tensor - VA = space(allVs[2:(end - 1), 1, 1, 2:(end - 1)]) - A = SparseBlockTensorMap{TA}(undef, VA) +where `A` is the bulk of interacting operators, `C`/`B` are the operators that start/finish +an interaction, `D` is the on-site term, and the diagonal `1`s are identities. + +## Representation - VB = removeunit(space(allVs[2:(end - 1), 1, 1, end]), 4) - B = SparseBlockTensorMap{TB}(undef, VB) +Rather than storing the dense block matrix, the genuine operators and the identities are kept separately: - VC = removeunit(space(allVs[1, 1, 1, 2:(end - 1)]), 1) - C = SparseBlockTensorMap{TC}(undef, VC) +- `tensors::SparseBlockTensorMap` holds the non-identity operators over the *full* virtual + space (so `A`, `B`, `C` and `D` all live at their `(row, 1, 1, col)` position). +- `scalars::Dict{CartesianIndex{4},T}` holds the scalar multiples of the identity, keyed by + their `(row, 1, 1, col)` virtual position; the diagonal corner `1`s are stored here as well. - VD = removeunit(removeunit(space(allVs[1, 1, 1, end:end]), 4), 1) - D = SparseBlockTensorMap{TD}(undef, VD) +An index is never present in both `tensors` and `scalars`. +This keeps the ubiquitous identity blocks free of dense storage and lets identities be materialized lazily only when needed. - return new{E, S, TA, TB, TC, TD}(V, A, B, C, D) +## Type parameters + +- `T <: Number`: the `scalartype` of the tensors. +- `S`: the `spacetype` of the tensors. +- `A <: DenseVector{T}`: the storage type of the underlying tensors. + +## Block accessors + +The reduced-leg `A`, `B`, `C` and `D` blocks are exposed as properties (`W.A`, `W.B`, `W.C`, `W.D`), +reconstructed on demand from `tensors` and `scalars`. +""" +struct JordanMPOTensor{ + T <: Number, S, A <: DenseVector{T}, + } <: AbstractBlockTensorMap{T, S, 2, 2} + tensors::SparseBlockTensorMap{TensorMap{T, S, 2, 2, A}, T, S, 2, 2, 4} + scalars::Dict{CartesianIndex{4}, T} + + # constructor from fields + function JordanMPOTensor{T, S, A}( + tensors::SparseBlockTensorMap{TensorMap{T, S, 2, 2, A}, T, S, 2, 2, 4}, + scalars::Dict{CartesianIndex{4}, T} + ) where {T, S, A} + return new{T, S, A}(tensors, scalars) end - # constructor from data - function JordanMPOTensor{E, S, TA, TB, TC, TD}( - V::TensorMapSumSpace, - A::SparseBlockTensorMap{<:Any, <:Any, S, 2, 2}, - B::SparseBlockTensorMap{<:Any, <:Any, S, 2, 1}, - C::SparseBlockTensorMap{<:Any, <:Any, S, 1, 2}, - D::SparseBlockTensorMap{<:Any, <:Any, S, 1, 1} - ) where {E, S, TA, TB, TC, TD} - return new{E, S, TA, TB, TC, TD}(V, A, B, C, D) + # uninitialized constructor + function JordanMPOTensor{T, S, A}( + ::UndefInitializer, V::TensorMapSumSpace{S, 2, 2} + ) where {T, S, A} + tensors = SparseBlockTensorMap{tensormaptype(S, 2, 2, A)}(undef, V) + scalars = Dict{CartesianIndex{4}, T}() + rows, cols = size(tensors, 1), size(tensors, 4) + cols > 1 && (scalars[CartesianIndex(1, 1, 1, 1)] = one(T)) + rows > 1 && (scalars[CartesianIndex(rows, 1, 1, cols)] = one(T)) + return new{T, S, A}(tensors, scalars) end end -const JordanMPOTensorMap{T, S, A <: DenseVector{T}} = JordanMPOTensor{ - T, S, - Union{TensorMap{T, S, 2, 2, A}, BraidingTensor{T, S, A}}, - TensorMap{T, S, 2, 1, A}, - TensorMap{T, S, 1, 2, A}, - TensorMap{T, S, 1, 1, A}, -} +function JordanMPOTensor( + tensors::SparseBlockTensorMap{TensorMap{T, S, 2, 2, A}, T, S, 2, 2, 4}, + scalars::Dict{CartesianIndex{4}, T} + ) where {T, S, A} + return JordanMPOTensor{T, S, A}(tensors, scalars) +end function JordanMPOTensor{E, S}(::UndefInitializer, V::TensorMapSumSpace{S}) where {E, S} return jordanmpotensortype(S, E)(undef, V) @@ -91,7 +98,27 @@ function JordanMPOTensor( VD = removeunit(removeunit(space(allVs[1, 1, 1, end:end]), 4), 1) VD == space(D) || throw(SpaceMismatch(lazy"D-block has incompatible spaces:\n$VD\n$(space(D))")) - return JordanMPOTensor{E, S, TA, TB, TC, TD}(V, A, B, C, D) + W = jordanmpotensortype(S, storagetype(A))(undef, V) + cols = size(W, 4) + # the trivial boundary legs inherit their duality from the virtual spaces of V + corner = allVs[1, 1, 1, end] + dual_left = isdual(corner[1]) + dual_right = !isdual(corner[4]) + for (I, v) in nonzero_pairs(A) + W[I[1] + 1, 1, 1, I[4] + 1] = v + end + for (I, v) in nonzero_pairs(B) + W[I[1] + 1, 1, 1, cols] = insertrightunit(v, 3; dual = dual_right) + end + for (I, v) in nonzero_pairs(C) + W[1, 1, 1, I[3] + 1] = insertleftunit(v, 1; dual = dual_left) + end + if nonzero_length(D) > 0 + W[1, 1, 1, cols] = insertrightunit( + insertleftunit(only(D), 1; dual = dual_left), 3; dual = dual_right + ) + end + return W end function JordanMPOTensor( V::TensorMapSumSpace{S, 2, 2}, @@ -107,28 +134,29 @@ function JordanMPOTensor( ) end +function _isbraiding(t::AbstractTensorMap) + t isa BraidingTensor && return true + τ = BraidingTensor{scalartype(t), spacetype(t), storagetype(t)}(space(t)) + return isapprox(t, τ) +end + function JordanMPOTensor(W::SparseBlockTensorMap{TT, E, S, 2, 2}) where {TT, E, S} - @assert W[1, 1, 1, 1] isa BraidingTensor && W[end, 1, 1, end] isa BraidingTensor + # the diagonal corners must be (close to) identity braidings + @assert _isbraiding(W[1, 1, 1, 1]) && _isbraiding(W[end, 1, 1, end]) # @assert all(I -> I[1] ≤ I[4], nonzero_keys(W)) - A = W[2:(end - 1), 1, 1, 2:(end - 1)] - B = W[2:(end - 1), 1, 1, end] - C = W[1, 1, 1, 2:(end - 1)] - D = W[1, 1, 1, end:end] # ensure still blocktensor to allow for sparse - - return JordanMPOTensor( - space(W), A, removeunit(B, 4), removeunit(C, 1), removeunit(removeunit(D, 4), 1) - ) + O = jordanmpotensortype(S, storagetype(W))(undef, space(W)) + corners = (CartesianIndex(1, 1, 1, 1), CartesianIndex(size(W, 1), 1, 1, size(W, 4))) + for (I, v) in nonzero_pairs(W) + I in corners && continue # already set as identity scalars by the constructor + O[I] = v + end + return O end function jordanmpotensortype(::Type{S}, ::Type{E}) where {S <: VectorSpace, E} - TA = tensormaptype(S, 2, 2, E) - T = scalartype(TA) - Tτ = BraidingTensor{T, S, storagetype(TA)} - TB = tensormaptype(S, 2, 1, E) - TC = tensormaptype(S, 1, 2, E) - TD = tensormaptype(S, 1, 1, E) - return JordanMPOTensor{T, S, Union{TA, Tτ}, TB, TC, TD} + TT = tensormaptype(S, 2, 2, E) + return JordanMPOTensor{scalartype(TT), S, storagetype(TT)} end function jordanmpotensortype(::Type{O}) where {O <: AbstractTensorMap} return jordanmpotensortype(spacetype(O), storagetype(O)) @@ -140,31 +168,86 @@ end # Properties # ---------- -TensorKit.space(W::JordanMPOTensor) = W.V -Base.eltype(::Type{JordanMPOTensor{E, S, TA, TB, TC, TD}}) where {E, S, TA, TB, TC, TD} = TA +TensorKit.space(W::JordanMPOTensor) = space(W.tensors) +function Base.eltype(::Type{<:JordanMPOTensor{T, S, A}}) where {T, S, A} + return tensormaptype(S, 2, 2, A) +end + +function Base.getproperty(W::JordanMPOTensor, sym::Symbol) + sym === :A && return _jordan_A(W) + sym === :B && return _jordan_B(W) + sym === :C && return _jordan_C(W) + sym === :D && return _jordan_D(W) + return getfield(W, sym) +end + +# materialize an identity scalar at virtual position `K` as a dense `TensorMap` +@inline function _identity_tensor(W::JordanMPOTensor, K::CartesianIndex{4}, c = one(scalartype(W))) + τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(eachspace(W)[K]) + return isone(c) ? TensorMap(τ) : scale!(TensorMap(τ), c) +end + +# reduced-leg blocks reconstructed from `tensors` and `scalars` +function _jordan_A(W::JordanMPOTensor) + rows, cols = size(W, 1), size(W, 4) + TA = tensormaptype(spacetype(W), 2, 2, storagetype(W)) + VA = space(eachspace(W.tensors)[2:(end - 1), 1, 1, 2:(end - 1)]) + A = SparseBlockTensorMap{TA}(undef, VA) + for (I, v) in nonzero_pairs(W.tensors) + (1 < I[1] < rows && 1 < I[4] < cols) || continue + A[I[1] - 1, 1, 1, I[4] - 1] = v + end + for (K, c) in W.scalars + (1 < K[1] < rows && 1 < K[4] < cols) || continue + A[K[1] - 1, 1, 1, K[4] - 1] = _identity_tensor(W, K, c) + end + return A +end +function _jordan_B(W::JordanMPOTensor) + rows, cols = size(W, 1), size(W, 4) + TB = tensormaptype(spacetype(W), 2, 1, storagetype(W)) + VB = removeunit(space(eachspace(W.tensors)[2:(end - 1), 1, 1, end]), 4) + B = SparseBlockTensorMap{TB}(undef, VB) + for (I, v) in nonzero_pairs(W.tensors) + (1 < I[1] < rows && I[4] == cols) || continue + B[I[1] - 1, 1, 1] = removeunit(v, 4) + end + for (K, c) in W.scalars + (1 < K[1] < rows && K[4] == cols) || continue + B[K[1] - 1, 1, 1] = removeunit(_identity_tensor(W, K, c), 4) + end + return B +end +function _jordan_C(W::JordanMPOTensor) + cols = size(W, 4) + TC = tensormaptype(spacetype(W), 1, 2, storagetype(W)) + VC = removeunit(space(eachspace(W.tensors)[1, 1, 1, 2:(end - 1)]), 1) + C = SparseBlockTensorMap{TC}(undef, VC) + for (I, v) in nonzero_pairs(W.tensors) + (I[1] == 1 && 1 < I[4] < cols) || continue + C[1, 1, I[4] - 1] = removeunit(v, 1) + end + for (K, c) in W.scalars + (K[1] == 1 && 1 < K[4] < cols) || continue + C[1, 1, K[4] - 1] = removeunit(_identity_tensor(W, K, c), 1) + end + return C +end +function _jordan_D(W::JordanMPOTensor) + cols = size(W, 4) + TD = tensormaptype(spacetype(W), 1, 1, storagetype(W)) + VD = removeunit(removeunit(space(eachspace(W.tensors)[1, 1, 1, end:end]), 4), 1) + D = SparseBlockTensorMap{TD}(undef, VD) + K = CartesianIndex(1, 1, 1, cols) + if haskey(W, K) + D[1, 1] = removeunit(removeunit(W[K], 4), 1) + end + return D +end function Base.haskey(W::JordanMPOTensor, I::CartesianIndex{4}) Base.checkbounds(W, I.I...) - # only has braiding tensors if sizes are large enough - sz = size(W) - ( - sz[1] > 1 && I == CartesianIndex(1, 1, 1, 1) || - sz[4] > 1 && I == CartesianIndex(sz[1], 1, 1, sz[4]) - ) && return true - - row, col = I.I[1], I.I[4] - - if row == 1 && col == sz[4] - return haskey(W.D, CartesianIndex(1, 1)) - elseif row == 1 - return haskey(W.C, CartesianIndex(1, 1, col - 1)) - elseif col == sz[4] - return haskey(W.B, CartesianIndex(row - 1, 1, 1)) - elseif 1 < row < sz[1] && 1 < col < sz[4] - return haskey(W.A, CartesianIndex(row - 1, 1, 1, col - 1)) - else - return false - end + return haskey(W.scalars, I) || haskey(W.tensors, I) end Base.parent(W::JordanMPOTensor) = parent(SparseBlockTensorMap(W)) @@ -174,51 +257,25 @@ BlockTensorKit.issparse(W::JordanMPOTensor) = true # Converters # ---------- function BlockTensorKit.SparseBlockTensorMap(W::JordanMPOTensor) - τ = BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(eachspace(W)[1]) - W′ = SparseBlockTensorMap{AbstractTensorMap{scalartype(W), spacetype(W), 2, 2}}( - undef_blocks, space(W) - ) - if size(W, 4) > 1 - W′[1, 1, 1, 1] = τ - end - if size(W, 1) > 1 - W′[end, 1, 1, end] = τ - end - - Ia = CartesianIndex(1, 0, 0, 1) - for (I, v) in nonzero_pairs(W.A) - W′[I + Ia] = v + W′ = SparseBlockTensorMap{eltype(W)}(undef_blocks, space(W)) + for (I, v) in nonzero_pairs(W.tensors) + W′[I] = v end - Ib = CartesianIndex(1, 0, 0) - for (I, v) in nonzero_pairs(W.B) - W′[I + Ib, size(W′, 4)] = insertrightunit(v, 3) + for (K, c) in W.scalars + W′[K] = _identity_tensor(W, K, c) end - Ic = CartesianIndex(0, 0, 1) - for (I, v) in nonzero_pairs(W.C) - W′[1, I + Ic] = insertleftunit(v, 1) - end - if nonzero_length(W.D) > 0 - W′[1, 1, 1, end] = insertrightunit(insertleftunit(only(W.D), 1), 3) - end - return W′ end for f in (:real, :complex) @eval function Base.$f(W::JordanMPOTensor) - E = $f(scalartype(W)) - W′ = similar(W, E) - for (I, v) in nonzero_pairs(W.A) - W′.A[I] = $f(v) - end - for (I, v) in nonzero_pairs(W.B) - W′.B[I] = $f(v) + W′ = similar(W, $f(scalartype(W))) + for (I, v) in nonzero_pairs(W.tensors) + W′.tensors[I] = $f(v) end - for (I, v) in nonzero_pairs(W.C) - W′.C[I] = $f(v) - end - for (I, v) in nonzero_pairs(W.D) - W′.D[I] = $f(v) + empty!(W′.scalars) + for (K, c) in W.scalars + W′.scalars[K] = $f(c) end return W′ end @@ -227,95 +284,68 @@ end # Indexing # -------- -@inline Base.getindex(W::JordanMPOTensor, I::CartesianIndex{4}) = W[I.I...] -@propagate_inbounds function Base.getindex(W::JordanMPOTensor, I::Vararg{Int, 4}) +@inline Base.getindex(W::JordanMPOTensor, I::Vararg{Int, 4}) = W[CartesianIndex(I)] +@propagate_inbounds function Base.getindex(W::JordanMPOTensor, I::CartesianIndex{4}) @assert I[2] == I[3] == 1 - i = I[1] - j = I[4] - if (size(W, 4) > 1 && i == 1 && j == 1) || - (size(W, 1) > 1 && i == size(W, 1) && j == size(W, 4)) - return BraidingTensor{scalartype(W), spacetype(W), storagetype(W)}(eachspace(W)[1]) - elseif i == 1 && j == size(W, 4) - return insertrightunit(insertleftunit(only(W.D), 1), 3) - elseif i == 1 - return insertleftunit(W.C[1, 1, j - 1], 1) - elseif j == size(W, 4) - return insertrightunit(W.B[i - 1, 1, 1], 3) - elseif 1 < i < size(W, 1) && 1 < j < size(W, 4) - return W.A[i - 1, 1, 1, j - 1] + c = get(W.scalars, I, nothing) + if c !== nothing + return _identity_tensor(W, I, c) + elseif haskey(W.tensors, I) + return W.tensors[I] else - return zeros(storagetype(W), eachspace(W)[i, 1, 1, j]) + return zeros(storagetype(W), eachspace(W)[I]) end end -@inline function Base.setindex!(W::JordanMPOTensor, v::MPOTensor, I::CartesianIndex{4}) - return setindex!(W, v, I.I...) +@inline function Base.setindex!(W::JordanMPOTensor, v::Union{MPOTensor, Number}, I::Vararg{Int, 4}) + return setindex!(W, v, CartesianIndex(I)) +end +@inline function Base.setindex!(W::JordanMPOTensor, v::Union{MPOTensor, Number}, I::Int) + return setindex!(W, v, CartesianIndices(W)[I]) end +# a scalar or a `BraidingTensor` is stored as a scalar identity; any other tensor is stored +# in `tensors` @propagate_inbounds function Base.setindex!( - W::JordanMPOTensor, v::MPOTensor, I::Vararg{Int, 4} + W::JordanMPOTensor, v::Union{MPOTensor, Number}, I::CartesianIndex{4} ) @assert I[2] == I[3] == 1 - i = I[1] - j = I[4] - if i == 1 && j == size(W, 4) - W.D[1] = removeunit(removeunit(v, 4), 1) - elseif i == 1 && 1 < j < size(W, 4) - W.C[1, 1, j - 1] = removeunit(v, 1) - elseif j == size(W, 4) && 1 < i < size(W, 1) - W.B[i - 1, 1, 1] = removeunit(v, 4) - elseif 1 < i < size(W, 1) && 1 < j < size(W, 4) - W.A[i - 1, 1, 1, j - 1] = v - elseif (size(W, 4) > 1 && i == 1 && j == 1) || - (size(W, 1) > 1 && i == size(W, 1) && j == size(W, 4)) - v isa BraidingTensor || throw(ArgumentError("Cannot set BraidingTensor")) + if v isa Number || v isa BraidingTensor + haskey(W.tensors, I) && delete!(W.tensors, I) + W.scalars[I] = v isa Number ? convert(scalartype(W), v) : one(scalartype(W)) else - throw(ArgumentError("Cannot set index ($i, 1, 1, $j)")) + delete!(W.scalars, I) + W.tensors[I] = v end return W end -@inline function Base.setindex!(W::JordanMPOTensor, v::MPOTensor, I::Int) - return setindex!(W, v, CartesianIndices(W)[I]) -end # Sparse functionality # -------------------- function BlockTensorKit.nonzero_keys(W::JordanMPOTensor) - nrows = size(W, 1) - ncols = size(W, 4) - p = CartesianIndex{4}[] - ncols > 1 && push!(p, CartesianIndex(1, 1, 1, 1)) - nrows > 1 && push!(p, CartesianIndex(nrows, 1, 1, ncols)) - - Ia = CartesianIndex(1, 0, 0, 1) - for I in nonzero_keys(W.A) - push!(p, I + Ia) - end - - Ib = CartesianIndex(1, 0, 0) - for I in nonzero_keys(W.B) - push!(p, CartesianIndex((I + Ib).I..., ncols)) - end - - Ic = CartesianIndex(0, 0, 1) - for I in nonzero_keys(W.C) - push!(p, CartesianIndex(1, (I + Ic).I...)) + p = collect(CartesianIndex{4}, nonzero_keys(W.tensors)) + for K in keys(W.scalars) + push!(p, K) end - - for I in nonzero_keys(W.D) - push!(p, CartesianIndex(1, 1, 1, ncols)) - end - return p end function BlockTensorKit.nonzero_values(W::JordanMPOTensor) - return Iterators.map(I -> W[I], nonzero_keys(W)) + return Iterators.flatten( + ( + nonzero_values(W.tensors), + (_identity_tensor(W, K, c) for (K, c) in W.scalars), + ) + ) end function BlockTensorKit.nonzero_pairs(W::JordanMPOTensor) - return Iterators.map(I -> I => W[I], nonzero_keys(W)) + return Iterators.flatten( + ( + nonzero_pairs(W.tensors), + (K => _identity_tensor(W, K, c) for (K, c) in W.scalars), + ) + ) end function BlockTensorKit.nonzero_length(W::JordanMPOTensor) - return nonzero_length(W.A) + nonzero_length(W.B) + nonzero_length(W.C) + - nonzero_length(W.D) + Int(size(W, 1) > 1) + Int(size(W, 4) > 1) + return nonzero_length(W.tensors) + length(W.scalars) end # linalg @@ -346,11 +376,15 @@ end function _conj_mpo(W::JordanMPOTensor) V = left_virtualspace(W)' ⊗ physicalspace(W) ← physicalspace(W) ⊗ right_virtualspace(W)' - A = _conj_mpo(W.A) - @plansor B[-1 -2; -3] ≔ conj(W.B[-1 -3; -2]) - @plansor C[-1; -2 -3] ≔ conj(W.C[-2; -1 -3]) - D = copy(adjoint(W.D)) - return JordanMPOTensor(V, A, B, C, D) + O = jordanmpotensortype(spacetype(W), storagetype(W))(undef, V) + for (I, v) in nonzero_pairs(W.tensors) + O.tensors[I] = _conj_mpo(v) + end + empty!(O.scalars) + for (K, c) in W.scalars + O.scalars[K] = conj(c) + end + return O end function add_physical_charge(O::JordanMPOTensor, charge::Sector) @@ -360,32 +394,31 @@ function add_physical_charge(O::JordanMPOTensor, charge::Sector) fuse(physicalspace(O), auxspace) ← fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O) Odst = JordanMPOTensor{scalartype(O)}(undef, Vdst) - for (I, v) in nonzero_pairs(O) - Odst[I] = add_physical_charge(v, charge) + for (I, v) in nonzero_pairs(O.tensors) + Odst.tensors[I] = add_physical_charge(v, charge) end + merge!(Odst.scalars, O.scalars) return Odst end # Utility # ------- function Base.copy(W::JordanMPOTensor) - return JordanMPOTensor(W.V, copy(W.A), copy(W.B), copy(W.C), copy(W.D)) + return typeof(W)(copy(W.tensors), copy(W.scalars)) end function Base.copy!(Wdst::JordanMPOTensor, Wsrc::JordanMPOTensor) space(Wdst) == space(Wsrc) || throw(SpaceMismatch()) - copy!(Wdst.A, Wsrc.A) - copy!(Wdst.B, Wsrc.B) - copy!(Wdst.C, Wsrc.C) - copy!(Wdst.D, Wsrc.D) + copy!(Wdst.tensors, Wsrc.tensors) + empty!(Wdst.scalars) + merge!(Wdst.scalars, Wsrc.scalars) return Wdst end # Avoid falling back to `norm(W1 - W2)` which has to convert to SparseBlockTensorMap function Base.isapprox(W1::JordanMPOTensor, W2::JordanMPOTensor; kwargs...) - return isapprox(W1.A, W2.A; kwargs...) && - isapprox(W1.B, W2.B; kwargs...) && - isapprox(W1.C, W2.C; kwargs...) && - isapprox(W1.D, W2.D; kwargs...) + isapprox(W1.tensors, W2.tensors; kwargs...) || return false + keys(W1.scalars) == keys(W2.scalars) || return false + return all(k -> isapprox(W1.scalars[k], W2.scalars[k]; kwargs...), keys(W1.scalars)) end function Base.showarg(io::IO, W::JordanMPOTensor, toplevel::Bool) @@ -394,6 +427,6 @@ function Base.showarg(io::IO, W::JordanMPOTensor, toplevel::Bool) return nothing end -function TensorKit.type_repr(::Type{<:JordanMPOTensor{E, S}}) where {E, S} - return "JordanMPOTensor{$E, " * TensorKit.type_repr(S) * ", …}" +function TensorKit.type_repr(::Type{JordanMPOTensor{E, S, A}}) where {E, S, A} + return "JordanMPOTensor{$E, " * TensorKit.type_repr(S) * ", $A}" end diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 401064493..fa957acc2 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -633,7 +633,7 @@ Base.copy(H::MPOHamiltonian) = MPOHamiltonian(map(copy, parent(H))) function Base.getproperty(H::MPOHamiltonian, sym::Symbol) if sym === :A - return map(h -> h[2:(end - 1), 1, 1, 2:(end - 1)], parent(H)) + return map(h -> h.A, parent(H)) elseif sym === :B return map(h -> h[2:(end - 1), 1, 1, end], parent(H)) elseif sym === :C @@ -649,9 +649,11 @@ function isidentitylevel(H::InfiniteMPOHamiltonian{<:JordanMPOTensor}, i::Int) if i == 1 || i == size(H[1], 1) return true else - return all(H.A) do A - return haskey(A, CartesianIndex(i - 1, 1, 1, i - 1)) && - A[i - 1, 1, 1, i - 1] isa BraidingTensor + # a diagonal level is an identity level iff every site stores a unit identity + # scalar there; pure identities live in `scalars` (genuine/scaled operators do not) + return all(parent(H)) do W + c = get(W.scalars, CartesianIndex(i, 1, 1, i), nothing) + return c !== nothing && isone(c) end end end @@ -782,12 +784,18 @@ Base.:-(H::MPOHamiltonian, λs::AbstractVector{<:Number}) = H + (-λs) Base.:-(λs::AbstractVector{<:Number}, H::MPOHamiltonian) = λs + (-H) Base.:-(H1::MPOHamiltonian, H2::MPOHamiltonian) = H1 + (-H2) +# scaling a Jordan MPO Hamiltonian scales every path exactly once, by scaling the +# transitions out of the starting level (the top row, excluding the identity corner) function VectorInterface.scale!( H::MPOHamiltonian{O}, λ::Number ) where {O <: JordanMPOTensor} - for i in 1:length(H) - scale!(H[i].C, λ) - scale!(H[i].D, λ) + for W in parent(H) + for (I, v) in nonzero_pairs(W.tensors) + I[1] == 1 && scale!(v, λ) + end + for K in collect(keys(W.scalars)) + (K[1] == 1 && K[4] != 1) && (W.scalars[K] *= λ) + end end return H end @@ -795,12 +803,15 @@ function VectorInterface.scale!( Hdst::MPOHamiltonian{<:JordanMPOTensor}, Hsrc::MPOHamiltonian{<:JordanMPOTensor}, λ::Number ) - N = check_length(Hdst, Hsrc) - for i in 1:N - scale!(Hdst[i].C, Hsrc[i].C, λ) - scale!(Hdst[i].D, Hsrc[i].D, λ) - copy!(Hdst[i].A, Hsrc[i].A) - copy!(Hdst[i].B, Hsrc[i].B) + check_length(Hdst, Hsrc) + for (Wd, Ws) in zip(parent(Hdst), parent(Hsrc)) + for (I, v) in nonzero_pairs(Ws.tensors) + Wd.tensors[I] = I[1] == 1 ? scale(v, λ) : copy(v) + end + empty!(Wd.scalars) + for (K, c) in Ws.scalars + Wd.scalars[K] = (K[1] == 1 && K[4] != 1) ? c * λ : c + end end return Hdst end diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index a91442213..ac0a01ffb 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -10,14 +10,18 @@ function left_canonicalize!( S = spacetype(W) d = sqrt(dim(P)) + # 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 + # orthogonalize second column against first Wi = W[1, 1, 1, 1] WI = removeunit(Wi, 1) - @plansor t[l; r] := conj(WI[p; p' l]) * W.C[p; p' r] + @plansor t[l; r] := conj(WI[p; p' l]) * WC[p; p' r] # TODO: the following is currently broken due to a TensorKit bug - # @plansor C′[p; p' r] := W.C[p; p' r] - WI[p; p' l] * t[l; r] + # @plansor C′[p; p' r] := WC[p; p' r] - WI[p; p' l] * t[l; r] @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r] - add!(C′, W.C) + add!(C′, WC) # QR of second column if size(W, 1) == 1 @@ -26,24 +30,24 @@ function left_canonicalize!( if dim(R) == 0 # fully truncated V = _rightunit ⊞ _rightunit - Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) - Q2 = typeof(W.C)(undef, P ← P ⊗ SumSpace{S}()) + Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) + Q2 = typeof(WC)(undef, P ← P ⊗ SumSpace{S}()) else V = ⊞(_rightunit, space(R, 1), _rightunit) scale!(Q, d) scale!(R, inv(d)) - Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ P ← P ⊗ space(R, 1)) + Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ space(R, 1)) Q2 = transpose(Q, ((2,), (1, 3))) end - H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D) + H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, WB, Q2, WD) else - tmp = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,))) + tmp = transpose(cat(insertleftunit(C′, 1), WA; dims = 1), ((3, 1, 2), (4,))) Q, R = left_orth!(tmp; alg) if dim(R) == 0 # fully truncated V = _rightunit ⊞ _rightunit - Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) - Q2 = typeof(W.C)(undef, P ← P ⊗ SumSpace{S}()) + Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) + Q2 = typeof(WC)(undef, P ← P ⊗ SumSpace{S}()) else scale!(Q, d) scale!(R, inv(d)) @@ -52,33 +56,34 @@ function left_canonicalize!( Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) V = ⊞(_rightunit, right_virtualspace(Q′), _rightunit) end - H[i] = JordanMPOTensor(codomain(W) ← P ⊗ V, Q1, W.B, Q2, W.D) + H[i] = JordanMPOTensor(codomain(W) ← P ⊗ V, Q1, WB, Q2, WD) end # absorb into next site W′ = H[i + 1] + W′A, W′B, W′C, W′D = W′.A, W′.B, W′.C, W′.D if size(W′, 4) > 1 && dim(R) != 0 - @plansor A′[l p; p' r] := R[l; r'] * W′.A[r' p; p' r] + @plansor A′[l p; p' r] := R[l; r'] * W′A[r' p; p' r] else - A′ = similar(W′.A, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′.A)) + A′ = similar(W′A, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′A)) end if size(W′, 4) > 1 - @plansor C′[l p; p' r] := t[l; r'] * W′.A[r' p; p' r] - C′ = add!(removeunit(C′, 1), W′.C) + @plansor C′[l p; p' r] := t[l; r'] * W′A[r' p; p' r] + C′ = add!(removeunit(C′, 1), W′C) else - C′ = W′.C # empty + C′ = W′C # empty end if dim(R) != 0 - @plansor B′[l p; p'] := R[l; r] * W′.B[r p; p'] + @plansor B′[l p; p'] := R[l; r] * W′B[r p; p'] else - B′ = similar(W′.B, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′.B)) + B′ = similar(W′B, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′B)) end - @plansor D′[l p; p'] := t[l; r] * W′.B[r p; p'] - D′ = add!(removeunit(D′, 1), W′.D) + @plansor D′[l p; p'] := t[l; r] * W′B[r p; p'] + D′ = add!(removeunit(D′, 1), W′D) H[i + 1] = JordanMPOTensor( right_virtualspace(H[i]) ⊗ physicalspace(W′) ← domain(W′), @@ -99,14 +104,18 @@ function right_canonicalize!( S = spacetype(W) d = sqrt(dim(P)) + # 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 + # orthogonalize second row against last Wi = W[end, 1, 1, end] WI = removeunit(Wi, 4) - @plansor t[l; r] := conj(WI[r p; p']) * W.B[l p; p'] + @plansor t[l; r] := conj(WI[r p; p']) * WB[l p; p'] # TODO: the following is currently broken due to a TensorKit bug - # @plansor B′[l p; p'] := W.B[l p; p'] - WI[r p; p'] * t[l; r] + # @plansor B′[l p; p'] := WB[l p; p'] - WI[r p; p'] * t[l; r] @plansor B′[l p; p'] := -WI[r p; p'] * t[l; r] - add!(B′, W.B) + add!(B′, WB) # LQ of second row if size(W, 4) == 1 @@ -115,24 +124,24 @@ function right_canonicalize!( if dim(R) == 0 V = _rightunit ⊞ _rightunit - Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) - Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ P ← P) + Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) + Q2 = typeof(WB)(undef, SumSpace{S}() ⊗ P ← P) else V = ⊞(_rightunit, space(Q, 1), _rightunit) scale!(Q, d) scale!(R, inv(d)) - Q1 = typeof(W.A)(undef, space(Q, 1) ⊗ P ← P ⊗ SumSpace{S}()) + Q1 = typeof(WA)(undef, space(Q, 1) ⊗ P ← P ⊗ SumSpace{S}()) Q2 = transpose(Q, ((1, 3), (2,))) end - H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, W.C, W.D) + H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, WC, WD) else B′′ = insertleftunit(B′, 4) - tmp = transpose(cat(B′′, W.A; dims = 4), ((1,), (3, 4, 2))) + tmp = transpose(cat(B′′, WA; dims = 4), ((1,), (3, 4, 2))) R, Q = right_orth!(tmp; alg) if dim(R) == 0 V = _rightunit ⊞ _rightunit - Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) - Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ P ← P) + Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) + Q2 = typeof(WB)(undef, SumSpace{S}() ⊗ P ← P) else scale!(Q, d) scale!(R, inv(d)) @@ -141,33 +150,34 @@ function right_canonicalize!( Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4) V = ⊞(_rightunit, left_virtualspace(Q′), _rightunit) end - H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, W.C, W.D) + H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, WC, WD) end # absorb into previous site W′ = H[i - 1] + W′A, W′B, W′C, W′D = W′.A, W′.B, W′.C, W′.D if size(W′, 1) > 1 && dim(R) != 0 - @plansor A′[l p; p' r] := W′.A[l p; p' r'] * R[r'; r] + @plansor A′[l p; p' r] := W′A[l p; p' r'] * R[r'; r] else - A′ = similar(W′.A, codomain(W′.A) ← physicalspace(W′.A) ⊗ left_virtualspace(H[i].A)) + A′ = similar(W′A, codomain(W′A) ← physicalspace(W′A) ⊗ left_virtualspace(H[i].A)) end if size(W′, 1) > 1 - @plansor B′[l p; p' r] := W′.A[l p; p' r'] * t[r'; r] - B′ = add!(removeunit(B′, 4), W′.B) + @plansor B′[l p; p' r] := W′A[l p; p' r'] * t[r'; r] + B′ = add!(removeunit(B′, 4), W′B) else - B′ = W′.B + B′ = W′B end if dim(R) != 0 - @plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r] + @plansor C′[p; p' r] := W′C[p; p' r'] * R[r'; r] else - C′ = similar(W′.C, codomain(W′.C) ← physicalspace(W′) ⊗ left_virtualspace(H[i].A)) + C′ = similar(W′C, codomain(W′C) ← physicalspace(W′) ⊗ left_virtualspace(H[i].A)) end - @plansor D′[p; p' r] := W′.C[p; p' r'] * t[r'; r] - D′ = add!(removeunit(D′, 3), W′.D) + @plansor D′[p; p' r] := W′C[p; p' r'] * t[r'; r] + D′ = add!(removeunit(D′, 3), W′D) H[i - 1] = JordanMPOTensor(codomain(W′) ← physicalspace(W′) ⊗ V, A′, B′, C′, D′) return H end diff --git a/test/setup/testsetup.jl b/test/setup/testsetup.jl index 68c99f096..053d272f3 100644 --- a/test/setup/testsetup.jl +++ b/test/setup/testsetup.jl @@ -63,16 +63,7 @@ function force_planar(x::SparseBlockTensorMap) return SparseBlockTensorMap{valtype(data)}(data, force_planar(space(x))) end function force_planar(W::JordanMPOTensor) - V = force_planar(space(W)) - TW = MPSKit.jordanmpotensortype(eltype(V[1]), scalartype(W)) - dst = TW(undef, V) - - for t in (:A, :B, :C, :D) - for (I, v) in nonzero_pairs(getproperty(W, t)) - getproperty(dst, t)[I] = force_planar(v) - end - end - return dst + return JordanMPOTensor(force_planar(W.tensors), copy(W.scalars)) end force_planar(mpo::MPOHamiltonian) = MPOHamiltonian(map(force_planar, parent(mpo))) force_planar(mpo::MPO) = MPO(map(force_planar, parent(mpo)))