diff --git a/Project.toml b/Project.toml
index ee80b10..c93b6e5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BlockTensorKit"
 uuid = "5f87ffc2-9cf1-4a46-8172-465d160bd8cd"
+version = "0.3.2"
 authors = ["Lukas Devos <ldevos98@gmail.com> and contributors"]
-version = "0.3.1"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
@@ -21,11 +21,11 @@ BlockArrays = "1"
 Combinatorics = "1"
 Compat = "4.13"
 LinearAlgebra = "1"
-MatrixAlgebraKit = "0.5"
+MatrixAlgebraKit = "0.6"
 Random = "1"
 SafeTestsets = "0.1"
 Strided = "2"
-TensorKit = "0.15.2"
+TensorKit = "0.16"
 TensorOperations = "5"
 Test = "1"
 TestExtras = "0.2, 0.3"
diff --git a/src/linalg/factorizations.jl b/src/linalg/factorizations.jl
index fcd0e80..347754f 100644
--- a/src/linalg/factorizations.jl
+++ b/src/linalg/factorizations.jl
@@ -1,7 +1,6 @@
 using MatrixAlgebraKit
 using MatrixAlgebraKit: AbstractAlgorithm, YALAPACK.BlasMat, Algorithm
 import MatrixAlgebraKit as MAK
-using TensorKit.Factorizations: @check_space, @check_scalar
 
 # Type piracy for defining the MAK rules on BlockArrays!
 # -----------------------------------------------------
@@ -16,70 +15,74 @@ function MatrixAlgebraKit.one!(A::BlockBlasMat)
     return A
 end
 
-for f in (
-        :svd_compact, :svd_full, :svd_vals, :qr_compact, :qr_full, :qr_null,
-        :lq_compact, :lq_full, :lq_null, :eig_full, :eig_vals, :eigh_full,
-        :eigh_vals, :left_polar, :right_polar,
-    )
+for f in
+    [
+        :svd_compact, :svd_full, :svd_vals,
+        :qr_compact, :qr_full, :qr_null,
+        :lq_compact, :lq_full, :lq_null,
+        :eig_full, :eig_vals, :eigh_full, :eigh_vals,
+        :left_polar, :right_polar,
+        :project_hermitian, :project_antihermitian, :project_isometric,
+    ]
     f! = Symbol(f, :!)
-    @eval MAK.$f!(t::BlockBlasMat, F, alg::AbstractAlgorithm) = $f!(Array(t), alg)
-    @eval MAK.$f!(t::BlockBlasMat, F, alg::MAK.DiagonalAlgorithm) = error("Not diagonal")
+    @eval MAK.default_algorithm(::typeof($f!), ::Type{T}; kwargs...) where {T <: AbstractBlockTensorMap} =
+        MAK.default_algorithm($f!, eltype(T); kwargs...)
 end
 
