diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 35f7547c7..04b807804 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -1,59 +1,37 @@
 env:
-  SECRET_CODECOV_TOKEN: "MH6hHjQi7vG2V1Yfotv5/z5Dkx1k5SdyGYlGTFXiQr22XksJgsXaBuvFKUrjC7JwcpBsOVU8103LuMKl3m7VJ35WzHZrOssYycVbdGcb2kloc6xvUOsN2R5BrhCQ4Pii0l6ZeVRjCnZVkcmb0Rf4glGFyfibCrqniry8RLhblsuFKFsijRK4OxiWYEs1IvUulN+ER8tEsEtw4+ZqC5nbLGMSnUG/saPkDQOVIBscvikbKEnBcCXBheGPktF+Y/cy/1Xa+FiBPoZcypwTeAjKG1g0MqyHXjaYekb/7fekaj+hukGaeJSCXxY8KEb2IZCh+Y36Tp6y6qsIp/AdtEnCpQ==;U2FsdGVkX18WQxvGLspPwzC4aDe+U7TXU+itebTbgh8LUkE6GukxxReHYiDZ6IrBiVvSGTVJMquW0c8KsOI1pw=="
+  SECRET_CODECOV_TOKEN: "yXWVfuSYtwx4ksmqquaMEj0TWi42b1QZxPEwJccOBQ0fSZqNdkqGS/Z2fZ4bEPOoUGEzFPlqn75ZVe9nwf4XlbZCN0mPyYcA3DQiqAwO+9SskfwolCYBfI2RsBsVlZy7YXDHRH5KTGvelU0/fuuCt/DAk2j+Xk9HOZr4kx5RxEG1dKvBzUGB8q5phgJjvm0DUQ3w12iJMVeQVWU02P6dmrDPdiUJzoTzF0Gqpo4ZaKBfK9u58WP5Ao6vwNcffVoOiy4fiH93oXExOKoc2dqlKeAyLBHVXAy6wtKxZvNOrigyqKp0UEdmhkIe4iNNzo3bg3o0AIaS3HbkvPaU0ttg7Q==;U2FsdGVkX18oOnMGdD2pf5EYcmGkA11S5cIoWyXhsCDO9HMj9r0sNxou8l5EkreJ2FZDJs1AaIMf1X9tHs3WnQ=="
 
 steps:
-  - label: "Julia v1 -- CUDA"
+  - label: "Julia {{matrix.julia}} -- {{matrix.queue}} / {{matrix.group}}"
     plugins:
       - JuliaCI/julia#v1:
-          version: "1"
-      - JuliaCI/julia-test#v1: ~
+          version: "{{matrix.julia}}"
+      - JuliaCI/julia-test#v1:
+          test_args: "{{matrix.group}}"
       - JuliaCI/julia-coverage#v1:
           dirs:
             - src
             - ext
     agents:
-      queue: "cuda"
+      queue: "{{matrix.queue}}"
     if: build.message !~ /\[skip tests\]/
-    timeout_in_minutes: 60
+    timeout_in_minutes: 120
+    matrix:
+      setup:
+        julia:
+          - "1.10"
+          - "1.13"
+        queue:
+          - "cuda"
+          - "rocm"
+        group:
+          - "boundarymps"
+          - "bondenv"
+          - "bp"
+          - "compress"
+          - "ctmrg"
+          - "gradients"
+          - "timeevol"
+          - "toolbox"
+          - "utility"
 
-  - label: "Julia LTS -- CUDA"
-    plugins:
-      - JuliaCI/julia#v1:
-          version: "1.10" # "lts" isn't valid
-      - JuliaCI/julia-test#v1: ~
-      - JuliaCI/julia-coverage#v1:
-          dirs:
-            - src
-            - ext
-    agents:
-      queue: "cuda"
-    if: build.message !~ /\[skip tests\]/
-    timeout_in_minutes: 60
-  
-  - label: "Julia v1 -- AMDGPU"
-    plugins:
-      - JuliaCI/julia#v1:
-          version: "1"
-      - JuliaCI/julia-test#v1: ~
-      - JuliaCI/julia-coverage#v1:
-          dirs:
-            - src
-            - ext
-    agents:
-      queue: "rocm"
-    if: build.message !~ /\[skip tests\]/
-    timeout_in_minutes: 60
-
-  - label: "Julia LTS -- AMDGPU"
-    plugins:
-      - JuliaCI/julia#v1:
-          version: "1.10" # "lts" isn't valid
-      - JuliaCI/julia-test#v1: ~
-      - JuliaCI/julia-coverage#v1:
-          dirs:
-            - src
-            - ext
-    agents:
-      queue: "rocm"
-    if: build.message !~ /\[skip tests\]/
-    timeout_in_minutes: 60
diff --git a/Project.toml b/Project.toml
index 89027f73f..e50f6f023 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,6 +23,8 @@ OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
+StridedViews = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
 TensorKitTensors = "41b62e7d-e9d1-4e23-942c-79a97adf954b"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
@@ -30,8 +32,23 @@ TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
+[weakdeps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
+
+[sources]
+MPSKit = {rev = "main", url = "https://github.com/QuantumKitHub/MPSKit.jl"}
+MatrixAlgebraKit = {rev = "ksh/gesvdx-rank1-guard", url = "https://github.com/QuantumKitHub/MatrixAlgebraKit.jl"}
+TensorKit = {rev = "ksh/batched_svd", url = "https://github.com/QuantumKitHub/TensorKit.jl"}
+
+[extensions]
+PEPSKitAdaptExt = "Adapt"
+PEPSKitGPUArraysExt = "GPUArrays"
+
 [compat]
 Accessors = "0.1"
+Adapt = "4"
+GPUArrays = "11"
 ChainRulesCore = "1.0"
 Compat = "3.46, 4.2"
 DocStringExtensions = "0.9.3"
@@ -49,8 +66,8 @@ Random = "1"
 Statistics = "1"
 TensorKit = "0.16.5, 0.17"
 TensorKitTensors = "0.3.1"
-TensorOperations = "5"
+TensorOperations = "5.8.1"
 TupleTools = "1.6.0"
-VectorInterface = "0.4, 0.5, 0.6"
+VectorInterface = "0.4, 0.5, 0.6, 0.7"
 Zygote = "0.6, 0.7"
 julia = "1.10"
diff --git a/ext/PEPSKitAdaptExt.jl b/ext/PEPSKitAdaptExt.jl
new file mode 100644
index 000000000..23df36a83
--- /dev/null
+++ b/ext/PEPSKitAdaptExt.jl
@@ -0,0 +1,32 @@
+module PEPSKitAdaptExt
+
+using PEPSKit
+using Adapt
+
+function Adapt.adapt_structure(to, x::PEPSKit.LocalOperator{T, S}) where {T, S}
+    terms′ = Dict(k => adapt(to, v) for (k, v) in x.terms)
+    return PEPSKit.LocalOperator{valtype(terms′)}(x.lattice, terms′)
+end
+
+function Adapt.adapt_structure(to, x::PEPSKit.InfinitePEPS{T}) where {T}
+    A′ = map(a -> adapt(to, a), x.A)
+    return InfinitePEPS{eltype(A′)}(A′)
+end
+
+function Adapt.adapt_structure(to, x::PEPSKit.InfinitePEPO{T}) where {T}
+    A′ = map(a -> adapt(to, a), x.A)
+    return InfinitePEPO{eltype(A′)}(A′)
+end
+
+function Adapt.adapt_structure(to, x::PEPSKit.InfinitePartitionFunction{T}) where {T}
+    A′ = map(a -> adapt(to, a), x.A)
+    return InfinitePartitionFunction{eltype(A′)}(A′)
+end
+
+function Adapt.adapt_structure(to, x::PEPSKit.CTMRGEnv{C, T}) where {C, T}
+    C′ = map(c -> adapt(to, c), x.corners)
+    T′ = map(t -> adapt(to, t), x.edges)
+    return CTMRGEnv{eltype(C′), eltype(T′)}(C′, T′)
+end
+
+end
diff --git a/ext/PEPSKitGPUArraysExt.jl b/ext/PEPSKitGPUArraysExt.jl
new file mode 100644
index 000000000..af694983b
--- /dev/null
+++ b/ext/PEPSKitGPUArraysExt.jl
@@ -0,0 +1,179 @@
+module PEPSKitGPUArraysExt
+
+using GPUArrays
+using GPUArrays: AnyGPUArray, AllocCache
+using PEPSKit
+using TensorKit
+using TensorKit: MatrixAlgebraKit as MAK
+
+# Each caller (such as `su_iter`) gets a pair of caches. This makes sense to do on a per-caller basis
+# because what is being cached varies between algorithms.
+# For each caller we also store several caches, for SimultaneousCTMRG and SU,
+# one for even iterations and one for odd, and for SequentialCTMRG, 5 (one "round" plus one extra)
+# This has to be done because we can't reuse a cache from iteration `i`
+# until iteration `i+n` is completely finished and its result handed off.
+const ALLOC_CACHES = Dict{Tuple{Symbol, Int}, Vector{AllocCache}}()
+const ALLOC_CACHES_LOCK = ReentrantLock()
+
+function _caches(site::Symbol, depth::Int)
+    return Base.@lock ALLOC_CACHES_LOCK begin
+        get!(() -> [AllocCache() for _ in 1:depth], ALLOC_CACHES, (site, depth))
+    end
+end
+
+function PEPSKit._with_alloc_cache(f, ::Type{<:AnyGPUArray}, site::Symbol, iter::Int, depth::Int)
+    cache = @inbounds _caches(site, depth)[mod1(iter + 1, depth)]
+    return GPUArrays.@cached cache f()
+end
+
+# Reduce into a 0-dimensional device array instead of returning a host scalar. `sdiag_pow` only
+# feeds this into a broadcast, and a 0-dim array broadcasts as a scalar, so the value never has to
+# come back to the host. Returning a number here would force a device sync on every call, and
+# `sdiag_pow` runs once per bond per weight absorption in simple update.
+function PEPSKit._maxabs(data::AnyGPUArray)
+    T = real(eltype(data))
+    acc = similar(data, T, ())
+    fill!(acc, zero(T))
+    Base.mapreducedim!(abs, max, acc, data)
+    return acc
+end
+
+PEPSKit._uncache(x, ::Type{<:AnyGPUArray}) = deepcopy(x)
+
+function PEPSKit.free_alloc_caches!(::Type{<:AnyGPUArray}, caller::Symbol)
+    Base.@lock ALLOC_CACHES_LOCK begin
+        # collect first: freeing mutates ALLOC_CACHES
+        stale = [key for key in keys(ALLOC_CACHES) if first(key) === caller]
+        for key in stale
+            for cache in ALLOC_CACHES[key]
+                GPUArrays.unsafe_free!(cache)
+            end
+            delete!(ALLOC_CACHES, key)
+        end
+    end
+    return nothing
+end
+
+function PEPSKit.free_alloc_caches!(::Type{<:AnyGPUArray})
+    Base.@lock ALLOC_CACHES_LOCK begin
+        for caches in values(ALLOC_CACHES), cache in caches
+            GPUArrays.unsafe_free!(cache)
+        end
+        empty!(ALLOC_CACHES)
+    end
+    return nothing
+end
+
+
+# Batched truncated SVD of a whole cluster's internal bonds to avoid multiple small kernel launches.
+function PEPSKit.bond_svds(
+        ::Type{<:AnyGPUArray}, rls::AbstractVector, truncs::AbstractVector
+    )
+    isempty(rls) && return map(_ -> nothing, rls)
+    # The different GPU libaries offer different batching algos,
+    # make sure we have one that actually works.
+    alg = MAK.default_algorithm(MAK.batched_svd_compact!, eltype(rls))
+    Fs = map(rl -> MAK.initialize_output(MAK.svd_compact!, rl, alg), rls)
+    balg = _cluster_batched_alg(rls)
+    if isnothing(balg)
+        for (rl, F) in zip(rls, Fs)
+            MAK.svd_compact!(rl, F, alg)
+        end
+    else
+        # Pool every (bond, sector) block into one ragged batch. MatrixAlgebraKit batches
+        # blocks of equal size together even across different bonds, since the decomposition
+        # does not care which bond a block came from, and zero-pads the leftovers.
+        items = [(i, c) for i in eachindex(rls) for c in blocksectors(rls[i])]
+        As = [block(rls[i], c) for (i, c) in items]
+        Us = [block(Fs[i][1], c) for (i, c) in items]
+        Ss = [TensorKit.diagview(block(Fs[i][2], c)) for (i, c) in items]
+        Vᴴs = [block(Fs[i][3], c) for (i, c) in items]
+        MAK.batched_svd_compact!(As, (Us, Ss, Vᴴs), balg)
+    end
+    return map(Fs, truncs) do F, trunc
+        (U, S, Vᴴ) = F
+        USVᴴtrunc, ind = MAK.truncate(MAK.svd_trunc!, (U, S, Vᴴ), trunc)
+        ϵ = MAK.truncation_error!(TensorKit.diagview(S), ind)
+        return (USVᴴtrunc..., ϵ)
+    end
+end
+
+"""
+    CLUSTER_BATCHED_SVD[]
+
+Whether simple update batches the SVDs of a cluster's internal bonds into one call.
+Off by default.
+"""
+const CLUSTER_BATCHED_SVD = Ref(false)
+
+# Which batched algorithm the backend offers for the cluster's blocks, or `nothing`.
+function _cluster_batched_alg(rls::AbstractVector)
+    CLUSTER_BATCHED_SVD[] || return nothing
+    for i in eachindex(rls), c in blocksectors(rls[i])
+        return _batched_spectra_alg(block(rls[i], c))
+    end
+    return nothing
+end
+
+"""
+    _batched_spectra_alg(proto) -> alg or nothing
+
+Default batched SVD algorithm this backend offers, or `nothing` if it has none.
+"""
+function _batched_spectra_alg(proto)
+    # TODO BAD MAKE THIS A MAK CALL
+    alg = try
+        MAK.default_svd_algorithm(typeof(similar(proto, 0, 0, 0)))
+    catch
+        return nothing
+    end
+    return alg isa MAK.AbstractAlgorithm ? alg : nothing
+end
+
+# Hook into the collection-level convergence API. Deliberately restricted to the generic
+# CTMRG algorithms: `C4vCTMRG` overrides `corner_spectrum` to `eigh_vals` (its corners are
+# diagonal), so a blanket override here would silently switch it back to `svd_vals`.
+function PEPSKit.corner_spectra(
+        Cs::AbstractArray{<:AbstractTensorMap},
+        ::Union{PEPSKit.SequentialCTMRG, PEPSKit.SimultaneousCTMRG},
+    )
+    return _batched_spectra(Cs)
+end
+function PEPSKit.edge_spectra(
+        Ts::AbstractArray{<:AbstractTensorMap},
+        ::Union{PEPSKit.SequentialCTMRG, PEPSKit.SimultaneousCTMRG},
+    )
+    return _batched_spectra(Ts)
+end
+
+# `calc_convergence` decomposes every corner and every edge of the environment
+# for regular CTMRG, which is expensive, at *least* 8 separate `svd_vals` calls.
+# Across *multiple tensors* the situation is much better than within *one*,
+# because the corners have to all share a space,
+# so for a given sector their blocks have identical sizes and batch with no padding at all.
+function _batched_spectra(ts::AbstractArray{T}) where {T <: AbstractTensorMap}
+    # TODO BAD FIND A BETTER DISPATCH HERE
+    (isempty(ts) || !(TensorKit.storagetype(T) <: AnyGPUArray)) && return map(svd_vals, ts)
+    items = [(i, c) for i in eachindex(ts) for c in blocksectors(ts[i])]
+    isempty(items) && return map(svd_vals, ts)
+    alg = _batched_spectra_alg(block(ts[first(items)[1]], first(items)[2]))
+    isnothing(alg) && return map(svd_vals, ts)
+
+    Ss = map(
+        t -> MAK.initialize_output(
+            MAK.svd_vals!, t,
+            MAK.default_algorithm(
+                MAK.svd_vals!, typeof(t)
+            )
+        ), ts
+    )
+    # `batched_svd_vals!` destroys the blocks it has to decompose one at a time, but `ts` is the
+    # live environment, and computing the convergence spectra must not damage the environment
+    # it is measuring. Copying whole tensors costs one copy per tensor instead of one per block.
+    ts′ = map(copy, ts)
+    As = [block(ts′[i], c) for (i, c) in items]
+    MAK.batched_svd_vals!(As, [block(Ss[i], c) for (i, c) in items], alg)
+    return Ss
+end
+
+end
diff --git a/src/Defaults.jl b/src/Defaults.jl
index 17c22c72d..088ca1c54 100644
--- a/src/Defaults.jl
+++ b/src/Defaults.jl
@@ -52,7 +52,7 @@ Module containing default algorithm parameter values and arguments.
 
 ## `eigh` forward & reverse
 
-* `eigh_fwd_alg=:$(Defaults.eigh_fwd_alg)` : `eigh` algorithm that is used in the forward pass.
+* `eigh_fwd_alg=:$(Defaults.eigh_fwd_alg)` : `eigh` algorithm that is used in the forward pass. **Note** that on GPU, `DivideAndConquer` is much more performant than `QRIteration`.
     - `:DefaultAlgorithm` : MatrixAlgebraKit's default Eigh algorithm for a given matrix type.
     - `:DivideAndConquer` : MatrixAlgebraKit's [`DivideAndConquer`](@extref MatrixAlgebraKit.DivideAndConquer)
     - `:QRIteration` : MatrixAlgebraKit's [`QRIteration`](@extref MatrixAlgebraKit.QRIteration)
diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 827a16689..f3a5c08b2 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -18,11 +18,14 @@ using TensorKit
 using TensorKit: AdjointTensorMap, SectorDict
 using TensorKit: throw_invalid_innerproduct, similarstoragetype
 using TensorKit.Factorizations: TruncationSpace, _notrunc_ind
+import TensorKit: storagetype
 
 using KrylovKit
 using KrylovKit: Lanczos, BlockLanczos
 
-using TensorOperations, OptimKit
+using TensorOperations
+using TensorOperations: AbstractBackend, DefaultBackend, DefaultAllocator
+using OptimKit
 using ChainRulesCore, Zygote
 using LoggingExtras
 import TupleTools
@@ -54,6 +57,7 @@ include("Defaults.jl")  # Include first to allow for docstring interpolation wit
 include("utility/util.jl")
 include("utility/contraction_labels.jl")
 include("utility/tensor_traces.jl")
+include("utility/alloc_cache.jl")
 include("utility/indexing.jl")
 include("utility/diffable_threads.jl")
 include("utility/twistdual.jl")
diff --git a/src/algorithms/bp/beliefpropagation.jl b/src/algorithms/bp/beliefpropagation.jl
index 5f64ed0b7..57c7e97b3 100644
--- a/src/algorithms/bp/beliefpropagation.jl
+++ b/src/algorithms/bp/beliefpropagation.jl
@@ -39,6 +39,7 @@ function leading_boundary(env₀::BPEnv, network::InfiniteSquareNetwork, alg::Be
         ϵ = Inf
         @infov 1 loginit!(log, ϵ)
         for iter in 1:(alg.maxiter)
+            # TODO investigate why caching doesn't help here and actually makes things worse
             env′ = bp_iteration(network, env, alg)
             ϵ = oftype(ϵ, tr_distance(env, env′))
             env = env′
diff --git a/src/algorithms/bp/gaugefix.jl b/src/algorithms/bp/gaugefix.jl
index 15693557a..a86599e04 100644
--- a/src/algorithms/bp/gaugefix.jl
+++ b/src/algorithms/bp/gaugefix.jl
@@ -129,7 +129,7 @@ function SUWeight(env::BPEnv)
         I = CartesianIndex(mod1(dir′ + 1, 2), row, col)
         sqrtM12, _, sqrtM21, _ = _sqrt_bp_messages(I, env)
         Λ = DiagonalTensorMap(svd_vals!(sqrtM12 * sqrtM21))
-        return isdual(space(sqrtM12, 1)) ? _fliptwist_s(Λ) : Λ
+        return isdual(space(sqrtM12, 1)) ? _fliptwist_s!(Λ) : Λ
     end
     return SUWeight(wts)
 end
diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl
index 3e6643325..19c982b35 100644
--- a/src/algorithms/contractions/bondenv/gaugefix.jl
+++ b/src/algorithms/contractions/bondenv/gaugefix.jl
@@ -21,9 +21,9 @@ function positive_approx(benv::AbstractTensorMap{T, S, N, N}) where {T, S, N}
     # If `benv` is negative (e.g. obtained approximately from CTMRG),
     # we can multiply it by (-1).
     data = D.data
-    @inbounds for i in eachindex(data)
-        d = (sgn == -1) ? -data[i] : data[i]
-        data[i] = (d > 0) ? sqrt(d) : zero(d)
+    map!(data, data) do d
+        d2 = (sgn < 0) ? -d : d
+        return (d2 > 0) ? sqrt(d2) : zero(d2)
     end
     Z = D * U'
     return Z
diff --git a/src/algorithms/contractions/ctmrg/characteristic_equations.jl b/src/algorithms/contractions/ctmrg/characteristic_equations.jl
index c554f95e1..666bce30f 100644
--- a/src/algorithms/contractions/ctmrg/characteristic_equations.jl
+++ b/src/algorithms/contractions/ctmrg/characteristic_equations.jl
@@ -272,12 +272,10 @@ function ChainRulesCore.rrule(
     return C, squareroot_pullback
 end
 function _squareroot_pullback(C::AbstractMatrix)
-    Fdata = similar(C)
-    for j in axes(Fdata, 2), i in axes(Fdata, 1)
-        # Taking the diagonal only is okay, when dA is diagonal anyway: Fdata[i, i] = 1 / (2 * conj(C[i, i]))
-        Fdata[i, j] = 1 / conj(C[i, i] + C[j, j])
-    end
-    return Fdata
+    # Taking the diagonal only is okay, when dA is diagonal anyway: Fdata[i, i] = 1 / (2 * conj(C[i, i]))
+    Cd = diagview(C)
+    Cdt = transpose(Cd)
+    return @. 1 / conj(Cd + Cdt)
 end
 
 # take fourth root of diagonal TensorMap, but with complex non-diagonal adjoint
diff --git a/src/algorithms/contractions/vumps_contractions.jl b/src/algorithms/contractions/vumps_contractions.jl
index 44d10e498..bee3cbd0f 100644
--- a/src/algorithms/contractions/vumps_contractions.jl
+++ b/src/algorithms/contractions/vumps_contractions.jl
@@ -4,27 +4,28 @@
 
 function MPSKit.transfer_left(
         GL::GenericMPSTensor{S, N}, O::Union{PEPSSandwich, PEPOSandwich},
-        A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N},
+        A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N}; kwargs...
     ) where {S, N}
     Ā = twistdual(Ā, 2:N)
-    return mps_transfer_left(GL, O, A, Ā)
+    return mps_transfer_left(GL, O, A, Ā; kwargs...)
 end
 
 function MPSKit.transfer_right(
         GR::GenericMPSTensor{S, N}, O::Union{PEPSSandwich, PEPOSandwich},
-        A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N},
+        A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N}; kwargs...
     ) where {S, N}
     Ā = twistdual(Ā, 2:N)
