diff --git a/src/algorithms/ctmrg_gauge_fix.jl b/src/algorithms/ctmrg_gauge_fix.jl
index 1c548d64..980ec185 100644
--- a/src/algorithms/ctmrg_gauge_fix.jl
+++ b/src/algorithms/ctmrg_gauge_fix.jl
@@ -177,7 +177,6 @@ function fix_global_phases(envprev::CTMRGEnv, envfix::CTMRGEnv)
     return CTMRGEnv(cornersgfix, edgesgfix)
 end
 
-
 function calc_convergence(envs, CSold, TSold)
     CSnew = map(x -> tsvd(x; alg=TensorKit.SVD())[2], envs.corners)
     ΔCS = maximum(zip(CSold, CSnew)) do (c_old, c_new)
diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/peps_opt.jl
index b8310a5e..1dca18ef 100644
--- a/src/algorithms/peps_opt.jl
+++ b/src/algorithms/peps_opt.jl
@@ -101,6 +101,9 @@ struct PEPSOptimize{G}
         if gradient_alg isa GradMode
             if S == :sequential && G.parameters[1] == :fixed
                 throw(ArgumentError(":sequential and :fixed are not compatible"))
+            elseif boundary_alg.projector_alg.svd_alg.fwd_alg isa IterSVD &&
+                G.parameters[1] == :fixed
+                throw(ArgumentError("IterSVD and :fixed are currently not compatible"))
             end
         end
         return new{G}(boundary_alg, optimizer, reuse_env, gradient_alg)
diff --git a/test/ctmrg/fixed_iterscheme.jl b/test/ctmrg/fixed_iterscheme.jl
index 57194674..1cfc9f83 100644
--- a/test/ctmrg/fixed_iterscheme.jl
+++ b/test/ctmrg/fixed_iterscheme.jl
@@ -13,7 +13,7 @@ using PEPSKit:
 # initialize parameters
 χbond = 2
 χenv = 16
-svd_algs = [SVDAdjoint(; fwd_alg=TensorKit.SVD()), SVDAdjoint(; fwd_alg=IterSVD())]
+svd_algs = [SVDAdjoint(; fwd_alg=TensorKit.SVD())]  #, SVDAdjoint(; fwd_alg=IterSVD())]
 unitcells = [(1, 1), (3, 4)]
 
 # test for element-wise convergence after application of fixed step
@@ -31,7 +31,7 @@ unitcells = [(1, 1), (3, 4)]
     # do extra iteration to get SVD
     env_conv2, info = ctmrg_iter(psi, env_conv1, ctm_alg)
     env_fix, signs = gauge_fix(env_conv1, env_conv2)
-    @test calc_elementwise_convergence(env_conv1, env_fix) ≈ 0 atol=1e-6
+    @test calc_elementwise_convergence(env_conv1, env_fix) ≈ 0 atol = 1e-6
 
     # fix gauge of SVD
     U_fix, V_fix = fix_relative_phases(info.U, info.V, signs)
@@ -43,36 +43,93 @@ unitcells = [(1, 1), (3, 4)]
     # do iteration with FixedSVD
     env_fixedsvd, = ctmrg_iter(psi, env_conv1, ctm_alg_fix)
     env_fixedsvd = fix_global_phases(env_conv1, env_fixedsvd)
-    @test calc_elementwise_convergence(env_conv1, env_fixedsvd) ≈ 0 atol=1e-6
+    @test calc_elementwise_convergence(env_conv1, env_fixedsvd) ≈ 0 atol = 1e-6
 end
 
