Allow various objects to accept a storagetype or pass it along

Project.toml

@@ -45,7 +45,7 @@ OptimKit = "0.5"
 Printf = "1"
 Random = "1"
 Statistics = "1"
-TensorKit = "0.16.5, 0.17"
+TensorKit = "0.17"
 TensorKitTensors = "0.3.1"
 TensorOperations = "5"
 TupleTools = "1.6.0"

docs/src/changelog.md

@@ -32,6 +32,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 
 ### Fixed
 
+- `initializePEPS` now returns a real-valued PEPS if given a real-valued input PEPO, rather than always returning a complex-valued PEPS
+
 ### Performance
 
 ## [0.8.1](https://github.com/quantumkithub/pepskit.jl/compare/v0.8.0...v0.8.1) - 2026-07-30

src/algorithms/ctmrg/gaugefix.jl

@@ -82,7 +82,7 @@ function compute_relative_phases(
 
         # Random MPS of same bond dimension
         M = map(Tsfinal) do t
-            randn(scalartype(t), codomain(t) ← domain(t))
+            randn(storagetype(T), codomain(t) ← domain(t))
         end
 
         # Find right fixed points of mixed transfer matrices
@@ -107,7 +107,7 @@ function compute_relative_phases(envfinal::CTMRGEnv{C, T}, envprev::CTMRGEnv{C,
 
     # Random Hermitian MPS of same bond dimension
     # (make Hermitian such that T-M transfer matrix has real eigenvalues)
-    M = _project_hermitian(randn(scalartype(Tfinal), space(Tfinal)))
+    M = _project_hermitian(randn(storagetype(Tfinal), space(Tfinal)))
 
     # Find right fixed points of mixed transfer matrices
     eigsolve_alg = Lanczos() # real eigenvalues
@@ -145,7 +145,7 @@ end
 
 function initialize_right_fixedpoint(tops, bottoms)
     ρ0 = randn(
-        scalartype(tops), space(tops[end], numind(tops[end]))' ← space(bottoms[end], numind(bottoms[end]))'
+        TensorKit.promote_storagetype(tops...), space(tops[end], numind(tops[end]))' ← space(bottoms[end], numind(bottoms[end]))'
     )
     return ρ0
 end

src/algorithms/ctmrg/initialization.jl

@@ -1,11 +1,11 @@
 """
-    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization, virtual_spaces...)
+    initialize_ctmrg_environment([elt::Type,] n::InfiniteSquareNetwork, alg::RandomInitialization, virtual_spaces...)
 
 Initialize a fully random `CTMRGEnv` using the given environment virtual spaces. See
 [`CTMRGEnv`](@ref) for details on the expected format of the virtual spaces.
 """
 function initialize_ctmrg_environment(
-        elt::Type{<:Number},
+        elt::Type,
         n::InfiniteSquareNetwork,
         alg::RandomInitialization,
         virtual_spaces... = oneunit(spacetype(n)),
@@ -14,13 +14,13 @@ function initialize_ctmrg_environment(
 end
 
 """
-    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization)
+    initialize_ctmrg_environment([elt::Type,] n::InfiniteSquareNetwork, alg::RandomInitialization)
 
 Initialize a `CTMRGEnv` corresponding to a product state with trivial virtual spaces and
 corners. The product state edge tensors are initialized as `alg.f(elt, V::ProductSpace)`.
 """
 function initialize_ctmrg_environment(
-        elt::Type{<:Number},
+        elt::Type,
         n::InfiniteSquareNetwork,
         alg::ProductStateInitialization,
     )
@@ -29,14 +29,14 @@ function initialize_ctmrg_environment(
 end
 
 """
-    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization, [env0])
+    initialize_ctmrg_environment([elt::Type,] n::InfiniteSquareNetwork, alg::RandomInitialization, [env0])
 
 Initialize a `CTMRGEnv` by applying a single untruncated iteration of
 [`SimultaneousCTMRG`](@ref) to a given initial environment. By default, the starting
 environment is chosen as a random product state.
 """
 function initialize_ctmrg_environment(
-        elt::Type{<:Number},
+        elt::Type,
         n::InfiniteSquareNetwork,
         alg::ApplicationInitialization,
         env0 = ProductStateEnv(alg.f, elt, n)
@@ -50,7 +50,7 @@ _check_two_layer(::InfiniteSquareNetwork) = false
 _check_two_layer(::InfiniteSquareNetwork{<:PEPSSandwich}) = true
 
 """
-    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization, [env0])
+    initialize_ctmrg_environment([elt::Type,] n::InfiniteSquareNetwork, alg::RandomInitialization, [env0])
 
 Initialize a `CTMRGEnv` corresponding to a product state acting as an identity between the
 virtual spaces of a two-layer network, for example
@@ -65,7 +65,7 @@ virtual spaces of a two-layer network, for example
 ```
 """
 function initialize_ctmrg_environment(
-        elt::Type{<:Number},
+        elt::Type,
         n::InfiniteSquareNetwork,
         ::IdentityInitialization,
     )
@@ -86,5 +86,5 @@ function initialize_ctmrg_environment(
         elt::Type{<:Number}, A::Union{InfinitePEPS, InfinitePartitionFunction}, args...;
         kwargs...
     )
-    return initialize_ctmrg_environment(elt, InfiniteSquareNetwork(A), args...; kwargs...)
+    return initialize_ctmrg_environment(similarstoragetype(storagetype(A), elt), InfiniteSquareNetwork(A), args...; kwargs...)
 end

src/algorithms/time_evolution/gaugefix_su.jl

@@ -17,7 +17,7 @@ $(TYPEDFIELDS)
     maxiter::Int = 100
 end
 
-function _trivial_gates(elt::Type{<:Number}, lattice::Matrix{S}) where {S <: ElementarySpace}
+function _trivial_gates(elt::Type, lattice::Matrix{S}) where {S <: ElementarySpace}
     Nr, Nc = size(lattice)
     gates = map(Iterators.product(1:2, 1:Nc, 1:Nr)) do (d, c, r)
         site1 = CartesianIndex(r, c)
@@ -38,7 +38,7 @@ Fix the gauge of `psi` using trivial simple update.
 """
 function gauge_fix(psi::InfiniteState, alg::SUGauge)
     time0 = time()
-    gates = _trivial_gates(scalartype(psi), physicalspace(psi))
+    gates = _trivial_gates(storagetype(psi), physicalspace(psi))
     trunc = _get_fixedspacetrunc(psi)
     su_alg = SimpleUpdate(; trunc, bipartite = _is_bipartite(psi))
     wts0 = SUWeight(psi)

src/algorithms/time_evolution/simpleupdate3site.jl

@@ -1,14 +1,14 @@
 function _fuse_physicalspaces(O::GenericMPSTensor{S, 5}) where {S <: ElementarySpace}
     V1, V2 = codomain(O, 2), codomain(O, 3)
-    F = isomorphism(Int, fuse(V1, V2), V1 ⊗ V2)
+    F = isomorphism(similarstoragetype(O, Int), fuse(V1, V2), V1 ⊗ V2)
     @plansor O_fused[-1 -2 -4 -5; -6] := F[-2; 2 3] * O[-1 2 3 -4 -5; -6]
     return O_fused, F
 end
 
 function _unfuse_physicalspace(
         O::GenericMPSTensor{S, 4}, Vout::ElementarySpace, Vin::ElementarySpace = Vout'
     ) where {S <: ElementarySpace}
-    F = isomorphism(Int, Vout ⊗ Vin, fuse(Vout ⊗ Vin))
+    F = isomorphism(similarstoragetype(O, Int), Vout ⊗ Vin, fuse(Vout ⊗ Vin))
     @plansor O_unfused[-1 -2 -3 -4 -5; -6] := F[-2 -3; 1] * O[-1 1 -4 -5; -6]
     return O_unfused, F
 end

src/algorithms/time_evolution/time_evolve.jl

@@ -69,7 +69,7 @@ function _timeevol_sanity_check(
 end
 
 function MPSKit.infinite_temperature_density_matrix(H::LocalOperator)
-    T = scalartype(H)
+    T = storagetype(H)
     A = map(physicalspace(H)) do Vp
         ψ = permute(TensorKit.id(T, Vp), (1, 2))
         Vv = oneunit(Vp) # trivial (1D) virtual space

src/algorithms/time_evolution/trotter_gate.jl

@@ -157,7 +157,9 @@ function _trotterize_nnn2site!(gates::Vector, H::LocalOperator, dt::Number)
             term = permute(term, ((2, 1), (4, 3)))
         end
         gate = gate_to_mpo(exp(term * -dt / 2))
-        b = TensorKit.BraidingTensor{T}(physicalspace(H, x2), left_virtualspace(gate[2]))
+        A = similarstoragetype(term, T)
+        S = spacetype(H)
+        b = TensorKit.BraidingTensor{T, S, A}(physicalspace(H, x2), left_virtualspace(gate[2]))
         insert!(gate, 2, TensorMap(b))
         push!(gates, [x1, x2, x3] => gate)
     end

src/algorithms/toolbox.jl

@@ -13,7 +13,7 @@ function edge_transfer_spectrum(
         sector = one(sectortype(E))
     ) where {E <: CTMRGEdgeTensor}
     init = randn(
-        scalartype(E),
+        storagetype(E),
         space(first(bot), numind(first(bot)))' ← ℂ[typeof(sector)](sector => 1)' ⊗ space(first(top), 1),
     )
 

src/environments/bp_environments.jl

@@ -108,7 +108,7 @@ function BPEnv(
         f, T, D_north::P, D_east::P;
         unitcell::Tuple{Int, Int} = (1, 1), posdef::Bool = true
     ) where {P <: ProductSpace}
-    return BPEnv(f, T, N, _fill_edge_physical_spaces(D_north, D_east; unitcell)...; posdef)
+    return BPEnv(f, T, _fill_edge_physical_spaces(D_north, D_east; unitcell)...; posdef)
 end
 
 """
@@ -119,13 +119,13 @@ Construct a BP environment by specifying a corresponding [`InfiniteSquareNetwork
 function BPEnv(f, T, network::InfiniteSquareNetwork; posdef::Bool = true)
     Ds_north = _north_edge_physical_spaces(network)
     Ds_east = _east_edge_physical_spaces(network)
-    return BPEnv(f, T, Ds_north, Ds_east; posdef)
+    return BPEnv(f, similarstoragetype(storagetype(network), eltype(T)), Ds_north, Ds_east; posdef)
 end
 function BPEnv(network::Union{InfiniteSquareNetwork, InfinitePartitionFunction, InfinitePEPS, InfinitePEPO}, args...; kwargs...)
-    return BPEnv(isomorphism, scalartype(network), network, args...; kwargs...)
+    return BPEnv(isomorphism, storagetype(network), network, args...; kwargs...)
 end
 function BPEnv(f, T, state::Union{InfinitePartitionFunction, InfinitePEPS, InfinitePEPO}, args...; kwargs...)
-    return BPEnv(f, T, InfiniteSquareNetwork(state), args...; kwargs...)
+    return BPEnv(f, similarstoragetype(eltype(state), eltype(T)), InfiniteSquareNetwork(state), args...; kwargs...)
 end
 
 Base.eltype(::Type{BPEnv{T}}) where {T} = T
@@ -173,7 +173,7 @@ function CTMRGEnv(bp_env::BPEnv)
         return insertleftunit(insertleftunit(M), 1)
     end
     corners = map(CartesianIndices(edges)) do _
-        return TensorKit.id(scalartype(bp_env), oneunit(spacetype(bp_env)))
+        return TensorKit.id(storagetype(bp_env), oneunit(spacetype(bp_env)))
     end
     return CTMRGEnv(corners, edges)
 end

src/environments/product_state_environments.jl

@@ -82,13 +82,13 @@ Construct a product state environment by specifying a corresponding [`InfiniteSq
 function ProductStateEnv(f, T, network::InfiniteSquareNetwork)
     Ds_north = _north_edge_physical_spaces(network)
     Ds_east = _east_edge_physical_spaces(network)
-    return ProductStateEnv(f, T, Ds_north, Ds_east)
+    return ProductStateEnv(f, similarstoragetype(storagetype(network), eltype(T)), Ds_north, Ds_east)
 end
 function ProductStateEnv(network::Union{InfiniteSquareNetwork, InfinitePartitionFunction, InfinitePEPS})
-    return ProductStateEnv(randn, scalartype(network), network)
+    return ProductStateEnv(randn, storagetype(network), network)
 end
 function ProductStateEnv(f, T, state::Union{InfinitePartitionFunction, InfinitePEPS}, args...)
-    return ProductStateEnv(f, T, InfiniteSquareNetwork(state), args...)
+    return ProductStateEnv(f, similarstoragetype(eltype(state), eltype(T)), InfiniteSquareNetwork(state), args...)
 end
 
 Base.eltype(::Type{ProductStateEnv{T}}) where {T} = T

src/environments/suweight.jl

@@ -58,21 +58,43 @@ function SUWeight(wts_mats::AbstractMatrix{E}...) where {E <: PEPSWeight}
 end
 
 """
-    SUWeight(Nspaces::M, [Espaces::M]) where {M<:AbstractMatrix{<:ElementarySpace}}
+    SUWeight([T = ComplexF64,] Nspaces::M, [Espaces::M]) where {M<:AbstractMatrix{<:ElementarySpace}}
 
 Create a trivial `SUWeight` by specifying the vertical (north) or horizontal (east) virtual bond spaces.
+`T` is either a scalar type (like `ComplexF32`) or an array type (like `Vector{Float64}`) that
+specifies the storage type. If `T <: Number` the default storage type is `Vector{T}`.
 """
 function SUWeight(
+        T,
         Nspaces::M, Espaces::M = Nspaces
     ) where {M <: AbstractMatrix{<:ElementarySpace}}
     @assert size(Nspaces) == size(Espaces)
     Nr, Nc = size(Nspaces)
     weights = map(Iterators.product(1:2, 1:Nr, 1:Nc)) do (d, r, c)
         V = (d == 1 ? Espaces[r, c] : Nspaces[r, c])
-        DiagonalTensorMap(ones(reduceddim(V)), V)
+        realT = similarstoragetype(similarstoragetype(T), real(eltype(T)))
+        diag = realT(undef, reduceddim(V))
+        fill!(diag, one(real(eltype(T))))
+        DiagonalTensorMap(diag, V)
     end
     return SUWeight(weights)
 end
+SUWeight(Nspaces::M, Espaces::M = Nspaces) where {M <: AbstractMatrix{<:ElementarySpace}} =
+    SUWeight(ComplexF64, Nspaces, Espaces)
+
+"""
+    SUWeight(T, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1, 1)) where {S<:ElementarySpace}
+
+Create a trivial `SUWeight` by specifying its vertical (north) and horizontal (east)
+as `ElementarySpace`s) and unit cell size. `T` is either a scalar type (like `ComplexF32`)
+or an array type (like `Vector{Float64}`) that specifies the storage type.
+If `T <: Number` the default storage type is `Vector{T}`.
+"""
+function SUWeight(
+        T, Nspace::S, Espace::S = Nspace; unitcell::Tuple{Int, Int} = (1, 1)
+    ) where {S <: ElementarySpace}
+    return SUWeight(T, fill(Nspace, unitcell), fill(Espace, unitcell))
+end
 
 """
     SUWeight(Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1, 1)) where {S<:ElementarySpace}
@@ -83,7 +105,7 @@ as `ElementarySpace`s) and unit cell size.
 function SUWeight(
         Nspace::S, Espace::S = Nspace; unitcell::Tuple{Int, Int} = (1, 1)
     ) where {S <: ElementarySpace}
-    return SUWeight(fill(Nspace, unitcell), fill(Espace, unitcell))
+    return SUWeight(ComplexF64, fill(Nspace, unitcell), fill(Espace, unitcell))
 end
 
 """
@@ -94,7 +116,7 @@ Create a trivial `SUWeight` for a given InfinitePEPS.
 function SUWeight(peps::InfinitePEPS)
     Nspaces = map(Base.Fix2(domain, NORTH), unitcell(peps))
     Espaces = map(Base.Fix2(domain, EAST), unitcell(peps))
-    return SUWeight(Nspaces, Espaces)
+    return SUWeight(storagetype(peps), Nspaces, Espaces)
 end
 
 """
@@ -106,7 +128,7 @@ function SUWeight(pepo::InfinitePEPO)
     @assert size(pepo, 3) == 1
     Nspaces = map(Base.Fix2(domain, NORTH), @view(unitcell(pepo)[:, :, 1]))
     Espaces = map(Base.Fix2(domain, EAST), @view(unitcell(pepo)[:, :, 1]))
-    return SUWeight(Nspaces, Espaces)
+    return SUWeight(storagetype(pepo), Nspaces, Espaces)
 end
 
 Random.rand!(wts::SUWeight) = rand!(Random.default_rng(), wts)
@@ -329,7 +351,7 @@ which has the same real scalartype as ``wts`.
 """
 function CTMRGEnv(wts::SUWeight)
     _, Nr, Nc = size(wts)
-    elt = scalartype(wts)
+    elt = storagetype(wts)
     V_env = oneunit(spacetype(wts))
     edges = map(Iterators.product(1:4, 1:Nr, 1:Nc)) do (d, r, c)
         wt_idx = if d == NORTH

src/environments/vumps_environments.jl

@@ -25,18 +25,20 @@ function MPSKit.allocate_GL(
         bra::InfiniteMPS, mpo::InfiniteTransferMatrix, ket::InfiniteMPS, i::Int
     )
     T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
+    TA = similarstoragetype(storagetype(mpo), T)
     V =
         left_virtualspace(bra, i) ⊗ _elementwise_dual(left_virtualspace(mpo, i)) ←
         left_virtualspace(ket, i)
-    TT = TensorMap{T}
+    TT = TensorKit.TensorMapWithStorage{T, TA}
     return TT(undef, V)
 end
 
 function MPSKit.allocate_GR(
         bra::InfiniteMPS, mpo::InfiniteTransferMatrix, ket::InfiniteMPS, i::Int
     )
     T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
+    TA = similarstoragetype(storagetype(mpo), T)
     V = right_virtualspace(ket, i) ⊗ right_virtualspace(mpo, i) ← right_virtualspace(bra, i)
-    TT = TensorMap{T}
+    TT = TensorKit.TensorMapWithStorage{T, TA}
     return TT(undef, V)
 end

src/networks/tensors.jl

@@ -24,10 +24,11 @@ const PFTensor = PartitionFunctionTensor
 
 """
     PartitionFunctionTensor(f, ::Type{T}, Pspace::S, Nspace::S,
-               [Espace::S], [Sspace::S], [Wspace::S]) where {T,S<:Union{Int,ElementarySpace}}
+               [Espace::S], [Sspace::S], [Wspace::S]) where {T, S<:Union{Int,ElementarySpace}}
 
-Construct a PartitionFunctionTensor tensor based on the north, east, west and south spaces.
-The tensor elements are generated based on `f` and the element type is specified in `T`.
+Construct a `PartitionFunctionTensor` tensor based on the north, east, west and south spaces.
+The tensor elements are generated based on `f` and the element type and storage type is
+specified in `T`.
 """
 function PartitionFunctionTensor(
         f, ::Type{T},
@@ -80,13 +81,15 @@ const PEPSTensor{S <: ElementarySpace} = AbstractTensorMap{<:Any, S, 1, 4}
 
 """
     PEPSTensor(f, ::Type{T}, Pspace::S, Nspace::S,
-               [Espace::S], [Sspace::S], [Wspace::S]) where {T,S<:Union{Int,ElementarySpace}}
+               [Espace::S], [Sspace::S], [Wspace::S]) where {T, S<:Union{Int,ElementarySpace}}
 
 Construct a PEPS tensor based on the physical, north, east, south and west spaces.
-The tensor elements are generated based on `f` and the element type is specified in `T`.
+The tensor elements are generated based on `f` and the element type and storage type
+is specified in `T`.
 """
 function PEPSTensor(
-        f, ::Type{T},
+        f,
+        ::Type{T},
         Pspace::S,
         Nspace::S, Espace::S = Nspace, Sspace::S = Nspace', Wspace::S = Espace',
     ) where {T, S <: ElementarySpace}
@@ -155,7 +158,7 @@ herm_depth(t::PEPOTensor) = permute(t', ((5, 6), (3, 2, 1, 4)))
 Fuse the physical indices of a PEPO tensor, obtaining a PEPS tensor.
 """
 function fuse_physicalspaces(O::PEPOTensor)
-    F = isomorphism(Int, fuse(codomain(O)), codomain(O))
+    F = isomorphism(TensorKit.similarstoragetype(O, Int), fuse(codomain(O)), codomain(O))
     return F * O, F
 end
 

src/operators/infinitepepo.jl

@@ -121,7 +121,7 @@ function initializePEPS(
     end
     Nspaces = repeat([vspace], size(T, 1), size(T, 2))
     Espaces = repeat([vspace], size(T, 1), size(T, 2))
-    return InfinitePEPS(Pspaces, Nspaces, Espaces)
+    return InfinitePEPS(randn, storagetype(T), Pspaces, Nspaces, Espaces)
 end
 
 ## Unit cell interface