-    return mps_transfer_right(GR, O, A, Ā)
+    return mps_transfer_right(GR, O, A, Ā; kwargs...)
 end
 
 ## PEPS
 
 function mps_transfer_left(
         GL::GenericMPSTensor{S, 3}, O::PEPSSandwich,
-        A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3},
+        A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3};
+        backend = DefaultBackend(), allocator = DefaultAllocator()
     ) where {S}
-    return @autoopt @tensor GL′[χ_SE D_E_above D_E_below; χ_NE] :=
+    return @autoopt @tensor backend = backend allocator = allocator GL′[χ_SE D_E_above D_E_below; χ_NE] :=
         GL[χ_SW D_W_above D_W_below; χ_NW] *
         conj(Ā[χ_SW D_S_above D_S_below; χ_SE]) *
         ket(O)[d; D_N_above D_E_above D_S_above D_W_above] *
@@ -34,9 +35,10 @@ end
 
 function mps_transfer_right(
         GR::GenericMPSTensor{S, 3}, O::PEPSSandwich,
-        A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3},
+        A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3};
+        backend = DefaultBackend(), allocator = DefaultAllocator()
     ) where {S}
-    return @autoopt @tensor GR′[χ_NW D_W_above D_W_below; χ_SW] :=
+    return @autoopt @tensor backend = backend allocator = allocator GR′[χ_NW D_W_above D_W_below; χ_SW] :=
         GR[χ_NE D_E_above D_E_below; χ_SE] *
         conj(Ā[χ_SW D_S_above D_S_below; χ_SE]) *
         ket(O)[d; D_N_above D_E_above D_S_above D_W_above] *
@@ -48,7 +50,8 @@ end
 
 @generated function mps_transfer_left(
         GL::GenericMPSTensor{S, N}, O::PEPOSandwich{H},
-        A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N},
+        A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N};
+        backend = DefaultBackend(), allocator = DefaultAllocator()
     ) where {S, N, H}
     # sanity check
     @assert H == N - 3
@@ -67,12 +70,13 @@ end
         pepo_es...,
     )
 
-    return macroexpand(@__MODULE__, :(return @autoopt @tensor $GL´_e := $rhs))
+    return macroexpand(@__MODULE__, :(return @autoopt @tensor backend = backend allocator = allocator $GL´_e := $rhs))
 end
 
 @generated function mps_transfer_right(
         GR::GenericMPSTensor{S, N}, O::PEPOSandwich{H},
-        A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N},
+        A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N};
+        backend = DefaultBackend(), allocator = DefaultAllocator()
     ) where {S, N, H}
     # sanity check
     @assert H == N - 3
@@ -91,7 +95,7 @@ end
         pepo_es...,
     )
 