-# TODO: Why doesn't fixed work with IterSVD?
-##
-# ctm_alg = CTMRG(;
-#     tol=1e-12,
-#     miniter=4,
-#     maxiter=100,
-#     verbosity=1,
-#     ctmrgscheme=:simultaneous,
-#     svd_alg=SVDAdjoint(; fwd_alg=IterSVD()),
-# )
-
-# # initialize states
-# Random.seed!(91283219347)
-# psi = InfinitePEPS(2, χbond)
-# env_conv1 = leading_boundary(CTMRGEnv(psi, ComplexSpace(χenv)), psi, ctm_alg);
-
-# # do extra iteration to get SVD
-# env_conv2, info = ctmrg_iter(psi, env_conv1, ctm_alg);
-# env_fix, signs = gauge_fix(env_conv1, env_conv2);
-# @test calc_elementwise_convergence(env_conv1, env_fix) ≈ 0 atol=1e-6
-
-# # fix gauge of SVD
-# U_fix, V_fix = fix_relative_phases(info.U, info.V, signs);
-# svd_alg_fix = SVDAdjoint(; fwd_alg=FixedSVD(U_fix, info.S, V_fix));
-# ctm_alg_fix = CTMRG(; svd_alg=svd_alg_fix, trscheme=notrunc(), ctmrgscheme=:simultaneous);
-
-# # do iteration with FixedSVD
-# env_fixedsvd, = ctmrg_iter(psi, env_conv1, ctm_alg_fix);
-# env_fixedsvd = fix_global_phases(env_conv1, env_fixedsvd);
-# @test calc_elementwise_convergence(env_conv1, env_fixedsvd) ≈ 0 atol=1e-6
\ No newline at end of file
+# TODO: Why doesn't FixedSVD work with previous U, S and V from IterSVD?
+@testset "Element-wise consistency of TensorKit.SVD and IterSVD" begin
+    ctm_alg_iter = CTMRG(;
+        tol=1e-12,
+        verbosity=2,
+        ctmrgscheme=:simultaneous,
+        svd_alg=SVDAdjoint(; fwd_alg=IterSVD()),
+    )
+    ctm_alg_full = CTMRG(;
+        tol=1e-12,
+        verbosity=2,
+        ctmrgscheme=:simultaneous,
+        svd_alg=SVDAdjoint(; fwd_alg=TensorKit.SVD()),
+    )
+
+    # initialize states
+    Random.seed!(91283219347)
+    psi = InfinitePEPS(2, χbond)
+    env_init = CTMRGEnv(psi, ComplexSpace(χenv))
+    env_conv1 = leading_boundary(env_init, psi, ctm_alg_full)
+
+    # do extra iteration to get SVD
+    env_conv2_iter, info_iter = ctmrg_iter(psi, env_conv1, ctm_alg_iter)
+    env_fix_iter, signs_iter = gauge_fix(env_conv1, env_conv2_iter)
+    @test calc_elementwise_convergence(env_conv1, env_fix_iter) ≈ 0 atol = 1e-6
+
+    env_conv2_full, info_full = ctmrg_iter(psi, env_conv1, ctm_alg_full)
+    env_fix_full, signs_full = gauge_fix(env_conv1, env_conv2_full)
+    @test calc_elementwise_convergence(env_conv1, env_fix_full) ≈ 0 atol = 1e-6
+
+    # fix gauge of SVD
+    U_fix_iter, V_fix_iter = fix_relative_phases(info_iter.U, info_iter.V, signs_iter)
+    svd_alg_fix_iter = SVDAdjoint(; fwd_alg=FixedSVD(U_fix_iter, info_iter.S, V_fix_iter))
+    ctm_alg_fix_iter = CTMRG(;
+        svd_alg=svd_alg_fix_iter, trscheme=notrunc(), ctmrgscheme=:simultaneous
+    )
+
+    U_fix_full, V_fix_full = fix_relative_phases(info_full.U, info_full.V, signs_full)
+    svd_alg_fix_full = SVDAdjoint(; fwd_alg=FixedSVD(U_fix_full, info_full.S, V_fix_full))
+    ctm_alg_fix_full = CTMRG(;
+        svd_alg=svd_alg_fix_full, trscheme=notrunc(), ctmrgscheme=:simultaneous
+    )
+
+    # do iteration with FixedSVD
+    env_fixedsvd_iter, = ctmrg_iter(psi, env_conv1, ctm_alg_fix_iter)
+    env_fixedsvd_iter = fix_global_phases(env_conv1, env_fixedsvd_iter)
+    # @test calc_elementwise_convergence(env_conv1, env_fixedsvd_iter) ≈ 0 atol = 1e-6  # This should work, but doesn't!
+
+    env_fixedsvd_full, = ctmrg_iter(psi, env_conv1, ctm_alg_fix_full)
+    env_fixedsvd_full = fix_global_phases(env_conv1, env_fixedsvd_full)
+    @test calc_elementwise_convergence(env_conv1, env_fixedsvd_full) ≈ 0 atol = 1e-6
+
+    # check matching decompositions
+    atol = 1e-12
+    decomposition_check = all(
+        zip(info_iter.U, info_iter.S, info_iter.V, info_full.U, info_full.S, info_full.V),
+    ) do (U_iter, S_iter, V_iter, U_full, S_full, V_full)
+        diff = U_iter * S_iter * V_iter - U_full * S_full * V_full
+        all(x -> isapprox(abs(x), 0; atol), diff.data)
+    end
+    @test decomposition_check
+
+    # check matching singular values
+    svalues_check = all(zip(info_iter.S, info_full.S)) do (S_iter, S_full)
+        diff = S_iter - S_full
+        all(x -> isapprox(abs(x), 0; atol), diff.data)
+    end
+    @test svalues_check
+
+    # check normalization of U's and V's
+    Us = [info_iter.U, U_fix_iter, info_full.U, U_fix_full]
+    Vs = [info_iter.V, V_fix_iter, info_full.V, V_fix_full]
+    for (U, V) in zip(Us, Vs)
+        U_check = all(U) do u
+            uu = u' * u
+            diff = uu - id(space(uu, 1))
+            all(x -> isapprox(abs(x), 0; atol), diff.data)
+        end
+        @test U_check
+        V_check = all(V) do v
+            vv = v * v'
+            diff = vv - id(space(vv, 1))
+            all(x -> isapprox(abs(x), 0; atol), diff.data)
+        end
+        @test V_check
+    end
+end
diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl
index 724b1261..122e42fc 100644
--- a/test/ctmrg/gaugefix.jl
+++ b/test/ctmrg/gaugefix.jl
@@ -30,7 +30,7 @@ schemes = [:simultaneous, :sequential]
     ctm = leading_boundary(ctm, psi, alg)
     ctm2, = ctmrg_iter(psi, ctm, alg)
     ctm_fixed, = gauge_fix(ctm, ctm2)
-    @test calc_elementwise_convergence(ctm, ctm_fixed) ≈ 0 atol=1e-6
+    @test calc_elementwise_convergence(ctm, ctm_fixed) ≈ 0 atol = 1e-6
 end
 
 @testset "Z2 symmetry ($T) - ($unitcell) - ($ctmrgscheme)" for (T, unitcell, ctmrgscheme) in
@@ -52,5 +52,5 @@ end
     ctm = leading_boundary(ctm, psi, alg)
     ctm2, = ctmrg_iter(psi, ctm, alg)
     ctm_fixed, = gauge_fix(ctm, ctm2)
-    @test calc_elementwise_convergence(ctm, ctm_fixed) ≈ 0 atol=1e-6
+    @test calc_elementwise_convergence(ctm, ctm_fixed) ≈ 0 atol = 1e-6
 end