-# disambiguations
-for (f!, Alg) in (
-        (:lq_compact!, :LAPACK_HouseholderLQ), (:lq_full!, :LAPACK_HouseholderLQ), (:lq_null!, :LAPACK_HouseholderLQ),
-        (:lq_compact!, :LQViaTransposedQR), (:lq_full!, :LQViaTransposedQR), (:lq_null!, :LQViaTransposedQR),
-        (:qr_compact!, :LAPACK_HouseholderQR), (:qr_full!, :LAPACK_HouseholderQR), (:qr_null!, :LAPACK_HouseholderQR),
-        (:svd_compact!, :LAPACK_SVDAlgorithm), (:svd_full!, :LAPACK_SVDAlgorithm), (:svd_vals!, :LAPACK_SVDAlgorithm),
-        (:eig_full!, :LAPACK_EigAlgorithm), (:eig_trunc!, :TruncatedAlgorithm), (:eig_vals!, :LAPACK_EigAlgorithm),
-        (:eigh_full!, :LAPACK_EighAlgorithm), (:eigh_trunc!, :TruncatedAlgorithm), (:eigh_vals!, :LAPACK_EighAlgorithm),
-        (:left_polar!, :PolarViaSVD), (:right_polar!, :PolarViaSVD),
+for f! in (
+        :qr_compact!, :qr_full!, :lq_compact!, :lq_full!,
+        :eig_full!, :eigh_full!, :svd_compact!, :svd_full!,
+        :left_polar!, :right_polar!,
     )
-    @eval MAK.$f!(t::BlockBlasMat, F, alg::MAK.$Alg) = $f!(Array(t), alg)
-end
-
-const GPU_QRAlgorithm = Union{MAK.CUSOLVER_HouseholderQR, MAK.ROCSOLVER_HouseholderQR}
-for f! in (:qr_compact!, :qr_full!, :qr_null!)
-    @eval MAK.$f!(t::BlockBlasMat, QR, alg::GPU_QRAlgorithm) = error()
+    @eval function MAK.$f!(t::AbstractBlockTensorMap, F, alg::AbstractAlgorithm)
+        TensorKit.foreachblock(t, F...) do _, (tblock, Fblocks...)
+            Fblocks′ = MAK.$f!(Array(tblock), alg)
+            # deal with the case where the output is not in-place
+            for (b′, b) in zip(Fblocks′, Fblocks)
+                b === b′ || copy!(b, b′)
+            end
+            return nothing
+        end
+        return F
+    end
 end
 
-for (f!, Alg) in (
-        (:eigh_full!, :GPU_EighAlgorithm), (:eigh_vals!, :GPU_EighAlgorithm),
-        (:eig_full!, :GPU_EigAlgorithm), (:eig_vals!, :GPU_EigAlgorithm),
-        (:svd_full!, :GPU_SVDAlgorithm), (:svd_compact!, :GPU_SVDAlgorithm), (:svd_vals!, :GPU_SVDAlgorithm),
+# Handle these separately because single output instead of tuple
+for f! in (
+        :qr_null!, :lq_null!,
+        :svd_vals!, :eig_vals!, :eigh_vals!,
+        :project_hermitian!, :project_antihermitian!, :project_isometric!,
     )
-    @eval MAK.$f!(t::BlockBlasMat, F, alg::MAK.$Alg) = error()
+    @eval function MAK.$f!(t::AbstractBlockTensorMap, N, alg::AbstractAlgorithm)
+        TensorKit.foreachblock(t, N) do _, (tblock, Nblock)
+            Nblock′ = MAK.$f!(Array(tblock), alg)
+            # deal with the case where the output is not the same as the input
+            Nblock === Nblock′ || copy!(Nblock, Nblock′)
+            return nothing
+        end
+        return N
+    end
 end
 
-
-for f in (:qr, :lq, :eig, :eigh, :gen_eig, :svd, :polar)
-    default_f_algorithm = Symbol(:default_, f, :_algorithm)
-    @eval MAK.$default_f_algorithm(::Type{<:BlockBlasMat{T}}; kwargs...) where {T} =
-        MAK.$default_f_algorithm(Matrix{T}; kwargs...)
+# 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
-
-# Make sure sparse blocktensormaps have dense outputs
-function MAK.check_input(::typeof(qr_full!), t::AbstractBlockTensorMap, QR, ::AbstractAlgorithm)
-    Q, R = QR
-
-    # type checks
-    @assert Q isa AbstractTensorMap
-    @assert R isa AbstractTensorMap
-
-    # scalartype checks
-    @check_scalar Q t
-    @check_scalar R t
-
-    # space checks
-    V_Q = ⊕(fuse(codomain(t)))
-    @check_space(Q, codomain(t) ← V_Q)
-    @check_space(R, V_Q ← domain(t))
-
-    return nothing
+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)`
 end
-MAK.check_input(::typeof(qr_full!), t::AbstractBlockTensorMap, QR, ::DiagonalAlgorithm) = error()
 
+# Make sure sparse blocktensormaps have dense outputs
 function MAK.initialize_output(::typeof(qr_full!), t::AbstractBlockTensorMap, ::AbstractAlgorithm)
     V_Q = ⊕(fuse(codomain(t)))
     Q = dense_similar(t, codomain(t) ← V_Q)
@@ -99,26 +102,6 @@ function MAK.initialize_output(::typeof(qr_null!), t::AbstractBlockTensorMap, ::
     return N
 end
 
-function MAK.check_input(::typeof(lq_full!), t::AbstractBlockTensorMap, LQ, ::AbstractAlgorithm)
-    L, Q = LQ
-
-    # type checks
-    @assert L isa AbstractTensorMap
-    @assert Q isa AbstractTensorMap
-
-    # scalartype checks
-    @check_scalar L t
-    @check_scalar Q t
-
-    # space checks
-    V_Q = ⊕(fuse(domain(t)))
-    @check_space(L, codomain(t) ← V_Q)
-    @check_space(Q, V_Q ← domain(t))
-
-    return nothing
-end
-MAK.check_input(::typeof(lq_full!), t::AbstractBlockTensorMap, LQ, ::DiagonalAlgorithm) = error()
-
 function MAK.initialize_output(::typeof(lq_full!), t::AbstractBlockTensorMap, ::AbstractAlgorithm)
     V_Q = ⊕(fuse(domain(t)))
     L = dense_similar(t, codomain(t) ← V_Q)
@@ -138,99 +121,18 @@ function MAK.initialize_output(::typeof(lq_null!), t::AbstractBlockTensorMap, ::
     return N
 end
 
-function MAK.check_input(::typeof(MAK.left_orth_polar!), t::AbstractBlockTensorMap, WP, ::AbstractAlgorithm)
-    codomain(t) ≿ domain(t) ||
-        throw(ArgumentError("Polar decomposition requires `codomain(t) ≿ domain(t)`"))
-
-    W, P = WP
-    @assert W isa AbstractTensorMap
-    @assert P isa AbstractTensorMap
-
-    # scalartype checks
-    @check_scalar W t
-    @check_scalar P t
-
-    # space checks
-    VW = ⊕(fuse(domain(t)))
-    @check_space(W, codomain(t) ← VW)
-    @check_space(P, VW ← domain(t))
-
-    return nothing
-end
 function MAK.initialize_output(::typeof(left_polar!), t::AbstractBlockTensorMap, ::AbstractAlgorithm)
     W = dense_similar(t, space(t))
     P = dense_similar(t, domain(t) ← domain(t))
     return W, P
 end
 
-function MAK.check_input(::typeof(MAK.right_orth_polar!), t::AbstractBlockTensorMap, PWᴴ, ::AbstractAlgorithm)
-    codomain(t) ≾ domain(t) ||
-        throw(ArgumentError("Polar decomposition requires `domain(t) ≿ codomain(t)`"))
-
-    P, Wᴴ = PWᴴ
-    @assert P isa AbstractTensorMap
-    @assert Wᴴ isa AbstractTensorMap
-
-    # scalartype checks
-    @check_scalar P t
-    @check_scalar Wᴴ t
-
-    # space checks
-    VW = ⊕(fuse(codomain(t)))
-    @check_space(P, codomain(t) ← VW)
-    @check_space(Wᴴ, VW ← domain(t))
-
-    return nothing
-end
 function MAK.initialize_output(::typeof(right_polar!), t::AbstractBlockTensorMap, ::AbstractAlgorithm)
     P = dense_similar(t, codomain(t) ← codomain(t))
     Wᴴ = dense_similar(t, space(t))
     return P, Wᴴ
 end
 
-function MAK.initialize_output(::typeof(left_null!), t::AbstractBlockTensorMap, ::AbstractAlgorithm)
-    V_Q = infimum(fuse(codomain(t)), fuse(domain(t)))
-    V_N = ⊖(fuse(codomain(t)), V_Q)
-    return dense_similar(t, codomain(t) ← V_N)
-end
-function MAK.initialize_output(::typeof(right_null!), t::AbstractBlockTensorMap, ::AbstractAlgorithm)
-    V_Q = infimum(fuse(codomain(t)), fuse(domain(t)))
-    V_N = ⊖(fuse(domain(t)), V_Q)
-    return dense_similar(t, V_N ← domain(t))
-end
-
-function MAK.check_input(::typeof(eigh_full!), t::AbstractBlockTensorMap, DV, ::AbstractAlgorithm)
-    domain(t) == codomain(t) ||
-        throw(ArgumentError("Eigenvalue decomposition requires square input tensor"))
-
-    D, V = DV
-
-    # type checks
-    @assert D isa DiagonalTensorMap
-    @assert V isa AbstractTensorMap
-
-    # scalartype checks
-    @check_scalar D t real
-    @check_scalar V t
-
-    # space checks
-    V_D = ⊕(fuse(domain(t)))
-    @check_space(D, V_D ← V_D)
-    @check_space(V, codomain(t) ← V_D)
-
-    return nothing
-end
-MAK.check_input(::typeof(eigh_full!), t::AbstractBlockTensorMap, DV, ::DiagonalAlgorithm) = error()
-
-function MAK.check_input(::typeof(eigh_vals!), t::AbstractBlockTensorMap, D, ::AbstractAlgorithm)
-    @check_scalar D t real
-    @assert D isa DiagonalTensorMap
-    V_D = ⊕(fuse(domain(t)))
-    @check_space(D, V_D ← V_D)
-    return nothing
-end
-MAK.check_input(::typeof(eigh_vals!), t::AbstractBlockTensorMap, D, ::DiagonalAlgorithm) = error()
-
 function MAK.initialize_output(::typeof(eigh_full!), t::AbstractBlockTensorMap, ::AbstractAlgorithm)
     V_D = ⊕(fuse(domain(t)))
     T = real(scalartype(t))
@@ -239,30 +141,6 @@ function MAK.initialize_output(::typeof(eigh_full!), t::AbstractBlockTensorMap,
     return D, V
 end
 
-
-function MAK.check_input(::typeof(eig_full!), t::AbstractBlockTensorMap, DV, ::AbstractAlgorithm)
-    domain(t) == codomain(t) ||
-        throw(ArgumentError("Eigenvalue decomposition requires square input tensor"))
-
-    D, V = DV
-
-    # type checks
-    @assert D isa DiagonalTensorMap
-    @assert V isa AbstractTensorMap
-
-    # scalartype checks
-    @check_scalar D t complex
-    @check_scalar V t complex
-
-    # space checks
-    V_D = ⊕(fuse(domain(t)))
-    @check_space(D, V_D ← V_D)
-    @check_space(V, codomain(t) ← V_D)
-
-    return nothing
-end
-MAK.check_input(::typeof(eig_full!), t::AbstractBlockTensorMap, DV, ::DiagonalAlgorithm) = error()
-
 function MAK.initialize_output(::typeof(eig_full!), t::AbstractBlockTensorMap, ::AbstractAlgorithm)
     V_D = ⊕(fuse(domain(t)))
     Tc = complex(scalartype(t))
@@ -271,28 +149,6 @@ function MAK.initialize_output(::typeof(eig_full!), t::AbstractBlockTensorMap, :
     return D, V
 end
 
-function MAK.check_input(::typeof(svd_full!), t::AbstractBlockTensorMap, USVᴴ, ::AbstractAlgorithm)
-    U, S, Vᴴ = USVᴴ
-
-    # type checks
-    @assert U isa AbstractTensorMap
-    @assert S isa AbstractTensorMap
-    @assert Vᴴ isa AbstractTensorMap
-
-    # scalartype checks
-    @check_scalar U t
-    @check_scalar S t real
-    @check_scalar Vᴴ t
-
-    # space checks
-    V_cod = ⊕(fuse(codomain(t)))
-    V_dom = ⊕(fuse(domain(t)))
-    @check_space(U, codomain(t) ← V_cod)
-    @check_space(S, V_cod ← V_dom)
-    @check_space(Vᴴ, V_dom ← domain(t))
-
-    return nothing
-end
 function MAK.initialize_output(::typeof(svd_full!), t::AbstractBlockTensorMap, ::AbstractAlgorithm)
     V_cod = ⊕(fuse(codomain(t)))
     V_dom = ⊕(fuse(domain(t)))
@@ -312,8 +168,13 @@ end
 # Disambiguate Diagonal implementations
 # -------------------------------------
 # these shouldn't ever happen as blocktensors aren't diagonal
-for f! in (:eig_full!, :eigh_full!, :lq_compact!, :lq_full!, :qr_compact!, :qr_full!, :svd_full!)
-    @eval function MAK.initialize_output(::typeof($f!), t::AbstractBlockTensorMap, alg::DiagonalAlgorithm)
+for f! in (
+        :eig_full!, :eig_vals!, :eigh_full!, :eigh_vals!,
+        :lq_compact!, :lq_full!, :qr_compact!, :qr_full!,
+        :svd_full!, :svd_compact!, :svd_vals!,
+    )
+    @eval MAK.initialize_output(::typeof($f!), t::AbstractBlockTensorMap, ::DiagonalAlgorithm) =
+        error("Blocktensors are incompatible with diagonal algorithm")
+    @eval MAK.$f!(::AbstractBlockTensorMap, x, ::DiagonalAlgorithm) =
         error("Blocktensors are incompatible with diagonal algorithm")
-    end
 end
diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl
index 00c2558..56ba6ee 100644
--- a/src/tensors/tensoroperations.jl
+++ b/src/tensors/tensoroperations.jl
@@ -154,6 +154,80 @@ function TO.tensoralloc(
     return C
 end
 
+# unfortunate overlaod until TK fix
+function TK.blas_contract!(
+        C::AbstractBlockTensorMap,
+        A::AbstractBlockTensorMap, pA::Index2Tuple,
+        B::AbstractBlockTensorMap, pB::Index2Tuple,
+        pAB::Index2Tuple, α, β,
+        backend, allocator
+    )
+    bstyle = BraidingStyle(sectortype(C))
+    bstyle isa SymmetricBraiding ||
+        throw(SectorMismatch("only tensors with symmetric braiding rules can be contracted; try `@planar` instead"))
+    TC = scalartype(C)
+
+    # check which tensors have to be permuted/copied
+    copyA = !(TO.isblascontractable(A, pA) && eltype(A) === TC)
+    copyB = !(TO.isblascontractable(B, pB) && eltype(B) === TC)
+
+    if bstyle isa Fermionic && any(isdual ∘ Base.Fix1(space, B), pB[1])
+        # twist smallest object if neither or both already have to be permuted
+        # otherwise twist the one that already is copied
+        if !(copyA ⊻ copyB)
+            twistA = dim(A) < dim(B)
+        else
+            twistA = copyA
+        end
+        twistB = !twistA
+        copyA |= twistA
+        copyB |= twistB
+    else
+        twistA = false
+        twistB = false
+    end
+
+    # Bring A in the correct form for BLAS contraction
+    if copyA
+        Anew = TO.tensoralloc_add(TC, A, pA, false, Val(true), allocator)
+        Anew = TO.tensoradd!(Anew, A, pA, false, One(), Zero(), backend, allocator)
+        twistA && twist!(Anew, filter(!isdual ∘ Base.Fix1(space, Anew), domainind(Anew)))
+    else
+        Anew = permute(A, pA)
+    end
+    pAnew = (codomainind(Anew), domainind(Anew))
+
+    # Bring B in the correct form for BLAS contraction
+    if copyB
+        Bnew = TO.tensoralloc_add(TC, B, pB, false, Val(true), allocator)
+        Bnew = TO.tensoradd!(Bnew, B, pB, false, One(), Zero(), backend, allocator)
+        twistB && twist!(Bnew, filter(isdual ∘ Base.Fix1(space, Bnew), codomainind(Bnew)))
+    else
+        Bnew = permute(B, pB)
+    end
+    pBnew = (codomainind(Bnew), domainind(Bnew))
+
+    # Bring C in the correct form for BLAS contraction
+    ipAB = TO.oindABinC(pAB, pAnew, pBnew)
+    copyC = !TO.isblasdestination(C, ipAB)
+
+    if copyC
+        Cnew = TO.tensoralloc_add(TC, C, ipAB, false, Val(true), allocator)
+        mul!(Cnew, Anew, Bnew)
+        TO.tensoradd!(C, Cnew, pAB, false, α, β, backend, allocator)
+        TO.tensorfree!(Cnew, allocator)
+    else
+        Cnew = permute(C, ipAB)
+        mul!(Cnew, Anew, Bnew, α, β)
+    end
+
+    copyA && TO.tensorfree!(Anew, allocator)
+    copyB && TO.tensorfree!(Bnew, allocator)
+
+    return C
+end
+
+
 # tensorfree!
 # -----------
 function TO.tensorfree!(t::BlockTensorMap, allocator = TO.DefaultAllocator())
diff --git a/test/linalg/factorizations.jl b/test/linalg/factorizations.jl
index 2321bc1..dce9da3 100644
--- a/test/linalg/factorizations.jl
+++ b/test/linalg/factorizations.jl
@@ -45,18 +45,18 @@ for V in spacelist
 
                 Q, R = @constinferred qr_compact(t)
                 @test Q * R ≈ t
-                @test isisometry(Q)
+                @test isisometric(Q)
 
-                Q, R = @constinferred left_orth(t; kind = :qr)
+                Q, R = @constinferred left_orth(t; alg = :qr)
                 @test Q * R ≈ t
-                @test isisometry(Q)
+                @test isisometric(Q)
 
                 N = @constinferred qr_null(t)
-                @test isisometry(N)
+                @test isisometric(N)
                 @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
 
-                N = @constinferred left_null(t; kind = :qr)
-                @test isisometry(N)
+                N = @constinferred left_null(t; alg = :qr)
+                @test isisometric(N)
                 @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
             end
 
@@ -71,12 +71,12 @@ for V in spacelist
 
                 Q, R = @constinferred qr_compact(t)
                 @test Q * R ≈ t
-                @test isisometry(Q)
+                @test isisometric(Q)
                 @test dim(Q) == dim(R) == dim(t)
 
-                Q, R = @constinferred left_orth(t; kind = :qr)
+                Q, R = @constinferred left_orth(t; alg = :qr)
                 @test Q * R ≈ t
-                @test isisometry(Q)
+                @test isisometric(Q)
                 @test dim(Q) == dim(R) == dim(t)
 
                 N = @constinferred qr_null(t)
@@ -97,14 +97,14 @@ for V in spacelist
 
                 L, Q = @constinferred lq_compact(t)
                 @test L * Q ≈ t
-                @test isisometry(Q; side = :right)
+                @test isisometric(Q; side = :right)
 
-                L, Q = @constinferred right_orth(t; kind = :lq)
+                L, Q = @constinferred right_orth(t; alg = :lq)
                 @test L * Q ≈ t
-                @test isisometry(Q; side = :right)
+                @test isisometric(Q; side = :right)
 
                 Nᴴ = @constinferred lq_null(t)
-                @test isisometry(Nᴴ; side = :right)
+                @test isisometric(Nᴴ; side = :right)
                 @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t))
             end
 
@@ -119,12 +119,12 @@ for V in spacelist
 
                 L, Q = @constinferred lq_compact(t)
                 @test L * Q ≈ t
-                @test isisometry(Q; side = :right)
+                @test isisometric(Q; side = :right)
                 @test dim(Q) == dim(L) == dim(t)
 
-                L, Q = @constinferred right_orth(t; kind = :lq)
+                L, Q = @constinferred right_orth(t; alg = :lq)
                 @test L * Q ≈ t
-                @test isisometry(Q; side = :right)
+                @test isisometric(Q; side = :right)
                 @test dim(Q) == dim(L) == dim(t)
 
                 Nᴴ = @constinferred lq_null(t)
@@ -142,12 +142,12 @@ for V in spacelist
                 @assert domain(t) ≾ codomain(t)
                 w, p = @constinferred left_polar(t)
                 @test w * p ≈ t
-                @test isisometry(w)
+                @test isisometric(w)
                 # @test isposdef(p)
 
-                w, p = @constinferred left_orth(t; kind = :polar)
+                w, p = @constinferred left_orth(t; alg = :polar)
                 @test w * p ≈ t
-                @test isisometry(w)
+                @test isisometric(w)
             end
 
             for T in eltypes,
@@ -156,12 +156,12 @@ for V in spacelist
                 @assert codomain(t) ≾ domain(t)
                 p, wᴴ = @constinferred right_polar(t)
                 @test p * wᴴ ≈ t
-                @test isisometry(wᴴ; side = :right)
+                @test isisometric(wᴴ; side = :right)
                 # @test isposdef(p)
 
-                p, wᴴ = @constinferred right_orth(t; kind = :polar)
+                p, wᴴ = @constinferred right_orth(t; alg = :polar)
                 @test p * wᴴ ≈ t
-                @test isisometry(wᴴ; side = :right)
+                @test isisometric(wᴴ; side = :right)
             end
         end
 
@@ -180,25 +180,25 @@ for V in spacelist
 
                 u, s, vᴴ = @constinferred svd_compact(t)
                 @test u * s * vᴴ ≈ t
-                @test isisometry(u)
+                @test isisometric(u)
                 # @test isposdef(s)
-                @test isisometry(vᴴ; side = :right)
+                @test isisometric(vᴴ; side = :right)
 
                 s′ = LinearAlgebra.diag(s)
-                for (c, b) in LinearAlgebra.svdvals(t)
+                for (c, b) in pairs(LinearAlgebra.svdvals(t))
                     @test b ≈ s′[c]
                 end
 
-                v, c = @constinferred left_orth(t; kind = :svd)
+                v, c = @constinferred left_orth(t; alg = :svd)
                 @test v * c ≈ t
-                @test isisometry(v)
+                @test isisometric(v)
 
-                N = @constinferred left_null(t; kind = :svd)
-                @test isisometry(N)
+                N = @constinferred left_null(t; alg = :svd)
+                @test isisometric(N)
                 @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t))
 
-                Nᴴ = @constinferred right_null(t; kind = :svd)
-                @test isisometry(Nᴴ; side = :right)
+                Nᴴ = @constinferred right_null(t; alg = :svd)
+                @test isisometric(Nᴴ; side = :right)
                 @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t))
             end
 
@@ -227,22 +227,22 @@ for V in spacelist
 
                 U, S, Vᴴ = @constinferred svd_trunc(t; trunc = notrunc())
                 @test U * S * Vᴴ ≈ t
-                @test isisometry(U)
-                @test isisometry(Vᴴ; side = :right)
+                @test isisometric(U)
+                @test isisometric(Vᴴ; side = :right)
 
                 trunc = truncrank(dim(domain(S)) ÷ 2)
                 U1, S1, Vᴴ1 = @constinferred svd_trunc(t; trunc)
                 @test t * Vᴴ1' ≈ U1 * S1
-                @test isisometry(U1)
-                @test isisometry(Vᴴ1; side = :right)
+                @test isisometric(U1)
+                @test isisometric(Vᴴ1; side = :right)
                 @test dim(domain(S1)) <= trunc.howmany
 
                 λ = minimum(minimum, values(LinearAlgebra.diag(S1)))
                 trunc = trunctol(; atol = λ - 10eps(λ))
                 U2, S2, Vᴴ2 = @constinferred svd_trunc(t; trunc)
                 @test t * Vᴴ2' ≈ U2 * S2
-                @test isisometry(U2)
-                @test isisometry(Vᴴ2; side = :right)
+                @test isisometric(U2)
+                @test isisometric(Vᴴ2; side = :right)
                 @test minimum(minimum, values(LinearAlgebra.diag(S1))) >= λ
                 @test U2 ≈ U1
                 @test S2 ≈ S1
@@ -251,15 +251,15 @@ for V in spacelist
                 trunc = truncspace(space(S2, 1))
                 U3, S3, Vᴴ3 = @constinferred svd_trunc(t; trunc)
                 @test t * Vᴴ3' ≈ U3 * S3
-                @test isisometry(U3)
-                @test isisometry(Vᴴ3; side = :right)
+                @test isisometric(U3)
+                @test isisometric(Vᴴ3; side = :right)
                 @test space(S3, 1) ≾ space(S2, 1)
 
                 trunc = truncerror(; atol = 0.5)
                 U4, S4, Vᴴ4 = @constinferred svd_trunc(t; trunc)
                 @test t * Vᴴ4' ≈ U4 * S4
-                @test isisometry(U4)
-                @test isisometry(Vᴴ4; side = :right)
+                @test isisometric(U4)
+                @test isisometric(Vᴴ4; side = :right)
                 @test norm(t - U4 * S4 * Vᴴ4) <= 0.5
             end
         end
@@ -274,7 +274,7 @@ for V in spacelist
                 @test t * v ≈ v * d
 
                 d′ = LinearAlgebra.diag(d)
-                for (c, b) in LinearAlgebra.eigvals(t)
+                for (c, b) in pairs(LinearAlgebra.eigvals(t))
                     @test sort(b; by = abs) ≈ sort(d′[c]; by = abs)
                 end
 
@@ -289,7 +289,7 @@ for V in spacelist
 
                 t2 = (t + t')
                 D, V = eigen(t2)
-                @test isisometry(V)
+                @test isisometric(V)
                 D̃, Ṽ = @constinferred eigh_full(t2)
                 @test D ≈ D̃
                 @test V ≈ Ṽ