-    return macroexpand(@__MODULE__, :(return @autoopt @tensor $GR´_e := $rhs))
+    return macroexpand(@__MODULE__, :(return @autoopt @tensor backend = backend allocator = allocator $GR´_e := $rhs))
 end
 
 @generated function environment_overlap(
@@ -120,46 +124,46 @@ end
 const PEPS_C_Hamiltonian{S, N} = MPSKit.MPO_C_Hamiltonian{
     <:GenericMPSTensor{S, N}, <:GenericMPSTensor{S, N},
 } # this one is technically type-piracy
-PEPS_C_Hamiltonian(GL, GR) = MPSKit.MPODerivativeOperator(GL, (), GR)
+PEPS_C_Hamiltonian(GL, GR, backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()) = MPSKit.MPODerivativeOperator(GL, (), GR, backend, allocator)
 
 const PEPS_AC_Hamiltonian{S, N} = MPSKit.MPO_AC_Hamiltonian{
     <:GenericMPSTensor{S, N}, <:PEPSSandwich, <:GenericMPSTensor{S, N},
 }
-PEPS_AC_Hamiltonian(GL, O, GR) = MPSKit.MPODerivativeOperator(GL, (O,), GR)
+PEPS_AC_Hamiltonian(GL, O, GR, backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()) = MPSKit.MPODerivativeOperator(GL, (O,), GR, backend, allocator)
 
 const PEPS_AC2_Hamiltonian{S, N} = MPSKit.MPO_AC2_Hamiltonian{
     <:GenericMPSTensor{S, N}, <:PEPSSandwich, <:PEPSSandwich, <:GenericMPSTensor{S, N},
 }
-PEPS_AC2_Hamiltonian(GL, O1, O2, GR) = MPSKit.MPODerivativeOperator(GL, (O1, O2), GR)
+PEPS_AC2_Hamiltonian(GL, O1, O2, GR, backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()) = MPSKit.MPODerivativeOperator(GL, (O1, O2), GR, backend, allocator)
 
 # Constructors
 #
-function MPSKit.C_hamiltonian(site::Int, below, ::InfiniteTransferMatrix, above, envs; kwargs...)
+function MPSKit.C_hamiltonian(site::Int, below, ::InfiniteTransferMatrix, above, envs; backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator(), kwargs...)
     GL = leftenv(envs, site + 1, below)
     GL = twistdual(GL, 1)
     GR = rightenv(envs, site, below)
     GR = twistdual(GR, numind(GR))
-    return PEPS_C_Hamiltonian(GL, GR)
+    return PEPS_C_Hamiltonian(GL, GR, backend, allocator)
 end
 
 function MPSKit.AC_hamiltonian(
-        site::Int, below, operator::InfiniteTransferPEPS, above, envs; kwargs...
+        site::Int, below, operator::InfiniteTransferPEPS, above, envs; backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator(), kwargs...
     )
     GL = leftenv(envs, site, below)
     GL = twistdual(GL, 1)
     GR = rightenv(envs, site, below)
     GR = twistdual(GR, numind(GR))
-    return PEPS_AC_Hamiltonian(GL, operator[site], GR)
+    return PEPS_AC_Hamiltonian(GL, operator[site], GR, backend, allocator)
 end
 
 function MPSKit.AC2_hamiltonian(
-        site::Int, below, operator::InfiniteTransferPEPS, above, envs; kwargs...
+        site::Int, below, operator::InfiniteTransferPEPS, above, envs; backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator(), kwargs...
     )
     GL = leftenv(envs, site, below)
     GL = twistdual(GL, 1)
     GR = rightenv(envs, site + 1, below)
     GR = twistdual(GR, numind(GR))
-    return PEPS_AC2_Hamiltonian(GL, operator[site], operator[site + 1], GR)
+    return PEPS_AC2_Hamiltonian(GL, operator[site], operator[site + 1], GR, backend, allocator)
 end
 
 # Actions
@@ -213,16 +217,16 @@ end
 const PEPO_AC_Hamiltonian{S, N, H} = MPSKit.MPO_AC_Hamiltonian{
     <:GenericMPSTensor{S, N}, <:PEPOSandwich{H}, <:GenericMPSTensor{S, N},
 }
-PEPO_AC_Hamiltonian(GL, O, GR) = MPSKit.MPODerivativeOperator(GL, (O,), GR)
+PEPO_AC_Hamiltonian(GL, O, GR, backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator()) = MPSKit.MPODerivativeOperator(GL, (O,), GR, backend, allocator)
 
 function MPSKit.AC_hamiltonian(
-        site::Int, below, operator::InfiniteTransferPEPO, above, envs; kwargs...
+        site::Int, below, operator::InfiniteTransferPEPO, above, envs; backend::AbstractBackend = DefaultBackend(), allocator = DefaultAllocator(), kwargs...
     )
     GL = leftenv(envs, site, below)
     GL = twistdual(GL, 1)
     GR = rightenv(envs, site, below)
     GR = twistdual(GR, numind(GR))
-    return PEPO_AC_Hamiltonian(GL, operator[site], GR)
+    return PEPO_AC_Hamiltonian(GL, operator[site], GR, backend, allocator)
 end
 
 @generated function (h::PEPO_AC_Hamiltonian{S, N, H})(AC::GenericMPSTensor{S, N}) where {S, N, H}
diff --git a/src/algorithms/ctmrg/c4v.jl b/src/algorithms/ctmrg/c4v.jl
index a4d0da08c..b30c00c3e 100644
--- a/src/algorithms/ctmrg/c4v.jl
+++ b/src/algorithms/ctmrg/c4v.jl
@@ -212,7 +212,7 @@ Initialize a C₄ᵥ-symmetric `CTMRGEnv` on virtual spaces `Venv` with random e
 by `f` and scalartype `T`.
 """
 function initialize_random_c4v_env(state, Venv::ElementarySpace)
-    return initialize_random_c4v_env(randn, scalartype(state), state, Venv)
+    return initialize_random_c4v_env(randn, storagetype(state), state, Venv)
 end
 function initialize_random_c4v_env(f, T, state::InfinitePEPS, Venv::ElementarySpace)
     Vpeps = north_virtualspace(state, 1, 1)'
@@ -223,28 +223,39 @@ function initialize_random_c4v_env(f, T, state::InfinitePartitionFunction, Venv:
     return initialize_random_c4v_env(f, T, Vpf, Venv)
 end
 function initialize_random_c4v_env(f, T, Vstate::VectorSpace, Venv::ElementarySpace)
-    corner₀ = DiagonalTensorMap(randn(real(T), Venv ← Venv))
+    corner₀ = DiagonalTensorMap(randn(similarstoragetype(T, real(eltype(T))), Venv ← Venv))
     edge₀ = f(T, Venv ⊗ Vstate ← Venv)
     edge₀ = _project_hermitian(edge₀)
     return CTMRGEnv(corner₀, edge₀)
 end
 
 """
-    initialize_singlet_c4v_env([T=scalartype(state)], state::InfinitePEPS, Venv::ElementarySpace)
+    initialize_singlet_c4v_env([T=storagetype(state)], state::InfinitePEPS, Venv::ElementarySpace)
 
 Initialize a C₄ᵥ-symmetric `CTMRGEnv` with a singlet corner of dimension `dim(Venv)` and an
 identity edge from `id(T, Venv ⊗ Vpeps)`.
 """
 function initialize_singlet_c4v_env(state::InfinitePEPS, Venv::ElementarySpace)
-    return initialize_singlet_c4v_env(scalartype(state), state, Venv)
+    return initialize_singlet_c4v_env(storagetype(state), state, Venv)
 end
 function initialize_singlet_c4v_env(T, state::InfinitePEPS, Venv::ElementarySpace)
     Vpeps = north_virtualspace(state, 1, 1)'
-    return initialize_singlet_c4v_env(T, Vpeps, Venv)
+    return initialize_singlet_c4v_env(similarstoragetype(storagetype(state), real(eltype(T))), Vpeps, Venv)
 end
-function initialize_singlet_c4v_env(T, Vpeps::ElementarySpace, Venv::ElementarySpace)
-    corner₀ = DiagonalTensorMap(zeros(real(T), Venv ← Venv))
+function initialize_singlet_c4v_env(T::Type{<:Number}, Vpeps::ElementarySpace, Venv::ElementarySpace)
+    realT = real(T)
+    diag = zeros(realT, dim(Venv))
+    corner₀ = DiagonalTensorMap(diag, Venv)
     corner₀.data[1] = one(real(T))
     edge₀ = permute(id(T, Venv ⊗ Vpeps), ((1, 2, 4), (3,)))
     return CTMRGEnv(corner₀, edge₀)
 end
+function initialize_singlet_c4v_env(T::Type{<:AbstractArray}, Vpeps::ElementarySpace, Venv::ElementarySpace)
+    realT = similarstoragetype(T, real(eltype(T)))
+    diag = realT(undef, dim(Venv))
+    fill!(diag, 1)
+    corner₀ = DiagonalTensorMap(diag, Venv)
+    corner₀.data[2:end] .= zero(real(eltype(T)))
+    edge₀ = permute(id(T, Venv ⊗ Vpeps), ((1, 2, 4), (3,)))
+    return CTMRGEnv(corner₀, edge₀)
+end
diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl
index f8cc6e53f..d61d0e432 100644
--- a/src/algorithms/ctmrg/ctmrg.jl
+++ b/src/algorithms/ctmrg/ctmrg.jl
@@ -49,6 +49,12 @@ Perform a single CTMRG iteration in which all directions are being grown and ren
 """
 function ctmrg_iteration(network, env, alg::CTMRGAlgorithm) end
 
+# Signature of the buffer sizes a CTMRG run allocates: the corner and edge spaces fix them,
+# and the edges carry the network bond dimension as well as the environment dimension.
+function _ctmrg_cache_signature(env::CTMRGEnv)
+    return hash((alloc_cache_signature(env.corners), alloc_cache_signature(env.edges)))
+end
+
 """
     leading_boundary(env₀, network; kwargs...) -> env, info
     # expert version:
@@ -111,6 +117,11 @@ function leading_boundary(
         env₀::CTMRGEnv, network::InfiniteSquareNetwork, alg::CTMRGAlgorithm
     )
     check_input(leading_boundary, network, env₀, alg)
+    # cached buffer sizes are set by the corner/edge spaces, and the edges carry the network
+    # bond dimension too, so a change here means every pooled buffer has gone stale
+    ignore_derivatives() do
+        free_stale_alloc_caches!(storagetype(env₀), :ctmrg, :enter, _ctmrg_cache_signature(env₀))
+    end
     log = ignore_derivatives(() -> MPSKit.IterLog("CTMRG"))
     return LoggingExtras.withlevel(; alg.verbosity) do
         env = deepcopy(env₀)
@@ -122,7 +133,10 @@ function leading_boundary(
         local info_iter
         converged = false
         for iter in 1:(alg.maxiter)
-            env, info_iter = ctmrg_iteration(network, env, alg)
+            env, info_iter = with_alloc_cache(storagetype(env), :ctmrg, iter, alloc_cache_depth(alg)) do
+                ctmrg_iteration(network, env, alg)
+            end
+            env = uncache(env, storagetype(env))
             η, CS, TS = calc_convergence(env, CS, TS, alg)
 
             if η ≤ alg.tol && iter ≥ alg.miniter
@@ -136,6 +150,15 @@ function leading_boundary(
                 ctmrg_logiter!(log, iter, η, network, env)
             end
         end
+        env = uncache(env, storagetype(env))
+        # a truncation that is not fixed-space grows the corner/edge spaces as the loop runs,
+        # leaving one pooled buffer set per intermediate shape. `env` is copied out by now, so
+        # nothing live is backed by the cache and it is safe to drop those here. Re-checking
+        # the signature keeps the pool warm for the repeated fixed-space calls of an
+        # optimization loop, where the spaces do not move.
+        ignore_derivatives() do
+            free_stale_alloc_caches!(storagetype(env), :ctmrg, :exit, _ctmrg_cache_signature(env))
+        end
         info = (;
             converged,
             convergence_error = η,
@@ -195,7 +218,6 @@ function _singular_value_distance(S₁::SV, S₂::SV) where {SV <: TensorKit.Sec
     for (c, b) in blocks(S₂)
         diff[c][1:length(b)] .-= b
     end
-
     return norm(diff)
 end
 _singular_value_distance(S₁::DiagonalTensorMap, S₂::DiagonalTensorMap) =
@@ -217,16 +239,38 @@ Spectra of the tensors that [`convergence_tensors`](@ref) selects.
 """
 function convergence_spectra(env::CTMRGEnv, alg)
     corners, edges = convergence_tensors(env, alg)
-    return map(C -> corner_spectrum(C, alg), corners), map(T -> edge_spectrum(T, alg), edges)
+    return corner_spectra(corners, alg), edge_spectra(edges, alg)
 end
 
+"""
+    corner_spectra(Cs, alg)
+    edge_spectra(Ts, alg)
+
+Spectra of a whole collection of corners or edges.
+
+Separate from [`corner_spectrum`](@ref) so that a backend can decompose the entire
+collection in one go rather than one tensor at a time. That matters on GPU: each tensor
+here carries only a handful of sector blocks, so decomposing them individually is dominated
+by per-call overhead, while a corner or edge array holds tens of tensors whose blocks share
+sizes and can be batched.
+
+The defaults just map [`corner_spectrum`](@ref) / [`edge_spectrum`](@ref) over the
+collection, so any algorithm that overrides those keeps its behaviour.
+"""
+corner_spectra(Cs, alg) = map(C -> corner_spectrum(C, alg), Cs)
+edge_spectra(Ts, alg) = map(T -> edge_spectrum(T, alg), Ts)
+
 """
     corner_spectrum(C, alg)
-    edge_spectrum(T, alg)
 
-The spectrum of a corner or edge, used to measure CTMRG convergence.
+The spectrum of a corner, used to measure CTMRG convergence.
 """
 corner_spectrum(C, alg) = svd_vals(C)
+"""
+    edge_spectrum(T, alg)
+
+The spectrum of an edge, used to measure CTMRG convergence.
+"""
 edge_spectrum(T, alg) = svd_vals(T)
 
 function calc_convergence(env, CS_old, TS_old, alg)
diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl
index 6ff44e554..9c5c73d78 100644
--- a/src/algorithms/ctmrg/gaugefix.jl
+++ b/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
diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index 235611e0b..c44f04826 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -5,7 +5,7 @@ Initialize a fully random `CTMRGEnv` using the given environment virtual spaces.
 [`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)),
@@ -20,7 +20,7 @@ Initialize a `CTMRGEnv` corresponding to a product state with trivial virtual sp
 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,
     )
@@ -36,7 +36,7 @@ Initialize a `CTMRGEnv` by applying a single untruncated iteration of
 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)
@@ -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
diff --git a/src/algorithms/optimization/implicit_differentiation.jl b/src/algorithms/optimization/implicit_differentiation.jl
index f00632504..1553df987 100644
--- a/src/algorithms/optimization/implicit_differentiation.jl
+++ b/src/algorithms/optimization/implicit_differentiation.jl
@@ -496,7 +496,7 @@ function _rrule(
     UL = left_null(U)
 
     # instantiate the differentiable variables corresponding to the intermediate projector of the contraction algorithm
-    u = zeros(scalartype(U), space(UL, numind(UL))' ← space(U, numind(U))')
+    u = zeros(storagetype(U), space(UL, numind(UL))' ← space(U, numind(U))')
 
     # prepare pullback of C4v CTMRG environment constructor (artefact of reusing asymmetric environment type for C4v symmetric contraction)
     _, c4v_env_vjp = rrule_via_ad(config, CTMRGEnv, C, E)
@@ -640,10 +640,10 @@ function PEPSKit._rrule(
 
     # instantiate the variables used in the characteristic equations
     u = map(zip(U, UL)) do (Uc, ULc)
-        return zeros(scalartype(Uc), space(ULc, numind(ULc))' ← space(Uc, numind(Uc))')
+        return zeros(storagetype(Uc), space(ULc, numind(ULc))' ← space(Uc, numind(Uc))')
     end
     v = map(zip(V, VR)) do (Vc, VRc)
-        return zeros(scalartype(Vc), space(Vc, 1) ← space(VRc, 1))
+        return zeros(storagetype(Vc), space(Vc, 1) ← space(VRc, 1))
     end
     is = sdiag_pow.(s, -1) # also treat them as general complex tensors
 
diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl
index 634292d3e..9043d05e6 100644
--- a/src/algorithms/optimization/peps_optimization.jl
+++ b/src/algorithms/optimization/peps_optimization.jl
@@ -299,7 +299,8 @@ function fixedpoint(
                 alg.reuse_env && update!(env, env′)
                 tracked_finalizer.contraction_metrics[end] = info.contraction_metrics
             end
-            return cost_function(ψ, env′, operator)
+            cf = cost_function(ψ, env′, operator)
+            return cf
         end
         g = only(gs)  # `withgradient` returns tuple of gradients `gs`
         tracked_finalizer.gradnorms_unitcell[end] = norm.(unitcell(g))
diff --git a/src/algorithms/select_algorithm.jl b/src/algorithms/select_algorithm.jl
index c121fb3fd..dfb30792e 100644
--- a/src/algorithms/select_algorithm.jl
+++ b/src/algorithms/select_algorithm.jl
@@ -155,6 +155,5 @@ function select_algorithm(
         rrule_alg = (; tol = 1.0e1tol, verbosity = verbosity - 2, krylovdim, decomposition_alg.rrule_alg...)
         decomposition_alg = (; rrule_alg, decomposition_alg...)
     end
-
     return CTMRGAlgorithm(; alg, tol, verbosity, decomposition_alg, kwargs...)
 end
diff --git a/src/algorithms/time_evolution/apply_gate.jl b/src/algorithms/time_evolution/apply_gate.jl
index 63f6e50b7..9b6f811aa 100644
--- a/src/algorithms/time_evolution/apply_gate.jl
+++ b/src/algorithms/time_evolution/apply_gate.jl
@@ -45,7 +45,7 @@ function _apply_gate(a::MPSTensor, b::MPSTensor, gate::NNGate, trunc::Truncation
     a, s, b, ϵ = svd_trunc!(a2b2; trunc)
     a, b = absorb_s(a, s, b)
     if need_flip
-        a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1)
+        a, s, b = flip(a, numind(a)), _fliptwist_s!(s), flip(b, 1)
     end
     b = permute(b, ((1, 2), (3,)))
     return a, s, b, ϵ
diff --git a/src/algorithms/time_evolution/apply_mpo.jl b/src/algorithms/time_evolution/apply_mpo.jl
index 25fd47974..05e6e6dba 100644
--- a/src/algorithms/time_evolution/apply_mpo.jl
+++ b/src/algorithms/time_evolution/apply_mpo.jl
@@ -214,11 +214,34 @@ function _proj_from_RL(
     @assert isdual(domain(l, 1)) == isdual(codomain(l, 1)) == false
     rl = r * l
     u, s, vh, ϵ = svd_trunc!(rl; trunc)
+    return _proj_from_svd(r, l, u, s, vh, ϵ)
+end
+
+# Second half of `_proj_from_RL`, split off so that the decomposition can be
+# done for the whole cluster at once (see `bond_svds`).
+function _proj_from_svd(r::MPSBondTensor, l::MPSBondTensor, u, s, vh, ϵ)
     sinv = sdiag_pow(s, -1 / 2)
     Pa, Pb = l * vh' * sinv, sinv * u' * r
     return Pa, s, Pb, ϵ
 end
 
+"""
+    bond_svds(rls, truncs)
+
+Truncated SVD of every internal bond of a cluster, returning `(u, s, vh, ϵ)` per bond.
+
+Separate from [`_proj_from_RL`](@ref) so that *all* bonds of the cluster can be decomposed 
+in a single batched call.
+"""
+function bond_svds(rls::AbstractVector, truncs::AbstractVector)
+    return bond_svds(storagetype(eltype(rls)), rls, truncs)
+end
+function bond_svds(::Type, rls::AbstractVector, truncs::AbstractVector)
+    return map(rls, truncs) do rl, trunc
+        return svd_trunc!(rl; trunc)
+    end
+end
+
 """
 Given a cluster `Ms`, find all projectors `Pa`, `Pb`
 and Schmidt weights `wts` on internal bonds.
@@ -229,8 +252,14 @@ function _get_allprojs(
     N = length(Ms)
     Rs, Ls = _get_allRLs(Ms)
     @assert length(truncs) == N - 1
-    projs_errs = map(Rs, Ls, truncs) do R, L, trunc
-        return _proj_from_RL(R, L; trunc)
+    for (R, L) in zip(Rs, Ls)
+        @assert isdual(domain(R, 1)) == isdual(codomain(R, 1)) == false
+        @assert isdual(domain(L, 1)) == isdual(codomain(L, 1)) == false
+    end
+    # decompose every bond of the cluster in one go, then finish each projector locally
+    svds = bond_svds(map(*, Rs, Ls), truncs)
+    projs_errs = map(Rs, Ls, svds) do R, L, (u, s, vh, ϵ)
+        return _proj_from_svd(R, L, u, s, vh, ϵ)
     end
     Pas = map(Base.Fix2(getindex, 1), projs_errs)
     wts = map(Base.Fix2(getindex, 2), projs_errs)
@@ -299,7 +328,7 @@ function _apply_gatempo!(
     fusers = map(Iterators.drop(Ms, 1), Iterators.drop(gs, 1)) do M, g
         V1, V2 = space(M, 1), space(g, 1)
         @assert !isdual(V1) && !isdual(V2)
-        return isomorphism(fuse(V1, V2) ← V1 ⊗ V2)
+        return isomorphism(storagetype(T1), fuse(V1, V2) ← V1 ⊗ V2)
     end
     #= gate on codomain of PEPS
            -3                         -3                          -3
@@ -336,7 +365,7 @@ function _apply_gatempo!(
     fusers = map(Iterators.drop(Ms, 1), Iterators.drop(gs, 1)) do M, g
         V1, V2 = space(M, 1), space(g, 1)
         @assert !isdual(V1) && !isdual(V2)
-        return isomorphism(fuse(V1, V2) ← V1 ⊗ V2)
+        return isomorphism(storagetype(M), fuse(V1, V2) ← V1 ⊗ V2)
     end
     #= gate on codomain of PEPO (gate_ax = 1)
 
diff --git a/src/algorithms/time_evolution/gaugefix_su.jl b/src/algorithms/time_evolution/gaugefix_su.jl
index 13a2782d6..a373d773d 100644
--- a/src/algorithms/time_evolution/gaugefix_su.jl
+++ b/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)
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index cd9516239..046cce397 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -155,10 +155,38 @@ function su_iter(
     return state2, env2, ϵ
 end
 
+"""
+    check_su_state(psi, iter)
+
+Check that a simple-update step did not produce a degenerate state.
+
+Without this check, later CTMRG can "converge" immediately because a `NaN`
+objective compares equal to itself, reports `converged = true`, and the run finishes with
+a `NaN` energy and a suspiciously fast wall time.
+
+Only vector-space dimensions are inspected, never tensor data, so this is quick and inexpensive.
+"""
+function check_su_state(psi, iter)
+    for (idx, t) in pairs(unitcell(psi))
+        dim(space(t)) > 0 || throw(
+            ErrorException(
+                "simple update produced a degenerate state at iteration $iter: tensor $idx \
+                has an empty space ($(space(t)))."
+            )
+        )
+    end
+    return nothing
+end
+
 function Base.iterate(it::TimeEvolver{<:SimpleUpdate}, state = it.state)
     iter, t = state.iter, state.t
     (iter == it.nstep) && return nothing
-    psi, env, ϵ = su_iter(state.psi, it.circuit, it.alg, state.env)
+    storage = storagetype(state.psi)
+    psi, env, ϵ = with_alloc_cache(storage, :su, iter) do
+        su_iter(state.psi, it.circuit, it.alg, state.env)
+    end
+    psi, env = uncache(psi, storage), uncache(env, storage)
+    check_su_state(psi, iter + 1)
     # update internal state
     iter += 1
     t += it.dt
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index 80947b876..2b5b87673 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -1,6 +1,6 @@
 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
@@ -8,7 +8,7 @@ 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
@@ -51,7 +51,7 @@ function _su_iter!(
         _nn_bondrev(site1, site2)
     end
     for (wt, (bond, rev), flip) in zip(wts, bond_revs, flips)
-        wt_new = flip ? _fliptwist_s(wt) : wt
+        wt_new = flip ? _fliptwist_s!(wt) : wt
         wt_new = rev ? transpose(wt_new) : wt_new
         env[CartesianIndex(bond)] = normalize!(wt_new, Inf)
     end
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index c1dc32233..d6065d20a 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/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
diff --git a/src/algorithms/time_evolution/trotter_gate.jl b/src/algorithms/time_evolution/trotter_gate.jl
index c4a52bfd5..7d20f40ab 100644
--- a/src/algorithms/time_evolution/trotter_gate.jl
+++ b/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(TensorKit.promote(physicalspace(H, x2), left_virtualspace(gate[2]))[1])
+        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
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index 533fac679..b825b2f80 100644
--- a/src/algorithms/toolbox.jl
+++ b/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),
     )
 
diff --git a/src/algorithms/truncation/bond_truncation.jl b/src/algorithms/truncation/bond_truncation.jl
index 8fbb9d1ca..81983dfd0 100644
--- a/src/algorithms/truncation/bond_truncation.jl
+++ b/src/algorithms/truncation/bond_truncation.jl
@@ -141,7 +141,7 @@ function bond_truncate(a::MPSTensor, b::MPSTensor, benv::BondEnv, alg::ALSTrunca
     a, b = absorb_s(a, s, b)
     b = permute(b, ((1, 2), (3,)))
     if need_flip
-        a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1)
+        a, s, b = flip(a, numind(a)), _fliptwist_s!(s), flip(b, 1)
     end
     return a, s, b, (; fid, Δfid, Δs)
 end
@@ -182,7 +182,7 @@ function bond_truncate(a::MPSTensor, b::MPSTensor, benv::BondEnv, alg::FullEnvTr
     @tensor a[-1 -2; -3] := Qa[-1 -2 3] * u[3 -3]
     @tensor b[-1 -2; -3] := vh[-1 1] * Qb[1 -2 -3]
     if need_flip
-        a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1)
+        a, s, b = flip(a, numind(a)), _fliptwist_s!(s), flip(b, 1)
     end
     return a, s, b, info
 end
diff --git a/src/environments/bp_environments.jl b/src/environments/bp_environments.jl
index bc0429e0a..53785faf4 100644
--- a/src/environments/bp_environments.jl
+++ b/src/environments/bp_environments.jl
@@ -27,6 +27,8 @@ struct BPEnv{T}
     "4 x rows x cols array of message tensors, where the first dimension specifies the spatial direction"
     messages::Array{T, 3}
 end
+TensorKit.storagetype(::Type{BPEnv{T}}) where {T} = storagetype(T)
+
 
 """
 Construct a message tensor on a certain bond of a network,
@@ -118,13 +120,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
@@ -172,7 +174,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
diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl
index d6f39aa7d..35f080fbe 100644
--- a/src/environments/ctmrg_environments.jl
+++ b/src/environments/ctmrg_environments.jl
@@ -61,7 +61,7 @@ end
 
 """
     CTMRGEnv(
-        [f=randn, T=ComplexF64], Ds_north::A, Ds_east::A, chis_north::B, [chis_east::B], [chis_south::B], [chis_west::B]
+    [f=randn, T=ComplexF64, TA=Matrix{ComplexF64},] Ds_north::A, Ds_east::A, chis_north::B, [chis_east::B], [chis_south::B], [chis_west::B]
     ) where {A<:AbstractMatrix{<:VectorSpace}, B<:AbstractMatrix{<:ElementarySpace}}
 
 Construct a CTMRG environment by specifying matrices of north and east virtual spaces of the
@@ -82,10 +82,11 @@ of a partition function defined in terms of local rank-4 tensors) or a `ProductS
 for the case of a network representing overlaps of PEPSs and PEPOs).
 """
 function CTMRGEnv(
-        f, T, Ds_north::A, Ds_east::A, chis_north::B, chis_east::B = chis_north,
+        f, ::Type{T}, Ds_north::A, Ds_east::A, chis_north::B, chis_east::B = chis_north,
         chis_south::B = chis_north, chis_west::B = chis_north,
     ) where {
         A <: AbstractMatrix{<:ProductSpace}, B <: AbstractMatrix{<:ElementarySpace},
+        T,
     }
     # check all of the sizes
     size(Ds_north) == size(Ds_east) == size(chis_north) == size(chis_east) ==
@@ -100,7 +101,6 @@ function CTMRGEnv(
     st = spacetype(first(Ds_north))
     C_type = tensormaptype(st, 1, 1, T)
     T_type = tensormaptype(st, N + 1, 1, T)
-
     # First index is direction
     corners = Array{C_type}(undef, 4, size(Ds_north)...)
     edges = Array{T_type}(undef, 4, size(Ds_north)...)
@@ -178,9 +178,9 @@ The environment virtual spaces for each site correspond to virtual space of the
 corresponding edge tensor for each direction.
 """
 function CTMRGEnv(
-        f, T,
+        f, ::Type{T},
         D_north::S, D_east::S, virtual_spaces...; unitcell::Tuple{Int, Int} = (1, 1),
-    ) where {S <: VectorSpace}
+    ) where {S <: VectorSpace, T}
     return CTMRGEnv(
         f, T,
         _fill_edge_physical_spaces(D_north, D_east; unitcell)...,
@@ -215,19 +215,22 @@ of the corresponding edge tensor for each direction. Specifically, for a given s
 `chis_south[r, c]` corresponds to the east space of the south edge tensor, and
 `chis_west[r, c]` corresponds to the north space of the west edge tensor.
 """
-function CTMRGEnv(f, T, network::InfiniteSquareNetwork, virtual_spaces...)
+function CTMRGEnv(f, ::Type{T}, network::N, virtual_spaces...) where {T, N <: InfiniteSquareNetwork}
     Ds_north = _north_edge_physical_spaces(network)
     Ds_east = _east_edge_physical_spaces(network)
     virtual_spaces = _fill_environment_virtual_spaces(virtual_spaces...; unitcell = size(network))
-    return CTMRGEnv(f, T, Ds_north, Ds_east, virtual_spaces...)
+    return CTMRGEnv(f, similarstoragetype(storagetype(network), eltype(T)), Ds_north, Ds_east, virtual_spaces...)
 end
-function CTMRGEnv(network::Union{InfiniteSquareNetwork, InfinitePartitionFunction, InfinitePEPS}, virtual_spaces...)
-    return CTMRGEnv(randn, scalartype(network), network, virtual_spaces...)
+function CTMRGEnv(network::InfiniteSquareNetwork{O}, virtual_spaces...) where {O}
+    return CTMRGEnv(randn, storagetype(O), network, virtual_spaces...)
+end
+function CTMRGEnv(network::Union{<:InfinitePartitionFunction{T}, <:InfinitePEPS{T}}, virtual_spaces...) where {T}
+    return CTMRGEnv(randn, storagetype(T), network, virtual_spaces...)
 end
 
 # allow constructing environments for implicitly defined contractible networks
-function CTMRGEnv(f, T, state::Union{InfinitePartitionFunction, InfinitePEPS}, args...)
-    return CTMRGEnv(f, T, InfiniteSquareNetwork(state), args...)
+function CTMRGEnv(f, ::Type{T}, state::Union{InfinitePartitionFunction, InfinitePEPS}, args...) where {T}
+    return CTMRGEnv(f, similarstoragetype(eltype(state), eltype(T)), InfiniteSquareNetwork(state), args...)
 end
 
 # copy-like constructor
@@ -235,6 +238,8 @@ CTMRGEnv(env::CTMRGEnv) = CTMRGEnv(env.corners, env.edges)
 
 @non_differentiable CTMRGEnv(state::Union{InfinitePartitionFunction, InfinitePEPS}, args...)
 
+TensorKit.storagetype(::Type{CTMRGEnv{C, E}}) where {C, E} = storagetype(C) == storagetype(E) ? storagetype(C) : promote_type(storagetype(C), storagetype(E))
+
 # Custom adjoint for CTMRGEnv constructor, needed for fixed-point differentiation
 function ChainRulesCore.rrule(
         ::Type{CTMRGEnv}, corners::Array{C, 3}, edges::Array{T, 3}
@@ -249,17 +254,18 @@ function ChainRulesCore.rrule(::typeof(getproperty), e::CTMRGEnv, name::Symbol)
     if name === :corners
         function corner_pullback(Δcorners_)
             Δcorners = unthunk(Δcorners_)
-            return NoTangent(), CTMRGEnv(Δcorners, zerovector.(e.edges)), NoTangent()
+            zvs = CTMRGEnv(Δcorners, zerovector.(e.edges))
+            return NoTangent(), zvs, NoTangent()
         end
         return result, corner_pullback
     elseif name === :edges
         function edge_pullback(Δedges_)
             Δedges = unthunk(Δedges_)
-            return NoTangent(), CTMRGEnv(zerovector.(e.corners), Δedges), NoTangent()
+            zvs = CTMRGEnv(zerovector.(e.corners), Δedges)
+            return NoTangent(), zvs, NoTangent()
         end
         return result, edge_pullback
-    else
-        # this should never happen because already errored in forwards pass
+    else # this should never happen because already errored in forwards pass
         throw(ArgumentError("No rrule for getproperty of $name"))
     end
 end
diff --git a/src/environments/product_state_environments.jl b/src/environments/product_state_environments.jl
index 394bc5751..93a1577b9 100644
--- a/src/environments/product_state_environments.jl
+++ b/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
diff --git a/src/environments/suweight.jl b/src/environments/suweight.jl
index dbf7b0059..3a91e4098 100644
--- a/src/environments/suweight.jl
+++ b/src/environments/suweight.jl
@@ -63,13 +63,21 @@ end
 Create a trivial `SUWeight` by specifying the vertical (north) or horizontal (east) virtual bond spaces.
 """
 function SUWeight(
+        ::Type{TorA},
         Nspaces::M, Espaces::M = Nspaces
-    ) where {M <: AbstractMatrix{<:ElementarySpace}}
+    ) where {M <: AbstractMatrix{<:ElementarySpace}, TorA}
     @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)
+        if TorA <: AbstractArray
+            realTorA = similarstoragetype(TorA, real(eltype(TorA)))
+            diag = realTorA(undef, reduceddim(V))
+            fill!(diag, 1)
+        else
+            diag = ones(real(TorA), reduceddim(V))
+        end
+        DiagonalTensorMap(diag, V)
     end
     return SUWeight(weights)
 end
@@ -81,9 +89,10 @@ Create a trivial `SUWeight` by specifying its vertical (north) and horizontal (e
 as `ElementarySpace`s) and unit cell size.
 """
 function SUWeight(
+        ::Type{TorA},
         Nspace::S, Espace::S = Nspace; unitcell::Tuple{Int, Int} = (1, 1)
-    ) where {S <: ElementarySpace}
-    return SUWeight(fill(Nspace, unitcell), fill(Espace, unitcell))
+    ) where {S <: ElementarySpace, TorA}
+    return SUWeight(TorA, fill(Nspace, unitcell), fill(Espace, unitcell))
 end
 
 """
