diff --git a/src/linalg/factorizations.jl b/src/linalg/factorizations.jl
index 94cf57f..ce887f5 100644
--- a/src/linalg/factorizations.jl
+++ b/src/linalg/factorizations.jl
@@ -74,22 +74,16 @@ for f in [
         MAK.$f!(t, out, MAK.select_algorithm(MAK.$f!, t, nothing; alg.kwargs...))
 end
 
-# specializations until fixes in base package
-function MAK.is_left_isometric(A::BlockMatrix; atol::Real = 0, rtol::Real = MAK.defaulttol(A), norm = LinearAlgebra.norm)
-    P = A' * A
-    nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
-    for I in MAK.diagind(P)
-        P[I] -= 1
-    end
-    return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
-end
-function MAK.is_right_isometric(A::BlockMatrix; atol::Real = 0, rtol::Real = MAK.defaulttol(A), norm = LinearAlgebra.norm)
-    P = A * A'
-    nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
-    for I in MAK.diagind(P)
-        P[I] -= 1
-    end
-    return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
+# specializations until fixes in base packages:
+# - `similar(::MulAdd{<:AbstractBlockLayout,...})` hardcodes `Array` blocks, so `A' * A` on
+#   non-CPU blocks mixes host and device storage (JuliaArrays/BlockArrays.jl#215).
+# - `MAK.diagview` silently does not alias into a `BlockMatrix`, so subtracting the identity there is a no-op.
+# Densifying sidesteps both, and preserves the storage type of the blocks.
+function MAK.is_left_isometric(A::BlockMatrix; kwargs...)
+    return MAK.is_left_isometric(copy_dense!(similar_dense(A), A); kwargs...)
+end
+function MAK.is_right_isometric(A::BlockMatrix; kwargs...)
+    return MAK.is_right_isometric(copy_dense!(similar_dense(A), A); kwargs...)
 end
 
 # Make sure sparse blocktensormaps have dense outputs
diff --git a/test/abstracttensor/braidingtensor.jl b/test/abstracttensor/braidingtensor.jl
index 14ac8c8..4ddaabc 100644
--- a/test/abstracttensor/braidingtensor.jl
+++ b/test/abstracttensor/braidingtensor.jl
@@ -39,10 +39,10 @@ end
 @testset "Issue #68" begin
     x = randn(ComplexF64, (⊞(ℙ^1) ⊗ ⊞(ℙ^1)) ← ⊞(ℙ^1))
     O = randn(ComplexF64, ((ℙ^1 ⊞ ℙ^4 ⊞ ℙ^16 ⊞ ℙ^4) ⊗ ⊞(ℙ^2)) ← (⊞(ℙ^2) ⊗ ⊞(ℙ^1)))
-    A = randn(ComplexF64, (ℙ^4 ⊗ F^2 ⊗ (ℙ^2)') ← ℙ^1)
+    A = randn(ComplexF64, (ℙ^4 ⊗ ℙ^2 ⊗ (ℙ^2)') ← ℙ^1)
     Ab = randn(ComplexF64, (ℙ^4 ⊗ ℙ^2 ⊗ (ℙ^2)') ← ℙ^1)
 
     @plansor y[-1 -2; -3] ≔ A[-1 4 2; 1] * O[-2 6; 4 5] * τ[5 7; 2 3] *
         conj(Ab[-3 6 7; 8]) * x[1 3; 8]
-    @test space(y) == (⊞(ℙ^4) ⊗ ⊞(ℙ^1)) ← ⊞(ℙ^4)
+    @test space(y) == ((⊞(ℙ^4) ⊗ (ℙ^1 ⊞ ℙ^4 ⊞ ℙ^16 ⊞ ℙ^4)) ← ⊞(ℙ^4))
 end