@@ -94,7 +103,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 +115,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)
@@ -140,6 +149,8 @@ TensorKit.spacetype(::Type{T}) where {E, T <: SUWeight{E}} = spacetype(E)
 TensorKit.sectortype(w::SUWeight) = sectortype(typeof(w))
 TensorKit.sectortype(::Type{<:SUWeight{T}}) where {T} = sectortype(spacetype(T))
 
+TensorKit.storagetype(::Type{SUWeight{T}}) where {T} = storagetype(T)
+
 ## Bipartite check
 function _is_bipartite(wts::SUWeight)
     (size(wts, 2) == size(wts, 3) == 2) || (return false)
@@ -328,7 +339,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
diff --git a/src/environments/vumps_environments.jl b/src/environments/vumps_environments.jl
index 48605f2f1..6ef0375d0 100644
--- a/src/environments/vumps_environments.jl
+++ b/src/environments/vumps_environments.jl
@@ -25,10 +25,11 @@ 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
 
@@ -36,7 +37,8 @@ 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
diff --git a/src/networks/infinitesquarenetwork.jl b/src/networks/infinitesquarenetwork.jl
index fe62377b3..46161327b 100644
--- a/src/networks/infinitesquarenetwork.jl
+++ b/src/networks/infinitesquarenetwork.jl
@@ -25,6 +25,7 @@ struct InfiniteSquareNetwork{O}
     end
 end
 InfiniteSquareNetwork(n::InfiniteSquareNetwork) = n
+TensorKit.storagetype(::Type{InfiniteSquareNetwork{O}}) where {O} = storagetype(O)
 
 ## Unit cell interface
 
diff --git a/src/networks/local_sandwich.jl b/src/networks/local_sandwich.jl
index d29f54ea1..4ce94e850 100644
--- a/src/networks/local_sandwich.jl
+++ b/src/networks/local_sandwich.jl
@@ -54,6 +54,8 @@ _isapprox_localsandwich(O1::PFTensor, O2::PFTensor; kwargs...) = isapprox(O1, O2
 ## PEPS
 
 const PEPSSandwich{T <: PEPSTensor} = Tuple{T, T}
+TensorKit.storagetype(::Type{PEPSSandwich{T}}) where {T} = T
+TensorKit.storagetype(S::PEPSSandwich{T}) where {T} = T
 
 ket(O::PEPSSandwich) = O[1]
 bra(O::PEPSSandwich) = O[2]
diff --git a/src/networks/tensors.jl b/src/networks/tensors.jl
index 6ac6411da..d014150b4 100644
--- a/src/networks/tensors.jl
+++ b/src/networks/tensors.jl
@@ -23,17 +23,17 @@ const PartitionFunctionTensor{S <: ElementarySpace} = AbstractTensorMap{<:Any, S
 const PFTensor = PartitionFunctionTensor
 
 """
-    PartitionFunctionTensor(f, ::Type{T}, Pspace::S, Nspace::S,
-               [Espace::S], [Sspace::S], [Wspace::S]) where {T,S<:Union{Int,ElementarySpace}}
+    PartitionFunctionTensor(f, ::Type{TorA}, Pspace::S, Nspace::S,
+               [Espace::S], [Sspace::S], [Wspace::S]) where {TorA, S<:Union{Int,ElementarySpace}}
                 
-Construct a PartitionFunctionTensor tensor based on the north, east, west and south spaces.
+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`.
 """
 function PartitionFunctionTensor(
-        f, ::Type{T},
+        f, ::Type{TorA},
         Nspace::S, Espace::S = Nspace, Sspace::S = Nspace, Wspace::S = Espace,
-    ) where {T, S <: ElementarySpace}
-    return f(T, Wspace ⊗ Sspace ← Nspace ⊗ Espace)
+    ) where {TorA, S <: ElementarySpace}
+    return f(TorA, Wspace ⊗ Sspace ← Nspace ⊗ Espace)
 end
 
 Base.rotl90(t::PFTensor) = permute(t, ((3, 1), (4, 2)))
@@ -79,18 +79,19 @@ respectively.
 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}}
+    PEPSTensor(f, ::Type{TorA}, Pspace::S, Nspace::S,
+               [Espace::S], [Sspace::S], [Wspace::S]) where {TorA, 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`.
 """
 function PEPSTensor(
-        f, ::Type{T},
+        f,
+        ::Type{TorA},
         Pspace::S,
         Nspace::S, Espace::S = Nspace, Sspace::S = Nspace', Wspace::S = Espace',
-    ) where {T, S <: ElementarySpace}
-    return f(T, Pspace ← Nspace ⊗ Espace ⊗ Sspace ⊗ Wspace)
+    ) where {TorA, S <: ElementarySpace}
+    return f(TorA, Pspace ← Nspace ⊗ Espace ⊗ Sspace ⊗ Wspace)
 end
 
 Base.rotl90(t::PEPSTensor) = permute(t, ((1,), (3, 4, 5, 2)))
@@ -155,7 +156,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
 
diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl
index f902a3874..65ec5194d 100644
--- a/src/operators/infinitepepo.jl
+++ b/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
@@ -169,6 +169,8 @@ function physicalspace(T::InfinitePEPO, r::Int, c::Int)
     return codomain_physicalspace(T, r, c)
 end
 
+TensorKit.storagetype(::Type{InfinitePEPO{T}}) where {T} = storagetype(T)
+
 ## InfiniteSquareNetwork interface
 
 function InfiniteSquareNetwork(top::InfinitePEPS, mid::InfinitePEPO, bot::InfinitePEPS = top)
diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl
index c35e5ddb9..b38ee1dc0 100644
--- a/src/operators/localoperator.jl
+++ b/src/operators/localoperator.jl
@@ -86,6 +86,7 @@ function add_term!(
 
     return operator
 end
+TensorKit.storagetype(lo::LocalOperator{T, S}) where {T, S} = storagetype(first(lo.terms)[2]) # horrible!
 
 
 """
diff --git a/src/operators/transfermatrix.jl b/src/operators/transfermatrix.jl
index 76cfc54e3..e2724b25b 100644
--- a/src/operators/transfermatrix.jl
+++ b/src/operators/transfermatrix.jl
@@ -14,6 +14,8 @@ function which corresponds to the overlap between 'ket' and 'bra' `InfinitePEPS`
 """
 const InfiniteTransferPEPS{T <: PEPSTensor} = InfiniteMPO{PEPSSandwich{T}}
 
+TensorKit.storagetype(::InfiniteTransferPEPS{T}) where {T} = storagetype(T)
+
 function InfiniteTransferPEPS(
         top::PeriodicArray{T, 1}, bot::PeriodicArray{T, 1}
     ) where {T <: PEPSTensor}
@@ -69,6 +71,8 @@ function which corresponds to the expectation value of an `InfinitePEPO` between
 """
 const InfiniteTransferPEPO{H, T <: PEPSTensor, O <: PEPOTensor} = InfiniteMPO{PEPOSandwich{H, T, O}}
 
+TensorKit.storagetype(::InfiniteTransferPEPO{H, T, O}) where {H, T, O} = storagetype(T)
+
 function InfiniteTransferPEPO(
         top::PeriodicArray{T, 1}, mid::PeriodicArray{O, 2}, bot::PeriodicArray{T, 1}
     ) where {T, O}
@@ -127,13 +131,13 @@ virtualspace(O::InfiniteTransferMatrix, i, dir) = virtualspace(O[i], dir)
 """
     initialize_mps(
         f=randn,
-        T=scalartype(O),
+        T=storagetype(O),
         O::Union{InfiniteTransferPEPS,InfiniteTransferPEPO},
         virtualspaces::AbstractArray{<:ElementarySpace,1}
     )
     initialize_mps(
         f=randn,
-        T=scalartype(O),
+        T=storagetype(O),
         O::Union{MultilineTransferPEPS,MultilineTransferPEPO},
         virtualspaces::AbstractArray{<:ElementarySpace,2}
     )
@@ -141,37 +145,37 @@ virtualspace(O::InfiniteTransferMatrix, i, dir) = virtualspace(O[i], dir)
 Inialize a boundary MPS for the transfer operator `O` by specifying an array of virtual
 spaces consistent with the unit cell.
 """
-function initialize_mps(O::Union{InfiniteTransferMatrix, MultilineTransferMatrix}, arg) # initialize(f=randn, T=scalartype(O), O, ...)
-    return initialize_mps(randn, scalartype(O), O, arg)
+function initialize_mps(O::Union{InfiniteTransferMatrix, MultilineTransferMatrix}, arg; kwargs...) # initialize(f=randn, T=scalartype(O), O, ...)
+    return initialize_mps(randn, storagetype(O), O, arg; kwargs...)
 end
 function initialize_mps(
-        f, T, O::InfiniteTransferMatrix, virtualspaces::AbstractArray{S, 1}
-    ) where {S}
+        f, ::Type{TorA}, O::InfiniteTransferMatrix, virtualspaces::AbstractArray{S, 1}; kwargs...
+    ) where {S, TorA}
     return InfiniteMPS(
         [
             f(
-                T,
+                TorA,
                 virtualspaces[_prev(i, end)] * _elementwise_dual(north_virtualspace(O, i)),
                 virtualspaces[mod1(i, end)],
             ) for i in 1:length(O)
-        ]
+        ]; kwargs...
     )
 end
 function initialize_mps(
-        f, T, O::MultilineTransferMatrix, virtualspaces::AbstractArray{S, 2}
-    ) where {S}
+        f, ::Type{TorA}, O::MultilineTransferMatrix, virtualspaces::AbstractArray{S, 2}; kwargs...
+    ) where {S, TorA}
     mpss = map(1:size(O, 1)) do r
-        return initialize_mps(f, T, O[r], virtualspaces[r, :])
+        return initialize_mps(f, TorA, O[r], virtualspaces[r, :]; kwargs...)
     end
     return MPSKit.Multiline(mpss)
 end
 function initialize_mps(
-        f, T, O::MultilineTransferMatrix, virtualspaces::AbstractArray{S, 1}
-    ) where {S}
-    return initialize_mps(f, T, O, repeat(virtualspaces, length(O), 1))
+        f, ::Type{TorA}, O::MultilineTransferMatrix, virtualspaces::AbstractArray{S, 1}; kwargs...
+    ) where {S, TorA}
+    return initialize_mps(f, TorA, O, repeat(virtualspaces, length(O), 1); kwargs...)
 end
-function initialize_mps(f, T, O::MultilineTransferMatrix, V::ElementarySpace)
-    return initialize_mps(f, T, O, repeat([V], length(O), length(O[1])))
+function initialize_mps(f, ::Type{TorA}, O::MultilineTransferMatrix, V::ElementarySpace; kwargs...) where {TorA}
+    return initialize_mps(f, TorA, O, repeat([V], length(O), length(O[1])); kwargs...)
 end
 
 @doc """
diff --git a/src/states/infinitepartitionfunction.jl b/src/states/infinitepartitionfunction.jl
index 392668f21..9104aa86a 100644
--- a/src/states/infinitepartitionfunction.jl
+++ b/src/states/infinitepartitionfunction.jl
@@ -25,6 +25,7 @@ struct InfinitePartitionFunction{T <: PartitionFunctionTensor}
         return new{T}(A)
     end
 end
+TensorKit.storagetype(::Type{InfinitePartitionFunction{T}}) where {T} = storagetype(T)
 
 const InfinitePF{T} = InfinitePartitionFunction{T}
 
@@ -50,8 +51,8 @@ of the PEPS tensor at each site in the unit cell as a matrix. Each individual sp
 specified as either an `Int` or an `ElementarySpace`.
 """
 function InfinitePartitionFunction(
-        f, T, Nspaces::M, Espaces::M = Nspaces
-    ) where {M <: AbstractMatrix{<:ElementarySpace}}
+        f, ::Type{T}, Nspaces::M, Espaces::M = Nspaces
+    ) where {M <: AbstractMatrix{<:ElementarySpace}, T <: Number}
     size(Nspaces) == size(Espaces) ||
         throw(ArgumentError("Input spaces should have equal sizes."))
 
diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl
index 3d3a82f0f..1fe16e5e8 100644
--- a/src/states/infinitepeps.jl
+++ b/src/states/infinitepeps.jl
@@ -46,8 +46,8 @@ Create an `InfinitePEPS` by specifying the physical, north virtual and east virt
 of the PEPS tensor at each site in the unit cell as a matrix.
 """
 function InfinitePEPS(
-        f, T::Type{<:Number}, Pspaces::M, Nspaces::M, Espaces::M = Nspaces
-    ) where {M <: AbstractMatrix{<:ElementarySpace}}
+        f, ::Type{TorA}, Pspaces::M, Nspaces::M, Espaces::M = Nspaces
+    ) where {M <: AbstractMatrix{<:ElementarySpace}, TorA}
     size(Pspaces) == size(Nspaces) == size(Espaces) ||
         throw(ArgumentError("Input spaces should have equal sizes."))
 
@@ -55,7 +55,7 @@ function InfinitePEPS(
     Wspaces = adjoint.(circshift(Espaces, (0, 1)))
 
     A = map(Pspaces, Nspaces, Espaces, Sspaces, Wspaces) do P, N, E, S, W
-        return PEPSTensor(f, T, P, N, E, S, W)
+        return PEPSTensor(f, TorA, P, N, E, S, W)
     end
 
     return InfinitePEPS(A)
@@ -63,8 +63,9 @@ end
 function InfinitePEPS(
         Pspaces::A, virtual_spaces...; kwargs...
     ) where {A <: Union{AbstractMatrix{<:ElementarySpace}, ElementarySpace}}
-    return InfinitePEPS(randn, ComplexF64, Pspaces, virtual_spaces...; kwargs...)
+    return InfinitePEPS(randn, Vector{ComplexF64}, Pspaces, virtual_spaces...; kwargs...)
 end
+TensorKit.storagetype(::Type{InfinitePEPS{T}}) where {T} = storagetype(T)
 
 """
     InfinitePEPS(A::PEPSTensor; unitcell=(1, 1))
@@ -108,8 +109,8 @@ end
 Create an InfinitePEPS by specifying its physical, north and east spaces and unit cell.
 """
 function InfinitePEPS(
-        f, T::Type{<:Number}, Pspace::S, vspaces...; unitcell::Tuple{Int, Int} = (1, 1)
-    ) where {S <: ElementarySpace}
+        f, ::Type{T}, Pspace::S, vspaces...; unitcell::Tuple{Int, Int} = (1, 1)
+    ) where {S <: ElementarySpace, T}
     return InfinitePEPS(
         f, T,
         _fill_state_physical_spaces(Pspace; unitcell),
@@ -126,7 +127,7 @@ Base.eltype(::Type{InfinitePEPS{T}}) where {T} = T
 Base.eltype(A::InfinitePEPS) = eltype(typeof(A))
 
 Base.copy(A::InfinitePEPS) = InfinitePEPS(copy(unitcell(A)))
-function Base.similar(A::InfinitePEPS, T::Type{TorA} = scalartype(A)) where {TorA}
+function Base.similar(A::InfinitePEPS, T::Type = scalartype(A))
     return InfinitePEPS(map(t -> similar(t, T), unitcell(A)))
 end
 Base.repeat(A::InfinitePEPS, counts...) = InfinitePEPS(repeat(unitcell(A), counts...))
diff --git a/src/utility/alloc_cache.jl b/src/utility/alloc_cache.jl
new file mode 100644
index 000000000..d17b45abe
--- /dev/null
+++ b/src/utility/alloc_cache.jl
@@ -0,0 +1,131 @@
+"""
+    with_alloc_cache(f, storage, caller::Symbol, iter::Int) -> f()
+
+Run one iteration `f()` of an iterative algorithm using a memory cache.
+
+**This is a no-op unless `GPUArrays` is loaded _and_ `storage` is a GPU array.**
+
+There seems to be no benefit to caching for CPU-side memory, but for the device,
+using a warm memory pool which recycles memory avoids everything being blocked
+while the GPU device driver allocates.
+
+`caller` ids the calling algo (`:su`, `:ctmrg`) so that algorithms allocating different
+buffer sizes avoid stepping on each others' caches. Not every algorithm needs an allocation cache,
+only the ones with repeated iterations.
+
+Belief propagation seems to not benefit from the caching as much, so it's currently unused there.
+
+`iter` selects between `depth` alternating caches per `caller`. A buffer allocated during
+iteration `i` can't become reusable until later iterations have completely used and discarded
+its output, so with too few caches it would be possible to overwrite state that is still in
+use. A caller that copies its result out with [`uncache`](@ref) before the next iteration
+begins needs only a single cache; simple update keeps 2, since it hands its result off after
+the following step.
+
+!!! note
+    Buffers are only returned to the cache when the enclosing block exits, so a single
+    block retains *everything* it allocated rather than just its maximum live block.
+    Keeping more caches in the rotation multiplies this effect, which matters for
+    algorithms like (esp. sequential) CTMRG whose iterations allocate many
+    short-lived temporaries.
+
+Caching is skipped while `Zygote.jl` is differentiating, because the reverse-mode tape holds references
+to intermediates, and recycling those could silently corrupt gradients.
+"""
+function with_alloc_cache(f, storage::Type, caller::Symbol, iter::Int, depth::Int = 2)
+    Zygote.isderiving() && return f()
+    return _with_alloc_cache(f, storage, caller, iter, depth)
+end
+
+"""
+    uncache(x, storage::Type) -> x
+
+Copy `x` out of any allocation-cache region, so that it stays live and is not overwriten
+once the enclosing [`with_alloc_cache`](@ref) block has exited. This is necessary to
+ensure `leading_boundary` and other functions which call back into AD handle cached
+memory correctly.
+
+Buffers allocated inside a cache block are handed back to the pool when the block exits,
+and the next cached call may hand them out again, which *silently overwrites* a result the
+caller is still holding. Anything that escapes such a block must be copied out.
+
+**This is a no-op unless `GPUArrays` is loaded _and_ `storage` is a GPU array**, and also
+while Zygote is differentiating, since caching is skipped in both of those cases.
+"""
+function uncache(x, storage::Type)
+    Zygote.isderiving() && return x
+    return _uncache(x, storage)
+end
+_uncache(x, ::Type) = x
+
+"""
+    alloc_cache_depth(alg)
+
+How many iterations a buffer must go unused for before it may be recycled.
+"""
+alloc_cache_depth(alg) = 1
+_with_alloc_cache(f, ::Type, ::Symbol, ::Int, ::Int) = f()
+
+"""
+    free_alloc_caches!(storage)
+    free_alloc_caches!(storage, caller::Symbol)
+
+Release memory held by the allocation caches for storage type `storage`, either for every
+`caller` (if none is provided) or only for the given one. Does nothing unless the
+`PEPSKitGPUArraysExt` extension is loaded and `storage` is a GPU array type.
+This should be called when the bond dimension changes, because the cache keys depend on
+buffer size and the caches can't be reused when the bond dimension has changed.
+
+!!! warning
+    This releases the underlying device memory, so it is only safe to call when nothing
+    still references a buffer that was allocated inside a [`with_alloc_cache`](@ref) block.
+    Results that escape such a block must have been copied out with [`uncache`](@ref)
+    first, otherwise freeing the cache leaves them undefined.
+"""
+free_alloc_caches!(::Type) = nothing
+free_alloc_caches!(::Type, ::Symbol) = nothing
+
+# last buffer-shape signature seen per caller, used to detect when cached buffer sizes
+# have gone stale because a bond dimension changed
+const ALLOC_CACHE_SIGNATURES = Dict{Tuple{Symbol, Symbol}, UInt}()
+const ALLOC_CACHE_SIGNATURES_LOCK = ReentrantLock()
+
+"""
+    free_stale_alloc_caches!(storage, caller::Symbol, phase::Symbol, signature::UInt)
+
+Release `caller`'s allocation caches when `signature` differs from the one seen at the same
+`phase` of the previous call, and record `signature` as the current one for that phase.
+
+The caches are keyed by buffer size, so that when the size changes the old, unusable caches
+can be freed and new ones allocated, corresponding to the new size. This avoids the cache
+size growing without bound.
+
+`phase` distinguishes the points a caller checks from. For example, CTMRG grows the
+corner and edge spaces as it converges, so its incoming and outgoing environments differ
+whenever the truncation is *not* fixed-space. Recording both under one key makes the stored
+signature alternate between them, so every check reports stale and the pool is freed on every
+call. Comparing each phase only against itself keeps the pool across repeated calls while
+still invalidating it when the spaces genuinely change.
+
+Like [`free_alloc_caches!`](@ref) this is a no-op without a GPU storage type, and it's
+skipped while `Zygote.jl` is differentiating, since caching is disabled there anyway.
+"""
+function free_stale_alloc_caches!(storage::Type, caller::Symbol, phase::Symbol, signature::UInt)
+    Zygote.isderiving() && return nothing
+    stale = Base.@lock ALLOC_CACHE_SIGNATURES_LOCK begin
+        key = (caller, phase)
+        previous = get(ALLOC_CACHE_SIGNATURES, key, nothing)
+        ALLOC_CACHE_SIGNATURES[key] = signature
+        !isnothing(previous) && previous != signature
+    end
+    stale && free_alloc_caches!(storage, caller)
+    return nothing
+end
+
+"""
+    alloc_cache_signature(tensors) -> UInt
+
+Hash the spaces of `tensors`, for use as a [`free_stale_alloc_caches!`](@ref) signature.
+Two states whose tensors live in the same spaces allocate the same buffer sizes.
+"""
+alloc_cache_signature(tensors) = hash(map(space, tensors))
diff --git a/src/utility/eigh.jl b/src/utility/eigh.jl
index 0db3bb7ef..e4eba3a61 100644
--- a/src/utility/eigh.jl
+++ b/src/utility/eigh.jl
@@ -224,6 +224,7 @@ function _compute_eighdata!(
     I = sectortype(f)
     dims = SectorDict{I, Int}()
 
+    Dtype = similarstoragetype(f, real(scalartype(f)))
     sectors = trunc isa NoTruncation ? blocksectors(f) : blocksectors(trunc.space)
     generator = Base.Iterators.map(sectors) do c
         b = block(f, c)
@@ -233,7 +234,7 @@ function _compute_eighdata!(
             D, V = eigh_full!(b)
             lm_ordering = sortperm(abs.(D.diag); rev = true) # order values and vectors consistently with eigsolve
             D = D.diag[lm_ordering] # extracts diagonal as Vector instead of Diagonal to make compatible with D of svdsolve
-            V = stack(eachcol(V)[lm_ordering])[:, 1:howmany]
+            V = V[:, view(lm_ordering, 1:howmany)]
         else
             x₀ = alg.start_vector(b)
             eig_alg = alg.alg
@@ -246,18 +247,17 @@ function _compute_eighdata!(
                 D, V = eigh_full!(b)
                 lm_ordering = sortperm(abs.(D.diag); rev = true)
                 D = D.diag[lm_ordering]
-                V = stack(eachcol(V)[lm_ordering])[:, 1:howmany]
+                V = V[:, view(lm_ordering, 1:howmany)]
             else  # Slice in case more values were converged than requested
                 V = stack(view(lvecs, 1:howmany))
             end
         end
-
         # make it deterministic-ish
         MatrixAlgebraKit.gaugefix!(eigh_full!, V)
 
         resize!(D, howmany)
         dims[c] = length(D)
-        return c => (D, V)
+        return c => (Dtype(D), V)
     end
 
     eigdata = SectorDict(generator)
@@ -314,7 +314,7 @@ function ChainRulesCore.rrule(
 
     function eigh_trunc!_full_pullback(ΔDV)
         Δt = eigh_pullback!(
-            zeros(scalartype(t), space(t)), t, (D, V), ΔDV, inds;
+            zeros(storagetype(t), space(t)), t, (D, V), ΔDV, inds;
             gauge_atol = gtol(ΔDV), degeneracy_atol = alg.rrule_alg.degeneracy_atol,
         )
         return NoTangent(), Δt, NoTangent()
@@ -338,7 +338,7 @@ function ChainRulesCore.rrule(
 
     function eigh_trunc!_trunc_pullback(ΔDV)
         Δf = eigh_trunc_pullback!(
-            zeros(scalartype(t), space(t)), t, (D, V), ΔDV;
+            zeros(storagetype(t), space(t)), t, (D, V), ΔDV;
             gauge_atol = gtol(ΔDV), degeneracy_atol = alg.rrule_alg.degeneracy_atol,
         )
         return NoTangent(), Δf, NoTangent()
diff --git a/src/utility/qr.jl b/src/utility/qr.jl
index a512821c8..bb67b43f8 100644
--- a/src/utility/qr.jl
+++ b/src/utility/qr.jl
@@ -96,7 +96,7 @@ function ChainRulesCore.rrule(
     gtol = _get_pullback_gauge_tol(alg.rrule_alg.verbosity)
 
     function left_orth!_pullback(ΔQR)
-        Δt = zeros(scalartype(t), space(t))
+        Δt = zeros(storagetype(t), space(t))
         MatrixAlgebraKit.qr_pullback!(Δt, t, QR, unthunk.(ΔQR); gauge_atol = gtol(ΔQR))
         return NoTangent(), Δt, NoTangent()
     end
diff --git a/src/utility/svd.jl b/src/utility/svd.jl
index 8e17913fe..9dd630cd8 100644
--- a/src/utility/svd.jl
+++ b/src/utility/svd.jl
@@ -186,7 +186,7 @@ end
 _default_svd_rrule_alg(::IterSVD) = :TruncPullback
 
 random_start_vector(t::AbstractMatrix) = randn(scalartype(t), size(t, 1))
-deterministic_start_vector(t::AbstractMatrix) = ones(scalartype(t), size(t, 1))
+deterministic_start_vector(t::AbstractMatrix) = fill!(similar(t, size(t, 1)), one(scalartype(t)))
 
 # Compute SVD data block-wise using KrylovKit algorithm
 # TODO: redefine _empty_svdtensors, _create_svdtensors
@@ -206,7 +206,7 @@ end
 function MatrixAlgebraKit.svd_trunc!(f, alg::TruncatedAlgorithm{<:IterSVD})
     U, S, Vᴴ = svd_trunc_no_error!(f, alg)
     truncation_error =
-        (trunc isa NoTruncation || isempty(blocksectors(f))) ? abs(zero(scalartype(f))) : norm(U * S * Vᴴ - f)
+        (alg.trunc isa NoTruncation || isempty(blocksectors(f))) ? abs(zero(scalartype(f))) : norm(U * S * Vᴴ - f)
     return U, S, Vᴴ, truncation_error
 end
 
@@ -239,6 +239,7 @@ function _compute_svddata!(
     I = sectortype(f)
     dims = SectorDict{I, Int}()
 
+    Stype = similarstoragetype(f, real(scalartype(f)))
     sectors = trunc isa NoTruncation ? blocksectors(f) : blocksectors(trunc.space)
     generator = Base.Iterators.map(sectors) do c
         b = block(f, c)
@@ -273,7 +274,7 @@ function _compute_svddata!(
 
         resize!(S, howmany)
         dims[c] = length(S)
-        return c => (U, S, V)
+        return c => (U, Stype(S), V)
     end
 
     SVDdata = SectorDict(generator)
@@ -303,7 +304,7 @@ function ChainRulesCore.rrule(
     function svd_trunc!_full_pullback(ΔUSV′)
         ΔUSV = unthunk.(ΔUSV′)
         Δt = svd_pullback!(
-            zeros(scalartype(t), space(t)), t, (U, S, V⁺), ΔUSV, inds;
+            zeros(storagetype(t), space(t)), t, (U, S, V⁺), ΔUSV, inds;
             gauge_atol = gtol(ΔUSV), degeneracy_atol = alg.rrule_alg.degeneracy_atol,
         )
         return NoTangent(), Δt, NoTangent()
@@ -333,7 +334,7 @@ function ChainRulesCore.rrule(
     function svd_trunc!_full_pullback(ΔUSV′)
         ΔUSV = unthunk.(ΔUSV′)
         Δt = svd_pullback!(
-            zeros(scalartype(t), space(t)), t, (U, S, V⁺), ΔUSV, inds;
+            zeros(storagetype(t), space(t)), t, (U, S, V⁺), ΔUSV, inds;
             gauge_atol = gtol(ΔUSV), degeneracy_atol = alg.rrule_alg.degeneracy_atol,
         )
         return NoTangent(), Δt, NoTangent()
@@ -358,7 +359,7 @@ function ChainRulesCore.rrule(
     function svd_trunc!_trunc_pullback(ΔUSV′)
         ΔUSV = unthunk.(ΔUSV′)
         Δf = svd_trunc_pullback!(
-            zeros(scalartype(t), space(t)), t, (U, S, V⁺), ΔUSV;
+            zeros(storagetype(t), space(t)), t, (U, S, V⁺), ΔUSV;
             gauge_atol = gtol(ΔUSV), degeneracy_atol = alg.rrule_alg.degeneracy_atol,
         )
         return NoTangent(), Δf, NoTangent()
@@ -383,7 +384,7 @@ function ChainRulesCore.rrule(
     function svd_trunc!_trunc_pullback(ΔUSV′)
         ΔUSV = unthunk.(ΔUSV′)
         Δf = svd_trunc_pullback!(
-            zeros(scalartype(t), space(t)), t, (U, S, V⁺), ΔUSV;
+            zeros(storagetype(t), space(t)), t, (U, S, V⁺), ΔUSV;
             gauge_atol = gtol(ΔUSV), degeneracy_atol = alg.rrule_alg.degeneracy_atol,
         )
         return NoTangent(), Δf, NoTangent()
@@ -395,6 +396,20 @@ function ChainRulesCore.rrule(
     return (U, S, V⁺), svd_trunc!_trunc_pullback
 end
 
+# `block(::AdjointTensorMap, c)` hands back a `LinearAlgebra.Adjoint`, and GPU array
+# packages only recognize a single layer of array wrappers (see the `WrappedArray` union in
+# Adapt). Slicing a doubly-wrapped array -- as `eachcol(V')` does -- therefore escapes the
+# device fast paths and falls back to scalar indexing on the host, so materialize the
+# wrapper before slicing.
+_materialize_block(m::AbstractMatrix) = m
+_materialize_block(m::Union{Adjoint, Transpose}) = copyto!(similar(m, size(m)), m)
+
+# Columns of `m` as a vector of `A`-typed vectors, keeping the data on its original device.
+function _column_vectors(::Type{A}, m::AbstractMatrix) where {A}
+    m̃ = _materialize_block(m)
+    return A[m̃[:, j] for j in axes(m̃, 2)]
+end
+
 # KrylovKit rrule compatible with TensorMaps & function handles
 function ChainRulesCore.rrule(
         ::typeof(svd_trunc!),
@@ -416,12 +431,15 @@ function ChainRulesCore.rrule(
         for (c, b) in blocks(Δf)
             Uc, Sc, Vc = block(U, c), block(S, c), block(V, c)
             ΔUc, ΔSc, ΔVc = block(ΔU, c), block(ΔS, c), block(ΔV, c)
-            Sdc = view(Sc, diagind(Sc))
-            ΔSdc = ΔSc isa AbstractZero ? ΔSc : view(ΔSc, diagind(ΔSc))
+            # `compute_svdsolve_pullback_data` reads the singular values element-wise and
+            # combines them with dense `n_vals × n_vals` host matrices, so keep these on
+            # the CPU; the bulk data (`lvecs`, `rvecs`, `block(f, c)`) stays on device.
+            Sdc = collect(diagview(Sc))
+            ΔSdc = ΔSc isa AbstractZero ? zero(Sdc) : collect(diagview(ΔSc))
 
             n_vals = length(Sdc)
-            lvecs = Vector{Vector{scalartype(f)}}(eachcol(Uc))
-            rvecs = Vector{Vector{scalartype(f)}}(eachcol(Vc'))
+            lvecs = _column_vectors(storagetype(f), Uc)
+            rvecs = _column_vectors(storagetype(f), Vc')
 
             # Dummy objects only used for warnings
             minimal_info = KrylovKit.ConvergenceInfo(n_vals, nothing, nothing, -1, -1)  # Only num. converged is used
@@ -431,12 +449,12 @@ function ChainRulesCore.rrule(
                 Δlvecs = fill(ZeroTangent(), n_vals)
                 Δrvecs = fill(ZeroTangent(), n_vals)
             else
-                Δlvecs = Vector{Vector{scalartype(f)}}(eachcol(ΔUc))
-                Δrvecs = Vector{Vector{scalartype(f)}}(eachcol(ΔVc'))
+                Δlvecs = _column_vectors(storagetype(f), ΔUc)
+                Δrvecs = _column_vectors(storagetype(f), ΔVc')
             end
 
             xs, ys = KrylovKitCRCExt.compute_svdsolve_pullback_data(
-                ΔSc isa AbstractZero ? fill(zero(Sc[1]), n_vals) : ΔSdc,
+                ΔSdc,
                 Δlvecs,
                 Δrvecs,
                 Sdc,
@@ -448,10 +466,14 @@ function ChainRulesCore.rrule(
                 minimal_alg,
                 rrule_alg,
             )
+            # `construct∂f_svd` hands back an `InplaceableThunk`; copying from it directly
+            # would fall back to iterating it element-wise, so unthunk it first.
             copyto!(
                 b,
-                KrylovKitCRCExt.construct∂f_svd(
-                    HasReverseMode(), block(f, c), lvecs, rvecs, xs, ys
+                unthunk(
+                    KrylovKitCRCExt.construct∂f_svd(
+                        HasReverseMode(), block(f, c), lvecs, rvecs, xs, ys
+                    )
                 ),
             )
         end
@@ -485,12 +507,15 @@ function ChainRulesCore.rrule(
         for (c, b) in blocks(Δf)
             Uc, Sc, Vc = block(U, c), block(S, c), block(V, c)
             ΔUc, ΔSc, ΔVc = block(ΔU, c), block(ΔS, c), block(ΔV, c)
-            Sdc = view(Sc, diagind(Sc))
-            ΔSdc = ΔSc isa AbstractZero ? ΔSc : view(ΔSc, diagind(ΔSc))
+            # `compute_svdsolve_pullback_data` reads the singular values element-wise and
+            # combines them with dense `n_vals × n_vals` host matrices, so keep these on
+            # the CPU; the bulk data (`lvecs`, `rvecs`, `block(f, c)`) stays on device.
+            Sdc = collect(diagview(Sc))
+            ΔSdc = ΔSc isa AbstractZero ? zero(Sdc) : collect(diagview(ΔSc))
 
             n_vals = length(Sdc)
-            lvecs = Vector{Vector{scalartype(f)}}(eachcol(Uc))
-            rvecs = Vector{Vector{scalartype(f)}}(eachcol(Vc'))
+            lvecs = _column_vectors(storagetype(f), Uc)
+            rvecs = _column_vectors(storagetype(f), Vc')
 
             # Dummy objects only used for warnings
             minimal_info = KrylovKit.ConvergenceInfo(n_vals, nothing, nothing, -1, -1)  # Only num. converged is used
@@ -500,12 +525,12 @@ function ChainRulesCore.rrule(
                 Δlvecs = fill(ZeroTangent(), n_vals)
                 Δrvecs = fill(ZeroTangent(), n_vals)
             else
-                Δlvecs = Vector{Vector{scalartype(f)}}(eachcol(ΔUc))
-                Δrvecs = Vector{Vector{scalartype(f)}}(eachcol(ΔVc'))
+                Δlvecs = _column_vectors(storagetype(f), ΔUc)
+                Δrvecs = _column_vectors(storagetype(f), ΔVc')
             end
 
             xs, ys = KrylovKitCRCExt.compute_svdsolve_pullback_data(
-                ΔSc isa AbstractZero ? fill(zero(Sc[1]), n_vals) : ΔSdc,
+                ΔSdc,
                 Δlvecs,
                 Δrvecs,
                 Sdc,
@@ -517,10 +542,14 @@ function ChainRulesCore.rrule(
                 minimal_alg,
                 rrule_alg,
             )
+            # `construct∂f_svd` hands back an `InplaceableThunk`; copying from it directly
+            # would fall back to iterating it element-wise, so unthunk it first.
             copyto!(
                 b,
-                KrylovKitCRCExt.construct∂f_svd(
-                    HasReverseMode(), block(f, c), lvecs, rvecs, xs, ys
+                unthunk(
+                    KrylovKitCRCExt.construct∂f_svd(
+                        HasReverseMode(), block(f, c), lvecs, rvecs, xs, ys
+                    )
                 ),
             )
         end
diff --git a/src/utility/symmetrization.jl b/src/utility/symmetrization.jl
index 1d98defc3..d5694ce21 100644
--- a/src/utility/symmetrization.jl
+++ b/src/utility/symmetrization.jl
@@ -52,7 +52,7 @@ function _fit_spaces(
     ) where {T, S <: IndexSpace, N₁, N₂}
     for i in 1:(N₁ + N₂)
         if space(x, i) ≠ space(y, i)
-            f = unitary(space(x, i) ← space(y, i))
+            f = unitary(TensorKit.promote_storagetype(y, x), space(x, i) ← space(y, i))
             y = permute(
                 ncon([f, y], [[-i, 1], [-(1:(i - 1))..., 1, -((i + 1):(N₁ + N₂))...]]),
                 (Tuple(1:N₁), Tuple((N₁ + 1):(N₁ + N₂))),
diff --git a/src/utility/util.jl b/src/utility/util.jl
index 603aec675..f2fc3620b 100644
--- a/src/utility/util.jl
+++ b/src/utility/util.jl
@@ -8,6 +8,24 @@ function _elementwise_mult(a₁::AbstractTensorMap, a₂::AbstractTensorMap)
 end
 
 _safe_pow(a::Number, pow::Real, tol::Real) = (pow < 0 && abs(a) < tol) ? zero(a) : a^pow
+# Same cutoff, but with the relative tolerance and the scale kept as separate arguments so
+# that the scale doesn't need to be copied back to the CPU memory.
+# `mx` is either a plain number or a 0-dimensional array, which broadcasts as a scalar.
+_safe_pow(a::Number, pow::Real, tol::Real, mx::Number) = _safe_pow(a, pow, tol * mx)
+
+"""
+    _maxabs(data)
+
+Largest absolute value in `data`, equal to `norm(_, Inf)` for the diagonal storage of a
+`DiagonalTensorMap`.
+
+Returns a scalar by default. GPU backends override this to return a 0-dimensional
+array, to avoid a copy back to the CPU memory.
+"""
+function _maxabs(data::AbstractArray)
+    isempty(data) && return zero(real(eltype(data)))
+    return LinearAlgebra.normInf(data)
+end
 
 """
     sdiag_pow(s, pow::Real; tol::Real=eps(real(scalartype(s)))^(3 / 4))
@@ -15,10 +33,7 @@ _safe_pow(a::Number, pow::Real, tol::Real) = (pow < 0 && abs(a) < tol) ? zero(a)
 Compute `s^pow` for a diagonal matrix `s`.
 """
 function sdiag_pow(s::DiagonalTensorMap, pow::Real; tol::Real = eps(real(scalartype(s)))^(3 / 4))
-    # Relative tol w.r.t. largest abs value of `s` (use norm(∘, Inf) to make differentiable)
-    tol *= norm(s, Inf)
-    spow = DiagonalTensorMap(_safe_pow.(s.data, pow, tol), space(s, 1))
-    return spow
+    return DiagonalTensorMap(_safe_pow.(s.data, pow, tol, _maxabs(s.data)), space(s, 1))
 end
 function sdiag_pow(
         s::AbstractTensorMap{T, S, 1, 1}, pow::Real; tol::Real = eps(real(scalartype(s)))^(3 / 4)
@@ -62,14 +77,23 @@ function absorb_s(U::AbstractTensorMap, S::DiagonalTensorMap, V::AbstractTensorM
     return U * sqrt_S, sqrt_S * V
 end
 
-_fliptwist_s(s::DiagonalTensorMap) = twist!(DiagonalTensorMap(flip(s, 1:2)), 1)
+# returns new but destroys old S
+function _fliptwist_s!(s::DiagonalTensorMap)
+    for (f₁, f₂) in fusiontrees(s)
+        data = s[f₁, f₂]
+        _, θ = only(flip((f₁, f₂), (1, 2)))
+        θ *= twist(f₁.uncoupled[1])
+        scale!(data, θ)
+    end
+    return DiagonalTensorMap(s.data, flip(s.domain))
+end
 
 # Check whether diagonals contain degenerate values up to absolute or relative tolerance
 function is_degenerate_spectrum(
         S; atol::Real = 0, rtol::Real = atol > 0 ? 0 : sqrt(eps(scalartype(S)))
     )
     for (_, b) in blocks(S)
-        s = real(diag(b))
+        s = real(collect(diag(b)))
         for i in 1:(length(s) - 1)
             isapprox(s[i], s[i + 1]; atol, rtol) && return true
         end
diff --git a/test/Project.toml b/test/Project.toml
index 86bdd9467..d4d8e0ec4 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -3,13 +3,19 @@ name = "PEPSKitTests"
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+CUPTI = "9e67e8f6-ba02-4b6c-a7db-3b11ae1e7ab7"
+CUDATools = "9ec180c6-1c07-47c7-9e6e-ebefa4d1f6d0"
+GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
 MPSKitModels = "ca635005-6f8c-4cd1-b51d-8491250ef2ab"
 MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
+NVML = "611af6d1-644e-4c5d-bd58-854d7d1254b9"
 OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
 PEPSKit = "52969e89-939e-4361-9b68-9bc7cde4bdeb"
 ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
@@ -21,13 +27,23 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+cuBLAS = "182d3088-87b7-4494-8cad-fc6afaa545bc"
+cuFFT = "533571aa-0936-420e-b4be-9c66f5f626ca"
+cuRAND = "20fd9a0b-12d5-4c2f-a8af-7c34e9e60431"
+cuSOLVER = "887afef0-6a32-4de5-add4-7827692ba8fc"
+cuSPARSE = "b26da814-b3bc-49ef-b0ee-c816305aa060"
 
 [sources]
 PEPSKit = {path = ".."}
+GPUArrays = {url = "https://github.com/JuliaGPU/GPUArrays.jl", rev = "main"}
+AMDGPU = {url = "https://github.com/JuliaGPU/AMDGPU.jl", rev = "main"}
 
 [compat]
 Adapt = "4"
+AMDGPU = "2"
 ChainRulesTestUtils = "1.13"
+CUDA = "6"
+CUDATools = "6"
 ParallelTestRunner = "2.6.0"
 QuadGK = "2.11.1"
 Test = "1"
diff --git a/test/bondenv/benv_ctm.jl b/test/bondenv/benv_ctm.jl
index b740fc9d2..08f559f8f 100644
--- a/test/bondenv/benv_ctm.jl
+++ b/test/bondenv/benv_ctm.jl
@@ -1,5 +1,4 @@
-using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -9,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.bondenv_ctm(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.bondenv_ctm(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.bondenv_ctm(ROCArray)
+end
diff --git a/test/bondenv/benv_gaugefix.jl b/test/bondenv/benv_gaugefix.jl
index 4d25ba037..8df16a076 100644
--- a/test/bondenv/benv_gaugefix.jl
+++ b/test/bondenv/benv_gaugefix.jl
@@ -1,5 +1,4 @@
-using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -9,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.bondenv_gaugefix(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.bondenv_gaugefix(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.bondenv_gaugefix(ROCArray)
+end
diff --git a/test/bondenv/bond_truncate.jl b/test/bondenv/bond_truncate.jl
index eb14e97cb..662dcbf27 100644
--- a/test/bondenv/bond_truncate.jl
+++ b/test/bondenv/bond_truncate.jl
@@ -1,5 +1,4 @@
-using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -9,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.bondenv_truncate(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.bondenv_truncate(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.bondenv_truncate(ROCArray)
+end
diff --git a/test/boundarymps/vumps.jl b/test/boundarymps/vumps.jl
index 1906704fb..4db249123 100644
--- a/test/boundarymps/vumps.jl
+++ b/test/boundarymps/vumps.jl
@@ -1,4 +1,4 @@
-Random.seed!(29384293742893)
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -11,3 +11,17 @@ if !is_buildkite
     TestSuite.boundary_mps_fermionic_peps(Vector)
     TestSuite.boundary_mps_pepo_runthrough(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.boundary_mps_one_one_peps(CuArray)
+    TestSuite.boundary_mps_two_two_peps(CuArray)
+    TestSuite.boundary_mps_fermionic_peps(CuArray)
+    TestSuite.boundary_mps_pepo_runthrough(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.boundary_mps_one_one_peps(ROCArray)
+    TestSuite.boundary_mps_two_two_peps(ROCArray)
+    TestSuite.boundary_mps_fermionic_peps(ROCArray)
+    TestSuite.boundary_mps_pepo_runthrough(ROCArray)
+end
diff --git a/test/bp/expvals.jl b/test/bp/expvals.jl
index 73a507472..13a5c7379 100644
--- a/test/bp/expvals.jl
+++ b/test/bp/expvals.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -8,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.bp_expvals(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.bp_expvals(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.bp_expvals(ROCArray)
+end
diff --git a/test/bp/gaugefix.jl b/test/bp/gaugefix.jl
index bff0d4cac..251792d15 100644
--- a/test/bp/gaugefix.jl
+++ b/test/bp/gaugefix.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -8,3 +8,13 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.bp_gaugefix_bp_vs_su(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.bp_gaugefix_bp_vs_su(CuArray)
+end
+
+if AMDGPU.functional()
+    # rocSOLVER doesn't offer a general eigensolver yet,
+    # only Hermitian
+    TestSuite.bp_gaugefix_bp_vs_su(ROCArray; posdef_msgs = [true])
+end
diff --git a/test/bp/rotation.jl b/test/bp/rotation.jl
index 984594105..3777b7ab0 100644
--- a/test/bp/rotation.jl
+++ b/test/bp/rotation.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -8,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.bp_rotations(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.bp_rotations(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.bp_rotations(ROCArray)
+end
diff --git a/test/bp/unitcell.jl b/test/bp/unitcell.jl
index 7e7e31a28..9276c9fd6 100644
--- a/test/bp/unitcell.jl
+++ b/test/bp/unitcell.jl
@@ -1,5 +1,4 @@
-using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -10,3 +9,13 @@ if !is_buildkite
     TestSuite.bp_unitcell_random_cartesian_spaces(Vector)
     TestSuite.bp_unitcell_specific_u1_spaces(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.bp_unitcell_random_cartesian_spaces(CuArray)
+    TestSuite.bp_unitcell_specific_u1_spaces(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.bp_unitcell_random_cartesian_spaces(ROCArray)
+    TestSuite.bp_unitcell_specific_u1_spaces(ROCArray)
+end
diff --git a/test/compress/local.jl b/test/compress/local.jl
index df0be6fcf..52fd94b6c 100644
--- a/test/compress/local.jl
+++ b/test/compress/local.jl
@@ -1,5 +1,4 @@
-using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -11,3 +10,15 @@ if !is_buildkite
     TestSuite.compress_cost_function(Vector)
     TestSuite.compress_virtual_space_matching(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.compress_fermionic_twists(CuArray)
+    TestSuite.compress_cost_function(CuArray)
+    TestSuite.compress_virtual_space_matching(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.compress_fermionic_twists(ROCArray)
+    TestSuite.compress_cost_function(ROCArray)
+    TestSuite.compress_virtual_space_matching(ROCArray)
+end
diff --git a/test/ctmrg/contractions.jl b/test/ctmrg/contractions.jl
index c47ef3865..8180e2d00 100644
--- a/test/ctmrg/contractions.jl
+++ b/test/ctmrg/contractions.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -10,3 +10,15 @@ if !is_buildkite
     TestSuite.ctmrg_contractions_specific_u1_spaces(Vector)
     TestSuite.ctmrg_contractions_random_fermionic_spaces(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.ctmrg_contractions_random_cartesian_spaces(CuArray)
+    TestSuite.ctmrg_contractions_specific_u1_spaces(CuArray)
+    TestSuite.ctmrg_contractions_random_fermionic_spaces(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.ctmrg_contractions_random_cartesian_spaces(ROCArray)
+    TestSuite.ctmrg_contractions_specific_u1_spaces(ROCArray)
+    TestSuite.ctmrg_contractions_random_fermionic_spaces(ROCArray)
+end
diff --git a/test/ctmrg/correlation_length.jl b/test/ctmrg/correlation_length.jl
index e3998e7b8..914b817b6 100644
--- a/test/ctmrg/correlation_length.jl
+++ b/test/ctmrg/correlation_length.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -8,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.ctmrg_correlation_length(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.ctmrg_correlation_length(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.ctmrg_correlation_length(ROCArray)
+end
diff --git a/test/ctmrg/fixed_iterscheme.jl b/test/ctmrg/fixed_iterscheme.jl
index 50edc6ab8..81f6a969c 100644
--- a/test/ctmrg/fixed_iterscheme.jl
+++ b/test/ctmrg/fixed_iterscheme.jl
@@ -1,5 +1,4 @@
-using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -11,3 +10,17 @@ if !is_buildkite
     TestSuite.ctmrg_fixed_iterscheme_c4v(Vector)
     TestSuite.ctmrg_fixed_iterscheme_divide_and_conquer(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.ctmrg_fixed_iterscheme_asymmetric(CuArray; svd_alg = :QRIteration)
+    # CUSOLVER doesn't provide heev
+    TestSuite.ctmrg_fixed_iterscheme_c4v(CuArray; eigh_alg = :DivideAndConquer)
+    # CUSOLVER doesn't provide gesdd
+    #TestSuite.ctmrg_fixed_iterscheme_divide_and_conquer(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.ctmrg_fixed_iterscheme_asymmetric(ROCArray)
+    TestSuite.ctmrg_fixed_iterscheme_c4v(ROCArray)
+    TestSuite.ctmrg_fixed_iterscheme_divide_and_conquer(ROCArray)
+end
diff --git a/test/ctmrg/flavors.jl b/test/ctmrg/flavors.jl
index 557c9d451..93811981c 100644
--- a/test/ctmrg/flavors.jl
+++ b/test/ctmrg/flavors.jl
@@ -1,3 +1,5 @@
+using PEPSKit, CUDA, AMDGPU
+
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
 
@@ -8,3 +10,16 @@ if !is_buildkite
     TestSuite.ctmrg_flavors_fixedspace_truncation(Vector)
     TestSuite.ctmrg_flavors_c4v(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.ctmrg_flavors_unitcells(CuArray; minimal = true)
+    TestSuite.ctmrg_flavors_fixedspace_truncation(CuArray; minimal = true)
+    # CUSOLVER doesn't provide heev
+    TestSuite.ctmrg_flavors_c4v(CuArray; eigh_alg = :DivideAndConquer, minimal = true)
+end
+
+if AMDGPU.functional()
+    TestSuite.ctmrg_flavors_unitcells(ROCArray; minimal = true)
+    TestSuite.ctmrg_flavors_fixedspace_truncation(ROCArray; minimal = true)
+    TestSuite.ctmrg_flavors_c4v(ROCArray; minimal = true)
+end
diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl
index a8e6e1afd..94c11885b 100644
--- a/test/ctmrg/gaugefix.jl
+++ b/test/ctmrg/gaugefix.jl
@@ -1,5 +1,4 @@
-using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -10,3 +9,13 @@ if !is_buildkite
     TestSuite.ctmrg_gaugefix_asymmetric(Vector)
     TestSuite.ctmrg_gaugefix_c4v(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.ctmrg_gaugefix_asymmetric(CuArray; minimal = true)
+    TestSuite.ctmrg_gaugefix_c4v(CuArray; minimal = true)
+end
+
+if AMDGPU.functional()
+    TestSuite.ctmrg_gaugefix_asymmetric(ROCArray; minimal = true)
+    TestSuite.ctmrg_gaugefix_c4v(ROCArray; minimal = true)
+end
diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index 4432e54da..d8444ce90 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -9,3 +9,13 @@ if !is_buildkite
     TestSuite.ctmrg_initialization_critical_ising(Vector)
     TestSuite.ctmrg_initialization_peps(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.ctmrg_initialization_critical_ising(CuArray)
+    TestSuite.ctmrg_initialization_peps(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.ctmrg_initialization_critical_ising(ROCArray)
+    TestSuite.ctmrg_initialization_peps(ROCArray)
+end
diff --git a/test/ctmrg/jacobian_real_linear.jl b/test/ctmrg/jacobian_real_linear.jl
index b9d66acf5..3298d651a 100644
--- a/test/ctmrg/jacobian_real_linear.jl
+++ b/test/ctmrg/jacobian_real_linear.jl
@@ -1,5 +1,4 @@
-using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -9,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.ctmrg_jacobian_real_linear(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.ctmrg_jacobian_real_linear(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.ctmrg_jacobian_real_linear(ROCArray)
+end
diff --git a/test/ctmrg/partition_function.jl b/test/ctmrg/partition_function.jl
index 50188bd3b..4b8ea309f 100644
--- a/test/ctmrg/partition_function.jl
+++ b/test/ctmrg/partition_function.jl
@@ -1,5 +1,4 @@
-using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -10,3 +9,13 @@ if !is_buildkite
     TestSuite.ctmrg_partition_function_spaces(Vector)
     TestSuite.ctmrg_partition_function(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.ctmrg_partition_function_spaces(CuArray)
+    TestSuite.ctmrg_partition_function(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.ctmrg_partition_function_spaces(ROCArray)
+    TestSuite.ctmrg_partition_function(ROCArray)
+end
diff --git a/test/ctmrg/pepo.jl b/test/ctmrg/pepo.jl
index be56f52c7..6fd7e6c5e 100644
--- a/test/ctmrg/pepo.jl
+++ b/test/ctmrg/pepo.jl
@@ -1,5 +1,4 @@
-using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -10,3 +9,13 @@ if !is_buildkite
     TestSuite.ctmrg_pepo_runthroughs(Vector)
     TestSuite.ctmrg_pepo_fixed_point(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.ctmrg_pepo_runthroughs(CuArray; minimal = true)
+    TestSuite.ctmrg_pepo_fixed_point(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.ctmrg_pepo_runthroughs(ROCArray; minimal = true)
+    TestSuite.ctmrg_pepo_fixed_point(ROCArray)
+end
diff --git a/test/ctmrg/suweight.jl b/test/ctmrg/suweight.jl
index fa1cc1984..67f0eedef 100644
--- a/test/ctmrg/suweight.jl
+++ b/test/ctmrg/suweight.jl
@@ -1,5 +1,4 @@
-using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -9,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.ctmrg_suweight(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.ctmrg_suweight(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.ctmrg_suweight(ROCArray)
+end
diff --git a/test/ctmrg/unitcell.jl b/test/ctmrg/unitcell.jl
index 6599ce32f..53c80b787 100644
--- a/test/ctmrg/unitcell.jl
+++ b/test/ctmrg/unitcell.jl
@@ -1,5 +1,5 @@
 using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -10,3 +10,13 @@ if !is_buildkite
     TestSuite.ctmrg_unitcell_random_cartesian_spaces(Vector)
     TestSuite.ctmrg_unitcell_specific_u1_spaces(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.ctmrg_unitcell_random_cartesian_spaces(CuArray; minimal = true)
+    TestSuite.ctmrg_unitcell_specific_u1_spaces(CuArray; minimal = true)
+end
+
+if AMDGPU.functional()
+    TestSuite.ctmrg_unitcell_random_cartesian_spaces(ROCArray; minimal = true)
+    TestSuite.ctmrg_unitcell_specific_u1_spaces(ROCArray; minimal = true)
+end
diff --git a/test/gradients/c4v_ctmrg_gradients.jl b/test/gradients/c4v_ctmrg_gradients.jl
index 9bc6f9905..bbb900dcf 100644
--- a/test/gradients/c4v_ctmrg_gradients.jl
+++ b/test/gradients/c4v_ctmrg_gradients.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -8,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.gradients_c4v(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.gradients_c4v(CuArray; minimal = true)
+end
+
+if AMDGPU.functional()
+    TestSuite.gradients_c4v(ROCArray; minimal = true)
+end
diff --git a/test/gradients/ctmrg_gradients.jl b/test/gradients/ctmrg_gradients.jl
index ce50580a8..7744a93c9 100644
--- a/test/gradients/ctmrg_gradients.jl
+++ b/test/gradients/ctmrg_gradients.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -9,3 +9,13 @@ if !is_buildkite
     TestSuite.gradients_asymmetric(Vector)
     TestSuite.gradients_asymmetric_276(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.gradients_asymmetric(CuArray; minimal = true)
+    TestSuite.gradients_asymmetric_276(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.gradients_asymmetric(ROCArray; minimal = true)
+    TestSuite.gradients_asymmetric_276(ROCArray)
+end
diff --git a/test/testsuite/bondenv/benv_ctm.jl b/test/testsuite/bondenv/benv_ctm.jl
index b740d58a0..8332d276b 100644
--- a/test/testsuite/bondenv/benv_ctm.jl
+++ b/test/testsuite/bondenv/benv_ctm.jl
@@ -13,7 +13,7 @@ trunc_state = truncerror(; atol = 1.0e-10) & truncrank(4)
 ctm_alg = SequentialCTMRG(; tol = 1.0e-10, verbosity = 2, trunc = truncerror(; atol = 1.0e-10) & truncrank(8))
 # create Hubbard iPEPS using simple update
 function get_hubbard_peps(AT, t::Float64 = 1.0, U::Float64 = 8.0)
-    H = hubbard_model(ComplexF64, Trivial, U1Irrep, InfiniteSquare(Nr, Nc); t, U, mu = U / 2)
+    H = adapt(AT, hubbard_model(ComplexF64, Trivial, U1Irrep, InfiniteSquare(Nr, Nc); t, U, mu = U / 2))
     Vphy = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1 // 2) => 1, (1, -1 // 2) => 1)
     peps = adapt(AT, InfinitePEPS(rand, ComplexF64, Vphy, Vphy; unitcell = (Nr, Nc)))
     wts = SUWeight(peps)
diff --git a/test/testsuite/boundarymps/vumps.jl b/test/testsuite/boundarymps/vumps.jl
index df07e1593..9f8545750 100644
--- a/test/testsuite/boundarymps/vumps.jl
+++ b/test/testsuite/boundarymps/vumps.jl
@@ -15,14 +15,21 @@ function boundary_mps_one_one_peps(AT)
     return @testset "(1, 1) PEPS ($AT)" begin
         Vpeps = ComplexSpace(2)
         psi = adapt(AT, InfinitePEPS(Vpeps, Vpeps))
+        @test storagetype(psi) <: AT
         T = PEPSKit.InfiniteTransferPEPS(psi, 1, 1)
+        @test storagetype(T) <: AT
         foreach(V -> (@test V == Vpeps ⊗ Vpeps'), physicalspace(T))
         mps = initialize_mps(T, [ComplexSpace(20)])
+        @test storagetype(mps) <: AT
 
         mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
         N = abs(sum(expectation_value(mps, T)))
 
-        mps2, = changebonds(mps, T, OptimalExpand(; trscheme = truncrank(30))) # TODO: update `trscheme` to `trunc` once MPSKit does
+        @static if VERSION < v"1.11.0-rc" # so [sources] isn't used
+            mps2, = changebonds(mps, T, OptimalExpand(; trscheme = truncrank(30))) # TODO: update `trscheme` to `trunc` once MPSKit does
+        else
+            mps2, = changebonds(mps, T, OptimalExpand(; trunc = truncrank(30))) # TODO: update `trunc` to `trunc` once MPSKit does
+        end
         mps2, env2, ϵ = leading_boundary(mps2, T, vumps_alg)
         N2 = abs(sum(expectation_value(mps2, T)))
         @test N ≈ N2 rtol = 1.0e-2
@@ -39,9 +46,12 @@ function boundary_mps_two_two_peps(AT)
     return @testset "(2, 2) PEPS ($AT)" begin
         Vpeps = ComplexSpace(2)
         psi = adapt(AT, InfinitePEPS(Vpeps, Vpeps; unitcell = (2, 2)))
+        @test storagetype(psi) <: AT
         T = PEPSKit.MultilineTransferPEPS(psi, 1)
+        @test storagetype(T) <: AT
         # foreach(V -> (@test V == Vpeps ⊗ Vpeps'), physicalspace(T)) # TODO: MPSKit.physicalspace(::MultilineMPO) isn't implemented...
-        mps = initialize_mps(rand, scalartype(T), T, fill(ComplexSpace(20), 2, 2))
+        mps = adapt(AT, initialize_mps(rand, scalartype(T), T, fill(ComplexSpace(20), 2, 2)))
+        @test storagetype(mps) <: AT
         mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
         N = abs(prod(expectation_value(mps, T)))
 
@@ -60,12 +70,16 @@ function boundary_mps_fermionic_peps(AT)
         χ = Vect[fℤ₂](0 => 10, 1 => 10)
 
         psi = adapt(AT, InfinitePEPS(D, d; unitcell = (1, 1)))
+        @test storagetype(psi) <: AT
         n = InfiniteSquareNetwork(psi)
+        @test storagetype(n) <: AT
         T = InfiniteTransferPEPS(psi, 1, 1)
+        @test storagetype(T) <: AT
         foreach(V -> (@test V == D ⊗ D'), physicalspace(T))
 
         # compare boundary MPS contraction to CTMRG contraction
-        mps = initialize_mps(T, [χ])
+        mps = adapt(AT, initialize_mps(T, [χ]))
+        @test storagetype(mps) <: AT
         mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
         N_vumps = abs(prod(expectation_value(mps, T)))
 
@@ -116,20 +130,25 @@ function boundary_mps_pepo_runthrough(AT)
         # single-layer PEPO
         O = ising_pepo(1)
         psi = adapt(AT, PEPSKit.initializePEPS(O, Vpeps))
-        T = InfiniteTransferPEPO(psi, O, 1, 1)
+        @test storagetype(psi) <: AT
+        T = InfiniteTransferPEPO(psi, adapt(AT, O), 1, 1)
+        @test storagetype(T) <: AT
         foreach(V -> (@test V == Vpeps ⊗ Vpepo ⊗ Vpeps'), physicalspace(T))
 
-        mps = initialize_mps(rand, scalartype(T), T, [ComplexSpace(10)])
+        mps = adapt(AT, initialize_mps(rand, scalartype(T), T, [ComplexSpace(10)]))
+        @test storagetype(mps) <: AT
         mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
         f = abs(prod(expectation_value(mps, T)))
 
         # double-layer PEPO
         O2 = repeat(O, 1, 1, 2)
-        psi2 = initializePEPS(O2, Vpeps)
+        psi2 = adapt(AT, initializePEPS(O2, Vpeps))
+        @test storagetype(psi2) <: AT
         T2 = InfiniteTransferPEPO(psi, O2, 1, 1)
         foreach(V -> (@test V == Vpeps ⊗ Vpepo ⊗ Vpepo ⊗ Vpeps'), physicalspace(T2))
 
-        mps2 = initialize_mps(rand, scalartype(T2), T2, [ComplexSpace(8)])
+        mps2 = adapt(AT, initialize_mps(rand, scalartype(T2), T2, [ComplexSpace(8)]))
+        @test storagetype(mps2) <: AT
         mps2, env2, ϵ = leading_boundary(mps2, T2, vumps_alg)
         f = abs(prod(expectation_value(mps2, T2)))
     end
diff --git a/test/testsuite/bp/gaugefix.jl b/test/testsuite/bp/gaugefix.jl
index 404784e1b..c0a4df34f 100644
--- a/test/testsuite/bp/gaugefix.jl
+++ b/test/testsuite/bp/gaugefix.jl
@@ -6,10 +6,10 @@ using PEPSKit, Adapt
 using PEPSKit: compare_weights, random_dual!, twistdual
 using PEPSKit: _next, _is_bipartite
 
-function bp_gaugefix_bp_vs_su(AT)
+function bp_gaugefix_bp_vs_su(AT; posdef_msgs = [true, false])
     return @testset "BP vs SU ($AT) ($S, bipartite = $(bipartite), posdef msgs = $h)" for
         (S, bipartite, h) in Iterators.product(
-            [U1Irrep, FermionParity], [true, false], [true, false]
+            [U1Irrep, FermionParity], [true, false], posdef_msgs
         )
         unitcell = bipartite ? (2, 2) : (2, 3)
         elt = ComplexF64
@@ -51,9 +51,11 @@ function bp_gaugefix_bp_vs_su(AT)
                 peps0[2, c] = copy(peps0[1, c + 1])
             end
         end
+        @test storagetype(peps0) <: AT
 
         # start by gauging with SU
         peps1, wts1 = gauge_fix(peps0, SUGauge(; maxiter, tol))
+        @test storagetype(peps1) <: AT
         for (a0, a1) in zip(peps0.A, peps1.A)
             @test space(a0) == space(a1)
         end
@@ -66,6 +68,7 @@ function bp_gaugefix_bp_vs_su(AT)
         # find BP fixed point and SUWeight
         bp_alg = BeliefPropagation(; maxiter, tol, bipartite, project_hermitian = h)
         env = BPEnv(randn, elt, peps1; posdef = h)
+        @test storagetype(env) <: AT
         env, err = leading_boundary(env, peps1, bp_alg)
         if bipartite
             @test _is_bipartite(env)
@@ -85,7 +88,7 @@ function bp_gaugefix_bp_vs_su(AT)
         for (X, Xinv) in XXinv
             # X, Xinv should contract to identity
             @tensor tmp[-1; -2] := X[-1; 1] * Xinv[1; -2]
-            @test tmp ≈ twistdual(TensorKit.id(space(X, 1)), 1)
+            @test tmp ≈ twistdual(adapt(AT, TensorKit.id(space(X, 1))), 1)
             # BP should differ from SU only by a unitary gauge transformation
             @test inv(X) ≈ adjoint(X) ≈ Xinv
         end
diff --git a/test/testsuite/bp/rotation.jl b/test/testsuite/bp/rotation.jl
index e38f9131e..778077790 100644
--- a/test/testsuite/bp/rotation.jl
+++ b/test/testsuite/bp/rotation.jl
@@ -33,7 +33,9 @@ function bp_rotations(AT)
         ψDNs = random_dual!(fill(D, unitcell))
         ψDEs = random_dual!(fill(D, unitcell))
         ψ = adapt(AT, InfinitePEPS(ψds, ψDNs, ψDEs))
+        @test storagetype(ψ) <: AT
         env = BPEnv(ψ)
+        @test storagetype(env) <: AT
 
         op = adapt(AT, randn(d → d))
         meas1 = meas_sites(op, ψ, env)
diff --git a/test/testsuite/bp/unitcell.jl b/test/testsuite/bp/unitcell.jl
index 9c88d281b..6c178c074 100644
--- a/test/testsuite/bp/unitcell.jl
+++ b/test/testsuite/bp/unitcell.jl
@@ -9,7 +9,9 @@ elt = ComplexF64
 
 function test_unitcell(AT, unitcell, Pspaces, Nspaces, Espaces)
     peps = adapt(AT, InfinitePEPS(randn, elt, Pspaces, Nspaces, Espaces))
+    @test storagetype(peps) <: AT
     env0 = BPEnv(ones, elt, peps)
+    @test storagetype(env0) <: AT
     alg = BeliefPropagation()
 
     # apply one BP iteration
@@ -21,7 +23,7 @@ function test_unitcell(AT, unitcell, Pspaces, Nspaces, Espaces)
     # compute random expecation value to test matching bonds
     random_op = LocalOperator(
         Pspaces, (
-            (c,) => randn(elt, Pspaces[c], Pspaces[c])
+            (c,) => adapt(AT, randn(elt, Pspaces[c], Pspaces[c]))
                 for c in CartesianIndices(unitcell)
         )...,
     )
diff --git a/test/testsuite/ctmrg/fixed_iterscheme.jl b/test/testsuite/ctmrg/fixed_iterscheme.jl
index a360aa70a..f0f982283 100644
--- a/test/testsuite/ctmrg/fixed_iterscheme.jl
+++ b/test/testsuite/ctmrg/fixed_iterscheme.jl
@@ -20,12 +20,12 @@ using PEPSKit.Defaults: ctmrg_tol
 # initialize parameters
 D = 2
 χ = 16
-svd_algs = [(; alg = :DivideAndConquer), (; alg = :GKL)]
 projector_algs_asymm = [:HalfInfiniteProjector] #, :FullInfiniteProjector]
 unitcells = [(1, 1), (3, 4)]
 atol = 1.0e-5
 
-function ctmrg_fixed_iterscheme_asymmetric(AT)
+function ctmrg_fixed_iterscheme_asymmetric(AT; svd_alg = :DivideAndConquer)
+    svd_algs = [(; alg = svd_alg), (; alg = :GKL)]
     # test for element-wise convergence after application of fixed step
     return @testset "$unitcell unit cell with $(decomposition_alg.alg) and $projector_alg ($AT)" for (
             unitcell, decomposition_alg, projector_alg,
@@ -61,12 +61,12 @@ function ctmrg_fixed_iterscheme_asymmetric(AT)
 end
 
 # test same thing for C4v CTMRG
-c4v_algs = [
-    (:C4vQRProjector, (; alg = :Householder)),
-    (:C4vEighProjector, (; alg = :QRIteration)),
-    (:C4vEighProjector, (; alg = :Lanczos)),
-]
-function ctmrg_fixed_iterscheme_c4v(AT)
+function ctmrg_fixed_iterscheme_c4v(AT; eigh_alg = :QRIteration)
+    c4v_algs = [
+        (:C4vQRProjector, (; alg = :Householder)),
+        (:C4vEighProjector, (; alg = eigh_alg)),
+        (:C4vEighProjector, (; alg = :Lanczos)),
+    ]
     return @testset "$(decomposition_alg.alg) and $projector_alg ($AT)" for
         (projector_alg, decomposition_alg) in c4v_algs
         # initialize states
diff --git a/test/testsuite/ctmrg/flavors.jl b/test/testsuite/ctmrg/flavors.jl
index ef8833b56..a849124a9 100644
--- a/test/testsuite/ctmrg/flavors.jl
+++ b/test/testsuite/ctmrg/flavors.jl
@@ -11,15 +11,17 @@ D = 2
 χ = 16
 unitcells = [(1, 1), (3, 4)]
 projector_algs_asymm = [:HalfInfiniteProjector, :FullInfiniteProjector]
-projector_algs_c4v = [
-    (:C4vQRProjector, :Householder),
-    (:C4vEighProjector, :QRIteration), (:C4vEighProjector, :Lanczos),
-]
 Ts = [Float64, ComplexF64]
 
-function ctmrg_flavors_unitcells(AT)
+# minimal subsets of the combinations tested below which still cover every option value at least once
+minimal_combinations_unitcells = [((1, 1), :HalfInfiniteProjector), ((3, 4), :FullInfiniteProjector)]
+minimal_combinations_fixedspace = [
+    (:SequentialCTMRG, :FullInfiniteProjector), (:SimultaneousCTMRG, :HalfInfiniteProjector),
+]
+
+function ctmrg_flavors_unitcells(AT; minimal::Bool = false)
     return @testset "$(unitcell) unit cell with $projector_alg ($AT)" for (unitcell, projector_alg) in
-        Iterators.product(unitcells, projector_algs_asymm)
+        (minimal ? minimal_combinations_unitcells : Iterators.product(unitcells, projector_algs_asymm))
         # compute environments
         Random.seed!(32350283290358)
         psi = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(D); unitcell))
@@ -52,10 +54,12 @@ function ctmrg_flavors_unitcells(AT)
     end
 end
 
-function ctmrg_flavors_fixedspace_truncation(AT)
+function ctmrg_flavors_fixedspace_truncation(AT; minimal::Bool = false)
     # test fixedspace actually fixes space
-    return @testset "Fixedspace truncation using $alg and $projector_alg ($AT)" for (alg, projector_alg) in
-        Iterators.product([:SequentialCTMRG, :SimultaneousCTMRG], projector_algs_asymm)
+    return @testset "Fixedspace truncation using $alg and $projector_alg ($AT)" for (alg, projector_alg) in (
+            minimal ? minimal_combinations_fixedspace :
+                Iterators.product([:SequentialCTMRG, :SimultaneousCTMRG], projector_algs_asymm)
+        )
         Ds = ComplexSpace.(fill(2, 3, 3))
         χs = ComplexSpace.([16 17 18; 15 20 21; 14 19 22])
         psi = adapt(AT, InfinitePEPS(Ds, Ds, Ds))
@@ -70,9 +74,16 @@ function ctmrg_flavors_fixedspace_truncation(AT)
     end
 end
 
-function ctmrg_flavors_c4v(AT)
-    return @testset "C4v with ($T) - ($projector_alg, $decomp_alg) ($AT)" for (T, (projector_alg, decomp_alg)) in
+function ctmrg_flavors_c4v(AT; eigh_alg = :QRIteration, minimal::Bool = false)
+    projector_algs_c4v = [
+        (:C4vQRProjector, :Householder),
+        (:C4vEighProjector, eigh_alg), (:C4vEighProjector, :Lanczos),
+    ]
+    # minimal: test each projector alg once, alternating between scalar types
+    combinations = minimal ? zip(Iterators.cycle(Ts), projector_algs_c4v) :
         Iterators.product(Ts, projector_algs_c4v)
+    return @testset "C4v with ($T) - ($projector_alg, $decomp_alg) ($AT)" for (T, (projector_alg, decomp_alg)) in
+        combinations
 
         Random.seed!(29358293829382)
         symm = RotateReflect()
diff --git a/test/testsuite/ctmrg/gaugefix.jl b/test/testsuite/ctmrg/gaugefix.jl
index ce6597734..82f6ebe23 100644
--- a/test/testsuite/ctmrg/gaugefix.jl
+++ b/test/testsuite/ctmrg/gaugefix.jl
@@ -15,6 +15,17 @@ projector_algs_asymm = [:HalfInfiniteProjector, :FullInfiniteProjector]
 projector_algs_c4v = [:C4vEighProjector, :C4vQRProjector]
 gauge_algs_asymm = [ScramblingEnvGauge()]
 gauge_algs_c4v = [ScramblingEnvGaugeC4v()]
+
+# minimal subsets of the combinations above which still cover every option value at least once
+minimal_combinations_asymm = [
+    (ComplexSpace, Float64, (1, 1), SequentialCTMRG, :HalfInfiniteProjector, ScramblingEnvGauge()),
+    (Z2Space, ComplexF64, (2, 2), SimultaneousCTMRG, :FullInfiniteProjector, ScramblingEnvGauge()),
+    (ComplexSpace, ComplexF64, (3, 2), SimultaneousCTMRG, :HalfInfiniteProjector, ScramblingEnvGauge()),
+]
+minimal_combinations_c4v = [
+    (ComplexSpace, Float64, :C4vEighProjector, ScramblingEnvGaugeC4v()),
+    (Z2Space, ComplexF64, :C4vQRProjector, ScramblingEnvGaugeC4v()),
+]
 tol = 1.0e-6  # large tol due to χ=6
 χ = 6
 atol = 1.0e-4
@@ -72,11 +83,14 @@ function _preconverged_env_c4v(S, ::Type{T}) where {T}
     end
 end
 
-function ctmrg_gaugefix_asymmetric(AT)
+function ctmrg_gaugefix_asymmetric(AT; minimal::Bool = false)
     return @testset "($S) - ($T) - ($unitcell) - ($ctmrg_alg) - ($projector_alg) - ($gauge_alg) - ($AT)" for (
             S, T, unitcell, ctmrg_alg, projector_alg, gauge_alg,
-        ) in Iterators.product(
-            spacetypes, scalartypes, unitcells, ctmrg_algs_asymm, projector_algs_asymm, gauge_algs_asymm
+        ) in (
+            minimal ? minimal_combinations_asymm :
+                Iterators.product(
+                    spacetypes, scalartypes, unitcells, ctmrg_algs_asymm, projector_algs_asymm, gauge_algs_asymm
+                )
         )
         alg = ctmrg_alg(; tol, projector_alg)
         env_pre, psi = _preconverged_env(S, T, unitcell)
@@ -92,11 +106,12 @@ function ctmrg_gaugefix_asymmetric(AT)
     end
 end
 
-function ctmrg_gaugefix_c4v(AT)
+function ctmrg_gaugefix_c4v(AT; minimal::Bool = false)
     return @testset "($S) - ($T) - ($projector_alg) - ($gauge_alg) - ($AT)" for (
             S, T, projector_alg, gauge_alg,
-        ) in Iterators.product(
-            spacetypes, scalartypes, projector_algs_c4v, gauge_algs_c4v
+        ) in (
+            minimal ? minimal_combinations_c4v :
+                Iterators.product(spacetypes, scalartypes, projector_algs_c4v, gauge_algs_c4v)
         )
         alg = C4vCTMRG(; tol, projector_alg)
         env_pre, psi = _preconverged_env_c4v(S, T)
diff --git a/test/testsuite/ctmrg/initialization.jl b/test/testsuite/ctmrg/initialization.jl
index 927786119..316bb3dd5 100644
--- a/test/testsuite/ctmrg/initialization.jl
+++ b/test/testsuite/ctmrg/initialization.jl
@@ -48,7 +48,7 @@ function ctmrg_initialization_critical_ising(AT)
         @test info.convergence_error ≤ tol
 
         # specific custom starting product state
-        p_data = ComplexF64[1; 0;;]
+        p_data = adapt(AT, ComplexF64[1; 0;;])
         p = Tensor(p_data, P)
         prod_env0 = ProductStateEnv(reshape([p, p, flip(p, 1), flip(p, 1)], 4, 1, 1))
         env0_custom = initialize_ctmrg_environment(n, ApplicationInitialization(), prod_env0)
diff --git a/test/testsuite/ctmrg/pepo.jl b/test/testsuite/ctmrg/pepo.jl
index 1ed914d57..0258e943d 100644
--- a/test/testsuite/ctmrg/pepo.jl
+++ b/test/testsuite/ctmrg/pepo.jl
@@ -57,10 +57,19 @@ beta = 0.2391 # slightly lower temperature than βc ≈ 0.2216544
 # cover all different flavors
 ctm_styles = [:SequentialCTMRG, :SimultaneousCTMRG]
 projector_algs = [:HalfInfiniteProjector, :FullInfiniteProjector]
-
-function ctmrg_pepo_runthroughs(AT)
-    return @testset "PEPO CTMRG runthroughs for unitcell=$(unitcell) ($AT)" for unitcell in
-        [(1, 1, 1), (1, 1, 2)]
+unitcells = [(1, 1, 1), (1, 1, 2)]
+
+# minimal subset of the combinations tested below which still covers every option value at least once
+minimal_combinations = [
+    ((1, 1, 1), [(:SequentialCTMRG, :FullInfiniteProjector)]),
+    ((1, 1, 2), [(:SimultaneousCTMRG, :HalfInfiniteProjector)]),
+]
+
+function ctmrg_pepo_runthroughs(AT; minimal::Bool = false)
+    return @testset "PEPO CTMRG runthroughs for unitcell=$(unitcell) ($AT)" for (unitcell, algs) in (
+            minimal ? minimal_combinations :
+                [(uc, Iterators.product(ctm_styles, projector_algs)) for uc in unitcells]
+        )
         Random.seed!(81812781144)
 
         O, M, E = three_dimensional_classical_ising(AT; beta)
@@ -77,7 +86,7 @@ function ctmrg_pepo_runthroughs(AT)
 
         @testset "PEPO CTMRG contraction using $alg with $projector_alg" for (
                 alg, projector_alg,
-            ) in Iterators.product(ctm_styles, projector_algs)
+            ) in algs
             env, = leading_boundary(env0, n; alg, maxiter = 150, projector_alg)
         end
     end
diff --git a/test/testsuite/ctmrg/unitcell.jl b/test/testsuite/ctmrg/unitcell.jl
index 848e1ed6e..fb7b24e70 100644
--- a/test/testsuite/ctmrg/unitcell.jl
+++ b/test/testsuite/ctmrg/unitcell.jl
@@ -12,6 +12,8 @@ ctm_algs = [
     SimultaneousCTMRG(; projector_alg = :HalfInfiniteProjector),
     SimultaneousCTMRG(; projector_alg = :FullInfiniteProjector),
 ]
+# minimal subset which still covers both CTMRG and both projector algs
+ctm_algs_minimal = ctm_algs[[1, 4]]
 
 function test_unitcell(
         AT, ctm_alg, unitcell,
@@ -47,9 +49,10 @@ function test_unitcell(
     return nothing
 end
 
-function ctmrg_unitcell_random_cartesian_spaces(AT)
+function ctmrg_unitcell_random_cartesian_spaces(AT; minimal::Bool = false)
     Random.seed!(91283219347)
-    return @testset "Random Cartesian spaces with $ctm_alg ($AT)" for ctm_alg in ctm_algs
+    return @testset "Random Cartesian spaces with $ctm_alg ($AT)" for ctm_alg in
+        (minimal ? ctm_algs_minimal : ctm_algs)
         unitcell = (3, 3)
 
         Pspaces = ComplexSpace.(rand(2:3, unitcell...))
@@ -67,9 +70,10 @@ function ctmrg_unitcell_random_cartesian_spaces(AT)
     end
 end
 
-function ctmrg_unitcell_specific_u1_spaces(AT)
+function ctmrg_unitcell_specific_u1_spaces(AT; minimal::Bool = false)
     Random.seed!(91283219347)
-    return @testset "Specific U1 spaces with $ctm_alg ($AT)" for ctm_alg in ctm_algs
+    return @testset "Specific U1 spaces with $ctm_alg ($AT)" for ctm_alg in
+        (minimal ? ctm_algs_minimal : ctm_algs)
         unitcell = (2, 2)
 
         PA = U1Space(-1 => 1, 0 => 1)
diff --git a/test/testsuite/gradients/c4v_ctmrg_gradients.jl b/test/testsuite/gradients/c4v_ctmrg_gradients.jl
index 4f6cd5db4..41cc4238d 100644
--- a/test/testsuite/gradients/c4v_ctmrg_gradients.jl
+++ b/test/testsuite/gradients/c4v_ctmrg_gradients.jl
@@ -29,6 +29,19 @@ gradient_algs = [[nothing, :FixedPointGradient, :ImplicitGradient]]
 gradient_solver_algs = [[:GeomSum, :ManualIter, :GMRES, :BiCGStab, :Arnoldi]]
 steps = -0.01:0.005:0.01
 
+# minimal subset of (ctmrg_alg, projector_alg, decomposition_rrule_alg, gradient_alg, gradient_solver_alg)
+# combinations which still covers every option value at least once per model
+minimal_combinations = [
+    [
+        (:C4vCTMRG, :C4vEighProjector, :FullPullback, nothing, nothing),
+        (:C4vCTMRG, :C4vQRProjector, :FullPullback, :ImplicitGradient, :GMRES),
+        (:C4vCTMRG, :C4vEighProjector, :TruncPullback, :FixedPointGradient, :GeomSum),
+        (:C4vCTMRG, :C4vQRProjector, :FullPullback, :FixedPointGradient, :ManualIter),
+        (:C4vCTMRG, :C4vEighProjector, :TruncPullback, :FixedPointGradient, :BiCGStab),
+        (:C4vCTMRG, :C4vQRProjector, :FullPullback, :FixedPointGradient, :Arnoldi),
+    ],
+]
+
 # record which rrule alg is compatible with which projector alg
 allowed_rrule_algs = Dict(
     :C4vEighProjector => keys(PEPSKit.EIGH_RRULE_SYMBOLS),
@@ -38,7 +51,7 @@ allowed_rrule_algs = Dict(
 # be selective on which configurations to test the naive gradient for
 naive_gradient_combinations = [(:C4vCTMRG, :C4vEighProjector, :FullPullback), (:C4vCTMRG, :C4vQRProjector, :FullPullback)]
 
-function gradients_c4v(AT)
+function gradients_c4v(AT; minimal::Bool = false)
     naive_gradient_done = Set()
     return @testset "AD C4v CTMRG energy gradients for $(names[i]) model ($AT)" verbose = true for i in
         eachindex(
@@ -54,8 +67,9 @@ function gradients_c4v(AT)
         gsalgs = gradient_solver_algs[i]
         @testset "ctmrg_alg=:$ctmrg_alg, projector_alg=:$projector_alg, decomposition_rrule_alg=:$decomposition_rrule_alg and gradient_alg=(alg = :$gradient_alg, solver_alg = :$gradient_solver_alg)" for (
                 ctmrg_alg, projector_alg, decomposition_rrule_alg, gradient_alg, gradient_solver_alg,
-            ) in Iterators.product(
-                calgs, palgs, dalgs, galgs, gsalgs
+            ) in (
+                minimal ? minimal_combinations[i] :
+                    Iterators.product(calgs, palgs, dalgs, galgs, gsalgs)
             )
 
             # check for allowed algorithm combinations when testing naive gradient
@@ -89,6 +103,7 @@ function gradients_c4v(AT)
             Random.seed!(sd)
             dir = adapt(AT, InfinitePEPS(Pspace, Vspace))
             psi = adapt(AT, InfinitePEPS(Pspace, Vspace))
+            model = adapt(AT, models[i])
             symmetrize!(psi, symmetry)
             symmetrize!(dir, symmetry)
             # instantiate to avoid having to type this twice...
@@ -123,7 +138,7 @@ function gradients_c4v(AT)
                         contrete_ctmrg_alg;
                         alg_rrule = concrete_gradient_alg,
                     )
-                    return cost_function(psi, env2, models[i])
+                    return cost_function(psi, env2, model)
                 end
                 g = only(g)
                 symmetrize!(g, symmetry)
diff --git a/test/testsuite/gradients/ctmrg_gradients.jl b/test/testsuite/gradients/ctmrg_gradients.jl
index c6c9e995b..19b8023f9 100644
--- a/test/testsuite/gradients/ctmrg_gradients.jl
+++ b/test/testsuite/gradients/ctmrg_gradients.jl
@@ -29,6 +29,23 @@ gradient_solver_algs = [
 ]
 steps = -0.01:0.005:0.01
 
+# minimal subset of (ctmrg_alg, projector_alg, svd_rrule_alg, gradient_alg, gradient_solver_alg)
+# combinations which still covers every option value at least once per model
+minimal_combinations = [
+    [
+        (:SimultaneousCTMRG, :HalfInfiniteProjector, :FullPullback, nothing, nothing),
+        (:SimultaneousCTMRG, :HalfInfiniteProjector, :FullPullback, :ImplicitGradient, :GMRES),
+        (:SequentialCTMRG, :FullInfiniteProjector, :TruncPullback, :FixedPointGradient, :GeomSum),
+        (:SimultaneousCTMRG, :FullInfiniteProjector, :Arnoldi, :FixedPointGradient, :ManualIter),
+        (:SequentialCTMRG, :HalfInfiniteProjector, :FullPullback, :FixedPointGradient, :BiCGStab),
+        (:SimultaneousCTMRG, :HalfInfiniteProjector, :Arnoldi, :FixedPointGradient, :Arnoldi),
+    ],
+    [
+        (:SimultaneousCTMRG, :HalfInfiniteProjector, :FullPullback, :ImplicitGradient, :GMRES),
+        (:SequentialCTMRG, :FullInfiniteProjector, :Arnoldi, :FixedPointGradient, :GeomSum),
+    ],
+]
+
 # don't check naive AD gradients for all algorithm combinations, since it's slow
 naive_gradient_combinations = [
     (:SimultaneousCTMRG, :HalfInfiniteProjector, :FullPullback),
@@ -46,7 +63,7 @@ function _check_disallowed_combination(
     return false
 end
 
-function gradients_asymmetric(AT)
+function gradients_asymmetric(AT; minimal::Bool = false)
     naive_gradient_done = Set()
     return @testset "AD CTMRG energy gradients for $(names[i]) model ($AT)" verbose = true for i in
         eachindex(
@@ -62,8 +79,9 @@ function gradients_asymmetric(AT)
         gsalgs = gradient_solver_algs[i]
         @testset "ctmrg_alg=:$ctmrg_alg, projector_alg=:$projector_alg, svd_rrule_alg=:$svd_rrule_alg, gradient_alg=(; alg = :$gradient_alg, solver_alg = (; alg = :$gradient_solver_alg))" for (
                 ctmrg_alg, projector_alg, svd_rrule_alg, gradient_alg, gradient_solver_alg,
-            ) in Iterators.product(
-                calgs, palgs, salgs, galgs, gsalgs
+            ) in (
+                minimal ? minimal_combinations[i] :
+                    Iterators.product(calgs, palgs, salgs, galgs, gsalgs)
             )
 
             # only run GMRES for the implicit gradient, and skip distinction between decomposition rrule algs
@@ -95,6 +113,7 @@ function gradients_asymmetric(AT)
 
             @info "optimtest of ctmrg_alg=:$ctmrg_alg, projector_alg=:$projector_alg, svd_rrule_alg=:$svd_rrule_alg and gradient_alg=(; alg = :$gradient_alg, solver_alg = (; alg = :$gradient_solver_alg)) on $(names[i])"
             Random.seed!(42039482030)
+            model = adapt(AT, models[i])
             dir = adapt(AT, InfinitePEPS(Pspace, Vspace))
             psi = adapt(AT, InfinitePEPS(Pspace, Vspace))
             # instantiate to avoid having to type this twice...
@@ -128,7 +147,7 @@ function gradients_asymmetric(AT)
                         concrete_ctmrg_alg;
                         alg_rrule = concrete_gradient_alg,
                     )
-                    return cost_function(psi, env2, models[i])
+                    return cost_function(psi, env2, model)
                 end
 
                 return E, only(g)
diff --git a/test/testsuite/timeevol/tf_ising_finiteT.jl b/test/testsuite/timeevol/tf_ising_finiteT.jl
index b6e2fce53..51ff21c07 100644
--- a/test/testsuite/timeevol/tf_ising_finiteT.jl
+++ b/test/testsuite/timeevol/tf_ising_finiteT.jl
@@ -19,11 +19,11 @@ function converge_env(state, χ::Int)
     return env
 end
 
-function measure_mag(pepo::InfinitePEPO, env::CTMRGEnv; purified::Bool = false)
+function measure_mag(AT, pepo::InfinitePEPO, env::CTMRGEnv; purified::Bool = false)
     r, c = 1, 1
     lattice = physicalspace(pepo)
-    Mx = LocalOperator(lattice, ((r, c),) => σˣ(Float64, Trivial))
-    Mz = LocalOperator(lattice, ((r, c),) => σᶻ(Float64, Trivial))
+    Mx = LocalOperator(lattice, ((r, c),) => adapt(AT, σˣ(Float64, Trivial)))
+    Mz = LocalOperator(lattice, ((r, c),) => adapt(AT, σᶻ(Float64, Trivial)))
     if purified
         magx = expectation_value(pepo, Mx, pepo, env)
         magz = expectation_value(pepo, Mz, pepo, env)
@@ -61,7 +61,7 @@ function timeevol_ising_finiteT(AT)
         pepo, = gauge_fix(pepo, BPGauge(), bp_env)
 
         env = converge_env(InfinitePartitionFunction(pepo), 16)
-        result_β = measure_mag(pepo, env)
+        result_β = measure_mag(AT, pepo, env)
         @info "tr(σ(x,z)ρ) at T = $(1 / β): $(result_β)."
         @test β ≈ info.t
         @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2)
@@ -70,7 +70,7 @@ function timeevol_ising_finiteT(AT)
         pepo2, = compress((pepo, pepo), LocalTruncation(trunc_pepo))
         normalize!.(pepo2.A)
         env2 = converge_env(InfinitePartitionFunction(pepo2), 16)
-        result_2β = measure_mag(pepo2, env2)
+        result_2β = measure_mag(AT, pepo2, env2)
         @info "tr(σ(x,z)ρ) at T = $(1 / (2β)): $(result_2β)."
         @test isapprox(abs.(result_2β), bm_2β, rtol = 5.0e-3)
 
@@ -78,7 +78,7 @@ function timeevol_ising_finiteT(AT)
         alg = SimpleUpdate(; trunc = trunc_pepo, purified = true, bipartite, force_mpo)
         pepo, wts, info = time_evolve(pepo0, ham, dt, 2 * nstep, alg, wts0; symmetrize_gates)
         env = converge_env(InfinitePEPS(pepo), 8)
-        result_2β′ = measure_mag(pepo, env; purified = true)
+        result_2β′ = measure_mag(AT, pepo, env; purified = true)
         @info "⟨ρ|σ(x,z)|ρ⟩ at T = $(1 / (2β)): $(result_2β′)."
         @test 2 * β ≈ info.t
         @test isapprox(abs.(result_2β′), bm_2β, rtol = 1.0e-2)
diff --git a/test/testsuite/utility/eigh_wrapper.jl b/test/testsuite/utility/eigh_wrapper.jl
index d41fafcdb..24357e49d 100644
--- a/test/testsuite/utility/eigh_wrapper.jl
+++ b/test/testsuite/utility/eigh_wrapper.jl
@@ -15,7 +15,7 @@ function lossfun(A, alg, R = randn(space(A)), trunc = notrunc())
     return real(dot(R, V * V')) + dot(D, D)  # Overlap with random tensor R is gauge-invariant and differentiable
 end
 
-function utility_eigh_wrapper(AT)
+function utility_eigh_wrapper(AT; default_alg = :QRIteration)
     return @testset "eigh_wrapper ($AT)" begin
         dtype = ComplexF64
         n = 20
@@ -28,8 +28,8 @@ function utility_eigh_wrapper(AT)
         R = adapt(AT, randn(space(r)))
         R = 0.5 * (R + R')
 
-        full_alg = EighAdjoint(; fwd_alg = (; alg = :QRIteration), rrule_alg = (; alg = :FullPullback))
-        trunc_alg = EighAdjoint(; fwd_alg = (; alg = :QRIteration), rrule_alg = (; alg = :TruncPullback))
+        full_alg = EighAdjoint(; fwd_alg = (; alg = default_alg), rrule_alg = (; alg = :FullPullback))
+        trunc_alg = EighAdjoint(; fwd_alg = (; alg = default_alg), rrule_alg = (; alg = :TruncPullback))
         iter_alg = EighAdjoint(; fwd_alg = (; alg = :Lanczos), rrule_alg = (; alg = :TruncPullback))
 
         @testset "Non-truncated eigh" begin
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 56da57e98..064b34fa6 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -10,3 +10,15 @@ if !is_buildkite
     TestSuite.timeevol_cluster_identity_gate(Vector)
     TestSuite.timeevol_cluster_hubbard(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.timeevol_cluster_bond_truncation(CuArray)
+    TestSuite.timeevol_cluster_identity_gate(CuArray)
+    TestSuite.timeevol_cluster_hubbard(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.timeevol_cluster_bond_truncation(ROCArray)
+    TestSuite.timeevol_cluster_identity_gate(ROCArray)
+    TestSuite.timeevol_cluster_hubbard(ROCArray)
+end
diff --git a/test/timeevol/j1j2_finiteT.jl b/test/timeevol/j1j2_finiteT.jl
index 856214d70..b9f46c24f 100644
--- a/test/timeevol/j1j2_finiteT.jl
+++ b/test/timeevol/j1j2_finiteT.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -8,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.timeevol_j1j2_finiteT(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.timeevol_j1j2_finiteT(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.timeevol_j1j2_finiteT(ROCArray)
+end
diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl
index 220da6dec..f2b618155 100644
--- a/test/timeevol/sitedep_truncation.jl
+++ b/test/timeevol/sitedep_truncation.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -9,3 +9,13 @@ if !is_buildkite
     TestSuite.timeevol_sitedep_rotation(Vector)
     TestSuite.timeevol_sitedep_su(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.timeevol_sitedep_rotation(CuArray)
+    TestSuite.timeevol_sitedep_su(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.timeevol_sitedep_rotation(ROCArray)
+    TestSuite.timeevol_sitedep_su(ROCArray)
+end
diff --git a/test/timeevol/tf_ising_finiteT.jl b/test/timeevol/tf_ising_finiteT.jl
index 6c6e1c74e..72347dd89 100644
--- a/test/timeevol/tf_ising_finiteT.jl
+++ b/test/timeevol/tf_ising_finiteT.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -8,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.timeevol_ising_finiteT(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.timeevol_ising_finiteT(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.timeevol_ising_finiteT(ROCArray)
+end
diff --git a/test/timeevol/timestep.jl b/test/timeevol/timestep.jl
index 8d426f5cd..a3692114c 100644
--- a/test/timeevol/timestep.jl
+++ b/test/timeevol/timestep.jl
@@ -1,4 +1,4 @@
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -8,3 +8,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.timeevol_timestep(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.timeevol_timestep(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.timeevol_timestep(ROCArray)
+end
diff --git a/test/toolbox/densitymatrices.jl b/test/toolbox/densitymatrices.jl
index 3030ec187..97f9199b2 100644
--- a/test/toolbox/densitymatrices.jl
+++ b/test/toolbox/densitymatrices.jl
@@ -1,5 +1,5 @@
 using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -12,3 +12,13 @@ if !is_buildkite
     TestSuite.toolbox_densitymatrix_too_many_layers(Vector)
     TestSuite.toolbox_densitymatrix_generic_fallback(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.toolbox_single_layer_densitymatrix(CuArray)
+    TestSuite.toolbox_double_layer_densitymatrix(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.toolbox_single_layer_densitymatrix(ROCArray)
+    TestSuite.toolbox_double_layer_densitymatrix(ROCArray)
+end
diff --git a/test/utility/correlator.jl b/test/utility/correlator.jl
index 94dec8796..829fa6df6 100644
--- a/test/utility/correlator.jl
+++ b/test/utility/correlator.jl
@@ -1,5 +1,5 @@
 using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -11,3 +11,15 @@ if !is_buildkite
     TestSuite.utility_correlator_purified_ipepo(Vector)
     TestSuite.utility_correlator_single_layer_ipepo(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.utility_correlator_infinite_peps(CuArray)
+    TestSuite.utility_correlator_purified_ipepo(CuArray)
+    TestSuite.utility_correlator_single_layer_ipepo(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.utility_correlator_infinite_peps(ROCArray)
+    TestSuite.utility_correlator_purified_ipepo(ROCArray)
+    TestSuite.utility_correlator_single_layer_ipepo(ROCArray)
+end
diff --git a/test/utility/eigh_wrapper.jl b/test/utility/eigh_wrapper.jl
index 74435d9e3..4ba4d8fe0 100644
--- a/test/utility/eigh_wrapper.jl
+++ b/test/utility/eigh_wrapper.jl
@@ -1,5 +1,5 @@
 using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -9,3 +9,12 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.utility_eigh_wrapper(Vector)
 end
+
+if CUDA.functional()
+    # CUSOLVER doesn't provide QRIteration for eigh
+    TestSuite.utility_eigh_wrapper(CuArray; default_alg = :DivideAndConquer)
+end
+
+if AMDGPU.functional()
+    TestSuite.utility_eigh_wrapper(ROCArray)
+end
diff --git a/test/utility/retractions.jl b/test/utility/retractions.jl
index 24ee96426..93e570310 100644
--- a/test/utility/retractions.jl
+++ b/test/utility/retractions.jl
@@ -1,5 +1,5 @@
 using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -9,3 +9,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.utility_retractions(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.utility_retractions(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.utility_retractions(ROCArray)
+end
diff --git a/test/utility/svd_wrapper.jl b/test/utility/svd_wrapper.jl
index 8abda1b93..b13497c59 100644
--- a/test/utility/svd_wrapper.jl
+++ b/test/utility/svd_wrapper.jl
@@ -1,5 +1,5 @@
 using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -9,3 +9,11 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true"
 if !is_buildkite
     TestSuite.utility_svd_wrapper(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.utility_svd_wrapper(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.utility_svd_wrapper(ROCArray)
+end
diff --git a/test/utility/symmetrization.jl b/test/utility/symmetrization.jl
index 6b22bce8b..143b0f2e8 100644
--- a/test/utility/symmetrization.jl
+++ b/test/utility/symmetrization.jl
@@ -1,5 +1,5 @@
 using Test
-using PEPSKit
+using PEPSKit, CUDA, AMDGPU
 
 @isdefined(TestSuite) || include("../testsuite/TestSuite.jl")
 using .TestSuite
@@ -12,3 +12,17 @@ if !is_buildkite
     TestSuite.utility_symmetrization_rotate(Vector)
     TestSuite.utility_symmetrization_rotate_reflect(Vector)
 end
+
+if CUDA.functional()
+    TestSuite.utility_symmetrization_reflect_depth(CuArray)
+    TestSuite.utility_symmetrization_reflect_width(CuArray)
+    TestSuite.utility_symmetrization_rotate(CuArray)
+    TestSuite.utility_symmetrization_rotate_reflect(CuArray)
+end
+
+if AMDGPU.functional()
+    TestSuite.utility_symmetrization_reflect_depth(ROCArray)
+    TestSuite.utility_symmetrization_reflect_width(ROCArray)
+    TestSuite.utility_symmetrization_rotate(ROCArray)
+    TestSuite.utility_symmetrization_rotate_reflect(ROCArray)
+end