From 66678b656586e2997665f750f4141e4b4c3809d4 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 18:30:30 -0400
Subject: [PATCH 01/17] Open the v0.21 breaking round

Accumulates the contract and matricize interface redesign. The version
stays at 0.21.0-DEV until the release PR strips the suffix.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 146f410..9cf7ac4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "TensorAlgebra"
 uuid = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
-version = "0.20.1"
+version = "0.21.0-DEV"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]

From 4a63a19b1ebd2fcc1277caee89e0a4aa1dc442af Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 18:36:30 -0400
Subject: [PATCH 02/17] Delete the out-of-place unmatricizeperm

It had no callers outside its own forwarding method, and its behaviour is
covered by the in-place form. Also raises the subproject compat bounds the
round-opening bump missed.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 docs/Project.toml     |  2 +-
 examples/Project.toml |  2 +-
 src/matricize.jl      | 29 +----------------------------
 test/Project.toml     |  2 +-
 test/test_basics.jl   | 14 +++-----------
 5 files changed, 7 insertions(+), 42 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index e02c595..d9d7449 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -11,4 +11,4 @@ path = ".."
 Documenter = "1.8.1"
 ITensorFormatter = "0.2.27"
 Literate = "2.20.1"
-TensorAlgebra = "0.20"
+TensorAlgebra = "0.21"
diff --git a/examples/Project.toml b/examples/Project.toml
index b7e5f73..194d5b6 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -5,4 +5,4 @@ TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 path = ".."
 
 [compat]
-TensorAlgebra = "0.20"
+TensorAlgebra = "0.21"
diff --git a/src/matricize.jl b/src/matricize.jl
index 972cb07..5ede4d8 100644
--- a/src/matricize.jl
+++ b/src/matricize.jl
@@ -226,7 +226,7 @@ end
 # domain groups, given codomain-facing (un-dualized), the same convention as `similar_map`. A
 # matricize style stores the domain axes dualized, so its overload re-dualizes them with `conj`
 # (a no-op on a dense axis). This is the primary overload point for new matricize styles.
-# Permutation is handled separately by `unmatricizeperm`, so `unmatricize` never has to
+# Permutation is handled separately by `unmatricizeperm!`, so `unmatricize` never has to
 # disambiguate axis tuples from permutation tuples regardless of how unconstrained `m` and the
 # axes are.
 function unmatricize(style::MatricizeStyle, m, axes_codomain, axes_domain)
@@ -245,33 +245,6 @@ function bipartition_axes(t::Tuple, split...)
     return axes_codomain, conj.(axes_domain)
 end
 
-# Inverse-bipermutation form: split `axes_dest` into codomain/domain groups reordered by the
-# inverse bipermutation, unmatricize in that order, then permute back.
-function unmatricizeperm(
-        m, axes_dest,
-        invperm_codomain::Tuple{Vararg{Int}}, invperm_domain::Tuple{Vararg{Int}}
-    )
-    return unmatricizeperm(
-        MatricizeStyle(m),
-        m,
-        axes_dest,
-        invperm_codomain,
-        invperm_domain
-    )
-end
-function unmatricizeperm(
-        style::MatricizeStyle, m, axes_dest,
-        invperm_codomain::Tuple{Vararg{Int}}, invperm_domain::Tuple{Vararg{Int}}
-    )
-    invbiperm = BiTuple(invperm_codomain, invperm_domain)
-    length(axes_dest) == length(invbiperm) ||
-        throw(ArgumentError("axes do not match permutation"))
-    axes_codomain, axes_domain = bipartition_axes(axes_dest, invbiperm)
-    a12 = unmatricize(style, m, axes_codomain, axes_domain)
-    biperm_dest = BiTuple(Tuple(invperm(invbiperm)), Val(length_codomain(invbiperm)))
-    return bipermutedims(a12, biperm_dest)
-end
-
 function unmatricizeperm!(
         a_dest, m,
         invperm_codomain::Tuple{Vararg{Int}}, invperm_domain::Tuple{Vararg{Int}}
diff --git a/test/Project.toml b/test/Project.toml
index 08c7d26..0ba4cc3 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -37,7 +37,7 @@ Random = "1.10"
 SafeTestsets = "0.1"
 StableRNGs = "1.0.2"
 Suppressor = "0.2"
-TensorAlgebra = "0.20"
+TensorAlgebra = "0.21"
 TensorKit = "0.17"
 TensorOperations = "5.1.4"
 Test = "1.10"
diff --git a/test/test_basics.jl b/test/test_basics.jl
index 9e285f8..d5f7686 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -3,7 +3,7 @@ using EllipsisNotation: var".."
 using StableRNGs: StableRNG
 using TensorAlgebra: BiTuple, ContractAlgorithm, bipermutedims, bipermutedims!, contract,
     contract!, contractadd!, length_codomain, length_domain, matricizeperm, unmatricize,
-    unmatricizeperm, unmatricizeperm!
+    unmatricizeperm!
 using TensorOperations: TensorOperations
 using Test: @test, @test_broken, @test_throws, @testset
 
@@ -130,25 +130,18 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         @test eltype(a) === elt
         @test a ≈ a0
 
-        a = unmatricizeperm(m, axes0, (1, 2), (3, 4))
-        @test eltype(a) === elt
-        @test a ≈ a0
-
         perm_codomain = (4, 2)
         perm_domain = (1, 3)
         invperm_codomain = (3, 2)
         invperm_domain = (4, 1)
         perm = (4, 2, 1, 3)
-        a = unmatricizeperm(m, map(i -> axes0[i], perm), invperm_codomain, invperm_domain)
-        @test eltype(a) === elt
-        @test a ≈ permutedims(a0, perm)
-
         a = similar(a0)
         unmatricizeperm!(a, m, (1, 2), (3, 4))
         @test a ≈ a0
 
         m1 = matricizeperm(a0, perm_codomain, perm_domain)
-        a = unmatricizeperm(m1, axes0, perm_codomain, perm_domain)
+        a = similar(a0)
+        unmatricizeperm!(a, m1, perm_codomain, perm_domain)
         @test a ≈ a0
 
         a1 = permutedims(a0, perm)
@@ -174,7 +167,6 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         @test a isa Array{elt, 0}
         @test a[] == m[1, 1]
 
-        @test_throws ArgumentError unmatricizeperm(m, (), (1, 2), (3,))
         @test_throws ArgumentError unmatricizeperm!(m, m, (1, 2), (3,))
     end
 

From 17521648f7bd6de2ff40b21ed270d4e59f32eed0 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 18:44:13 -0400
Subject: [PATCH 03/17] Name the unmatricize bipermutation for what it is

The argument maps the destination's dimension order to the matrix's; it is
not intrinsically an inverse, so invperm_ described one caller's derivation
rather than the parameter.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/matricize.jl            | 25 +++++++++++++++----------
 test/test_matricizestyle.jl |  4 ++--
 2 files changed, 17 insertions(+), 12 deletions(-)

diff --git a/src/matricize.jl b/src/matricize.jl
index 5ede4d8..10465a5 100644
--- a/src/matricize.jl
+++ b/src/matricize.jl
@@ -215,10 +215,10 @@ end
 ismatricizeview(::MatricizeStyle, a, ndims_codomain::Val) = false
 function ismatricizeview(
         style::MatricizeStyle, a,
-        invperm_codomain::Tuple{Vararg{Int}}, invperm_domain::Tuple{Vararg{Int}}
+        perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
     )
-    isidentityperm((invperm_codomain..., invperm_domain...)) || return false
-    return ismatricizeview(style, a, Val(length(invperm_codomain)))
+    isidentityperm((perm_codomain..., perm_domain...)) || return false
+    return ismatricizeview(style, a, Val(length(perm_codomain)))
 end
 
 # ====================================  unmatricize  =======================================
@@ -245,22 +245,27 @@ function bipartition_axes(t::Tuple, split...)
     return axes_codomain, conj.(axes_domain)
 end
 
+# The bipermutation maps the destination's dimension order to the matrix's: `axes(a_dest)` grouped
+# by it gives the legs in `m`'s order, and the result is permuted back by its inverse. It is not
+# intrinsically an inverse permutation — the matricized-contraction destination path happens to
+# derive it as `invperm(biperm_dest)`, while a `matricizeperm`/`unmatricizeperm!` round trip passes
+# the same forward bipermutation to both.
 function unmatricizeperm!(
         a_dest, m,
-        invperm_codomain::Tuple{Vararg{Int}}, invperm_domain::Tuple{Vararg{Int}}
+        perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
     )
-    return unmatricizeperm!(MatricizeStyle(m), a_dest, m, invperm_codomain, invperm_domain)
+    return unmatricizeperm!(MatricizeStyle(m), a_dest, m, perm_codomain, perm_domain)
 end
 function unmatricizeperm!(
         style::MatricizeStyle, a_dest, m,
-        invperm_codomain::Tuple{Vararg{Int}}, invperm_domain::Tuple{Vararg{Int}}
+        perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
     )
-    invbiperm = BiTuple(invperm_codomain, invperm_domain)
-    ndims(a_dest) == length(invbiperm) ||
+    biperm_src = BiTuple(perm_codomain, perm_domain)
+    ndims(a_dest) == length(biperm_src) ||
         throw(ArgumentError("destination does not match permutation"))
-    axes_codomain, axes_domain = bipartition_axes(axes(a_dest), invbiperm)
+    axes_codomain, axes_domain = bipartition_axes(axes(a_dest), biperm_src)
     a_perm = unmatricize(style, m, axes_codomain, axes_domain)
-    biperm_dest = BiTuple(Tuple(invperm(invbiperm)), Val(length_codomain(invbiperm)))
+    biperm_dest = BiTuple(Tuple(invperm(biperm_src)), Val(length_codomain(biperm_src)))
     return bipermutedims!(a_dest, a_perm, biperm_dest)
 end
 
diff --git a/test/test_matricizestyle.jl b/test/test_matricizestyle.jl
index c43627f..c7851e5 100644
--- a/test/test_matricizestyle.jl
+++ b/test/test_matricizestyle.jl
@@ -18,10 +18,10 @@ module MatricizeStyleTestUtils
     end
     function TA.unmatricizeperm!(
             ::MyArrayMatricize, a_dest::MyArray, m,
-            invperm_codomain::Tuple{Vararg{Int}}, invperm_domain::Tuple{Vararg{Int}}
+            perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
         )
         TA.unmatricizeperm!(
-            TA.ReshapeMatricize(), a_dest.parent, m, invperm_codomain, invperm_domain
+            TA.ReshapeMatricize(), a_dest.parent, m, perm_codomain, perm_domain
         )
         return a_dest
     end

From 3a9139dc5244fa0e35a338eca4092de73a25066d Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 18:50:18 -0400
Subject: [PATCH 04/17] Merge unmatricizeperm! into unmatricize!

The two shapes are disjoint at the trailing argument, a Val split spec
against a pair of permutation tuples, so one name carries both and the
perm marker stops earning its place.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/contract/contract_matricize.jl |  4 ++--
 src/matricize.jl                   | 17 +++++++++--------
 test/test_basics.jl                | 10 +++++-----
 test/test_matricizestyle.jl        |  4 ++--
 4 files changed, 18 insertions(+), 17 deletions(-)

diff --git a/src/contract/contract_matricize.jl b/src/contract/contract_matricize.jl
index 9e04f5f..1593d73 100644
--- a/src/contract/contract_matricize.jl
+++ b/src/contract/contract_matricize.jl
@@ -34,13 +34,13 @@ function contractopadd!(
         # overwrites `a_dest` in full.
         a_dest_mat = a1_mat * a2_mat
         isone(α) || scale!(a_dest_mat, α)
-        unmatricizeperm!(output_style, a_dest, a_dest_mat, invperm_codomain, invperm_domain)
+        unmatricize!(output_style, a_dest, a_dest_mat, invperm_codomain, invperm_domain)
     else
         # `a_dest`'s data contributes through `β`, so gather it, multiply into the gathered
         # copy, and scatter back.
         a_dest_mat = matricizecopy(output_style, a_dest, invperm_codomain, invperm_domain)
         mul!(a_dest_mat, a1_mat, a2_mat, α, β)
-        unmatricizeperm!(output_style, a_dest, a_dest_mat, invperm_codomain, invperm_domain)
+        unmatricize!(output_style, a_dest, a_dest_mat, invperm_codomain, invperm_domain)
     end
     return a_dest
 end
diff --git a/src/matricize.jl b/src/matricize.jl
index 10465a5..b7e6711 100644
--- a/src/matricize.jl
+++ b/src/matricize.jl
@@ -226,7 +226,8 @@ end
 # domain groups, given codomain-facing (un-dualized), the same convention as `similar_map`. A
 # matricize style stores the domain axes dualized, so its overload re-dualizes them with `conj`
 # (a no-op on a dense axis). This is the primary overload point for new matricize styles.
-# Permutation is handled separately by `unmatricizeperm!`, so `unmatricize` never has to
+# Permutation is handled by the bipermutation form of `unmatricize!`, so out-of-place `unmatricize`
+# never has to
 # disambiguate axis tuples from permutation tuples regardless of how unconstrained `m` and the
 # axes are.
 function unmatricize(style::MatricizeStyle, m, axes_codomain, axes_domain)
@@ -248,15 +249,15 @@ end
 # The bipermutation maps the destination's dimension order to the matrix's: `axes(a_dest)` grouped
 # by it gives the legs in `m`'s order, and the result is permuted back by its inverse. It is not
 # intrinsically an inverse permutation — the matricized-contraction destination path happens to
-# derive it as `invperm(biperm_dest)`, while a `matricizeperm`/`unmatricizeperm!` round trip passes
+# derive it as `invperm(biperm_dest)`, while a `matricizeperm`/`unmatricize!` round trip passes
 # the same forward bipermutation to both.
-function unmatricizeperm!(
+function unmatricize!(
         a_dest, m,
         perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
     )
-    return unmatricizeperm!(MatricizeStyle(m), a_dest, m, perm_codomain, perm_domain)
+    return unmatricize!(MatricizeStyle(m), a_dest, m, perm_codomain, perm_domain)
 end
-function unmatricizeperm!(
+function unmatricize!(
         style::MatricizeStyle, a_dest, m,
         perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
     )
@@ -269,14 +270,14 @@ function unmatricizeperm!(
     return bipermutedims!(a_dest, a_perm, biperm_dest)
 end
 
-# In-place split-axes counterpart of `unmatricize`, as `unmatricizeperm!` is of `unmatricizeperm`:
+# In-place counterpart of `unmatricize`:
 # scatter the fused matrix `m` back into `a_dest`'s existing storage across the codomain/domain
-# split at `ndims_codomain`. The split applies no permutation, so this is `unmatricizeperm!` at the
+# split at `ndims_codomain`. The split applies no permutation, so this is the bipermutation form at the
 # trivial bipermutation, reusing its in-place block scatter (no intermediate `unmatricize` copy).
 function unmatricize!(style::MatricizeStyle, a_dest, m, ndims_codomain::Val)
     K = unval(ndims_codomain)
     N = ndims(a_dest)
-    return unmatricizeperm!(
+    return unmatricize!(
         style,
         a_dest,
         m,
diff --git a/test/test_basics.jl b/test/test_basics.jl
index d5f7686..64d4eee 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -3,7 +3,7 @@ using EllipsisNotation: var".."
 using StableRNGs: StableRNG
 using TensorAlgebra: BiTuple, ContractAlgorithm, bipermutedims, bipermutedims!, contract,
     contract!, contractadd!, length_codomain, length_domain, matricizeperm, unmatricize,
-    unmatricizeperm!
+    unmatricize!
 using TensorOperations: TensorOperations
 using Test: @test, @test_broken, @test_throws, @testset
 
@@ -136,17 +136,17 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         invperm_domain = (4, 1)
         perm = (4, 2, 1, 3)
         a = similar(a0)
-        unmatricizeperm!(a, m, (1, 2), (3, 4))
+        unmatricize!(a, m, (1, 2), (3, 4))
         @test a ≈ a0
 
         m1 = matricizeperm(a0, perm_codomain, perm_domain)
         a = similar(a0)
-        unmatricizeperm!(a, m1, perm_codomain, perm_domain)
+        unmatricize!(a, m1, perm_codomain, perm_domain)
         @test a ≈ a0
 
         a1 = permutedims(a0, perm)
         a = similar(a1)
-        unmatricizeperm!(a, m, invperm_codomain, invperm_domain)
+        unmatricize!(a, m, invperm_codomain, invperm_domain)
         @test a ≈ a1
 
         a = unmatricize(reshape(a0, 1, 120), (), axes0)
@@ -167,7 +167,7 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         @test a isa Array{elt, 0}
         @test a[] == m[1, 1]
 
-        @test_throws ArgumentError unmatricizeperm!(m, m, (1, 2), (3,))
+        @test_throws ArgumentError unmatricize!(m, m, (1, 2), (3,))
     end
 
     @testset "contract eltype widens like a matrix product" begin
diff --git a/test/test_matricizestyle.jl b/test/test_matricizestyle.jl
index c7851e5..038ac12 100644
--- a/test/test_matricizestyle.jl
+++ b/test/test_matricizestyle.jl
@@ -16,11 +16,11 @@ module MatricizeStyleTestUtils
     function TA.matricizecopy(::MyArrayMatricize, a::MyArray, ndims_codomain::Val)
         return TA.matricizecopy(TA.ReshapeMatricize(), a.parent, ndims_codomain)
     end
-    function TA.unmatricizeperm!(
+    function TA.unmatricize!(
             ::MyArrayMatricize, a_dest::MyArray, m,
             perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
         )
-        TA.unmatricizeperm!(
+        TA.unmatricize!(
             TA.ReshapeMatricize(), a_dest.parent, m, perm_codomain, perm_domain
         )
         return a_dest

From 00da2dc4d40224a6ff44423110ef49bd6b5c029b Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 19:00:41 -0400
Subject: [PATCH 05/17] Report an unusable contraction keyword by name

The entry points collect trailing keywords and forward them to the
resolver, whose methods declared none, so any unrecognized keyword
surfaced as a MethodError on an internal function.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/contract/contractalgorithm.jl | 35 ++++++++++++++++++++++++-------
 test/test_basics.jl               | 18 ++++++++++++++++
 2 files changed, 46 insertions(+), 7 deletions(-)

diff --git a/src/contract/contractalgorithm.jl b/src/contract/contractalgorithm.jl
index 02fb018..c172ccf 100644
--- a/src/contract/contractalgorithm.jl
+++ b/src/contract/contractalgorithm.jl
@@ -23,17 +23,38 @@ Base.@kwdef struct TensorOperationsAlgorithm{Backend, Allocator} <: ContractAlgo
     allocator::Allocator = nothing
 end
 
-function select_contract_algorithm(algorithm, a1, a2)
-    return error("Not implemented.")
+# The contraction entry points collect trailing keywords and forward them here, so these accept
+# `kwargs...` even though no `ContractAlgorithm` is configurable by keyword yet. Without it an
+# unrecognized keyword surfaces as a `MethodError` on this internal function rather than as a
+# complaint about the keyword the caller actually passed.
+function reject_algorithm_kwargs(algorithm; kwargs...)
+    isempty(kwargs) && return nothing
+    names = join(map(k -> "`$k`", collect(keys(kwargs))), ", ")
+    return throw(
+        ArgumentError(
+            "unsupported keyword argument(s) $names for contraction algorithm `$(nameof(typeof(algorithm)))`"
+        )
+    )
 end
-function select_contract_algorithm(algorithm::ContractAlgorithm, a1, a2)
+
+function select_contract_algorithm(algorithm, a1, a2; kwargs...)
+    return throw(
+        ArgumentError(
+            "`$algorithm` is not a contraction algorithm; pass a `ContractAlgorithm` as `alg`"
+        )
+    )
+end
+function select_contract_algorithm(algorithm::ContractAlgorithm, a1, a2; kwargs...)
+    reject_algorithm_kwargs(algorithm; kwargs...)
     return algorithm
 end
-function select_contract_algorithm(algorithm::DefaultContractAlgorithm, a1, a2)
-    return default_contract_algorithm(a1, a2)
+function select_contract_algorithm(algorithm::DefaultContractAlgorithm, a1, a2; kwargs...)
+    return default_contract_algorithm(a1, a2; kwargs...)
 end
-function default_contract_algorithm(a1, a2)
-    return default_contract_algorithm(typeof(a1), typeof(a2))
+function default_contract_algorithm(a1, a2; kwargs...)
+    algorithm = default_contract_algorithm(typeof(a1), typeof(a2))
+    reject_algorithm_kwargs(algorithm; kwargs...)
+    return algorithm
 end
 function default_contract_algorithm(A1::Type{<:AbstractArray}, A2::Type{<:AbstractArray})
     return Matricize(MatricizeStyle(MatricizeStyle(A1), MatricizeStyle(A2)))
diff --git a/test/test_basics.jl b/test/test_basics.jl
index 64d4eee..0d789da 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -170,6 +170,24 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         @test_throws ArgumentError unmatricize!(m, m, (1, 2), (3,))
     end
 
+    @testset "contraction algorithm selection rejects unusable keywords" begin
+        a1 = randn(2, 3)
+        a2 = randn(3, 4)
+        # A keyword no algorithm can consume must name itself, not surface as a `MethodError`
+        # from inside the resolver.
+        @test_throws ArgumentError contract((1, 3), a1, (1, 2), a2, (2, 3); nonsense = 1)
+        @test_throws ArgumentError contract(
+            (1, 3), a1, (1, 2), a2, (2, 3); alg = TensorAlgebra.Matricize(), nonsense = 1
+        )
+        # A non-algorithm passed as `alg` says so rather than erroring with "Not implemented".
+        @test_throws ArgumentError TensorAlgebra.select_contract_algorithm(:nope, a1, a2)
+        # The supported spellings still work.
+        @test contract((1, 3), a1, (1, 2), a2, (2, 3)) ≈ a1 * a2
+        @test contract(
+            (1, 3), a1, (1, 2), a2, (2, 3); alg = TensorAlgebra.Matricize()
+        ) ≈ a1 * a2
+    end
+
     @testset "contract eltype widens like a matrix product" begin
         a1 = ones(Bool, (2, 2))
         a2 = ones(Bool, (2, 2))

From 5df5d28a982d64e24d1bc841578a82b5a13a6be1 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 19:28:50 -0400
Subject: [PATCH 06/17] Drop the unused Ellipsis bipermutation spelling

Nothing passed `..`, and supporting it cost a dependency plus a
typed/untyped method tier whose only job was normalizing it. Names the
joint predicate `isbiperm` and routes the three duplicate validations
through `check_biperm`.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 Project.toml                    |  2 --
 src/bituple.jl                  |  5 ++++
 src/contract/allocate_output.jl |  2 +-
 src/matricize.jl                | 48 +++------------------------------
 test/Project.toml               |  2 --
 test/test_basics.jl             |  1 -
 6 files changed, 10 insertions(+), 50 deletions(-)

diff --git a/Project.toml b/Project.toml
index 9cf7ac4..30d47b0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,7 +7,6 @@ authors = ["ITensor developers <support@itensor.org> and contributors"]
 projects = ["benchmark", "dev", "docs", "examples", "test"]
 
 [deps]
-EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -29,7 +28,6 @@ TensorAlgebraTensorKitSectorsExt = "TensorKitSectors"
 TensorAlgebraTensorOperationsExt = "TensorOperations"
 
 [compat]
-EllipsisNotation = "1.8"
 LinearAlgebra = "1.10"
 MatrixAlgebraKit = "0.6"
 Mooncake = "0.4.202, 0.5"
diff --git a/src/bituple.jl b/src/bituple.jl
index 17acf4c..0501193 100644
--- a/src/bituple.jl
+++ b/src/bituple.jl
@@ -70,3 +70,8 @@ function bipartition(t::Tuple, group1::Tuple, group2::Tuple)
 end
 # Split `t` by the two groups of a `BiTuple`.
 bipartition(t::Tuple, bt::BiTuple) = bipartition(t, bt.t1, bt.t2)
+
+# Whether `perm_codomain` and `perm_domain` together permute `1:n`, i.e. whether they are a valid
+# bipartitioned permutation. The two halves only make sense jointly, so this takes them as a pair
+# rather than leaving every caller to splat and call `isperm`.
+isbiperm(perm_codomain, perm_domain) = isperm((perm_codomain..., perm_domain...))
diff --git a/src/contract/allocate_output.jl b/src/contract/allocate_output.jl
index c09f671..13b392c 100644
--- a/src/contract/allocate_output.jl
+++ b/src/contract/allocate_output.jl
@@ -1,7 +1,7 @@
 function check_biperm(a, perm_codomain, perm_domain)
     ndims(a) == length(perm_codomain) + length(perm_domain) ||
         throw(ArgumentError("Invalid bipartitioned permutation"))
-    isperm((perm_codomain..., perm_domain...)) ||
+    isbiperm(perm_codomain, perm_domain) ||
         throw(ArgumentError("Invalid bipartitioned permutation"))
     return nothing
 end
diff --git a/src/matricize.jl b/src/matricize.jl
index b7e6711..4b88f0c 100644
--- a/src/matricize.jl
+++ b/src/matricize.jl
@@ -1,4 +1,3 @@
-using EllipsisNotation: Ellipsis
 using LinearAlgebra: Diagonal
 
 # =====================================  MatricizeStyle  ======================================
@@ -38,9 +37,7 @@ Non-mutating version of `bipermutedimsopadd!`: returns
 function permutedimsop(op, src, perm_codomain, perm_domain)
     # Validate against `src` here: `bipermutedimsopadd!`'s `check_input` compares against `dest`,
     # which `allocate_output` builds from the same perms, so it cannot catch a non-covering perm.
-    perm = (perm_codomain..., perm_domain...)
-    (ndims(src) == length(perm) && isperm(perm)) ||
-        throw(ArgumentError("Invalid bipermutation"))
+    check_biperm(src, perm_codomain, perm_domain)
     dest = allocate_output(permutedimsop, op, src, perm_codomain, perm_domain)
     return bipermutedimsopadd!(dest, op, src, perm_codomain, perm_domain, true, false)
 end
@@ -132,39 +129,6 @@ function matricizeperm(
     return matricizeopperm(style, identity, a, perm_codomain, perm_domain)
 end
 
-# Process inputs such as `EllipsisNotation.Ellipsis`.
-function to_permblocks(a, permblocks::NTuple{2, Tuple{Vararg{Int}}})
-    isperm((permblocks[1]..., permblocks[2]...)) ||
-        throw(ArgumentError("Invalid bipermutation"))
-    return permblocks
-end
-# Like `setcomplement` is like `setdiff` but assumes t2 ⊆ t1.
-function tuplesetcomplement(t1::NTuple{N1}, t2::NTuple{N2}) where {N1, N2}
-    t2 ⊆ t1 || throw(ArgumentError("t2 must be a subset of t1"))
-    return NTuple{N1 - N2}(setdiff(t1, t2))
-end
-function to_permblocks(
-        a, permblocks::Tuple{Tuple{Ellipsis}, Tuple{Vararg{Int}}}
-    )
-    permblocks1 = tuplesetcomplement(ntuple(identity, ndims(a)), permblocks[2])
-    return (permblocks1, permblocks[2])
-end
-function to_permblocks(
-        a, permblocks::Tuple{Tuple{Vararg{Int}}, Tuple{Ellipsis}}
-    )
-    permblocks2 = tuplesetcomplement(ntuple(identity, ndims(a)), permblocks[1])
-    return (permblocks[1], permblocks2)
-end
-
-function matricizeperm(a, perm_codomain, perm_domain)
-    return matricizeperm(MatricizeStyle(a), a, perm_codomain, perm_domain)
-end
-function matricizeperm(
-        style::MatricizeStyle, a, perm_codomain, perm_domain
-    )
-    return matricizeperm(style, a, to_permblocks(a, (perm_codomain, perm_domain))...)
-end
-
 # ==================================  matricizeopperm  =====================================
 
 """
@@ -177,13 +141,10 @@ Has "maybe alias" semantics: the result may be a view/wrapper aliasing `a` or a
 copy, depending on the matricize style and array type. The caller should treat the result
 as read-only.
 """
-function matricizeopperm(op, a, perm_codomain, perm_domain)
-    return matricizeopperm(MatricizeStyle(a), op, a, perm_codomain, perm_domain)
-end
 function matricizeopperm(
-        style::MatricizeStyle, op, a, perm_codomain, perm_domain
+        op, a, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
     )
-    return matricizeopperm(style, op, a, to_permblocks(a, (perm_codomain, perm_domain))...)
+    return matricizeopperm(MatricizeStyle(a), op, a, perm_codomain, perm_domain)
 end
 # Whether `perm` is the identity permutation `(1, …, n)`.
 isidentityperm(perm::Tuple{Vararg{Int}}) = perm == ntuple(identity, length(perm))
@@ -196,8 +157,7 @@ function matricizeopperm(
         style::MatricizeStyle, op, a,
         perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
     )
-    ndims(a) == length(perm_codomain) + length(perm_domain) ||
-        throw(ArgumentError("Invalid bipermutation"))
+    check_biperm(a, perm_codomain, perm_domain)
     op === identity && isidentityperm((perm_codomain..., perm_domain...)) &&
         return matricize(style, a, Val(length(perm_codomain)))
     a_perm_op = permutedimsop(op, a, perm_codomain, perm_domain)
diff --git a/test/Project.toml b/test/Project.toml
index 0ba4cc3..a2924b9 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -3,7 +3,6 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 ITensorPkgSkeleton = "3d388ab1-018a-49f4-ae50-18094d5f71ea"
 JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -27,7 +26,6 @@ Adapt = "4"
 Aqua = "0.8.9"
 BenchmarkTools = "1"
 BlockArrays = "1.6.1"
-EllipsisNotation = "1.8"
 ITensorPkgSkeleton = "0.3.42"
 JLArrays = "0.3"
 LinearAlgebra = "<0.0.1, 1"
diff --git a/test/test_basics.jl b/test/test_basics.jl
index 0d789da..218ce0f 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -1,5 +1,4 @@
 import TensorAlgebra
-using EllipsisNotation: var".."
 using StableRNGs: StableRNG
 using TensorAlgebra: BiTuple, ContractAlgorithm, bipermutedims, bipermutedims!, contract,
     contract!, contractadd!, length_codomain, length_domain, matricizeperm, unmatricize,

From 1ec949ffdeda1c67d4f4b61cb591609e6f7756e7 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 19:34:32 -0400
Subject: [PATCH 07/17] Leave the bipermutation arguments untyped

A hook dispatches on the array or the style, never on the permutation, so
annotating it narrowed what a backend may pass without buying any dispatch.
Base leaves `permutedims`' perm untyped and validates at runtime.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/factorizations.jl       | 68 ++++++++++++++++++-------------------
 src/matricize.jl            | 16 ++++-----
 src/matrixfunctions.jl      |  4 +--
 test/test_matricizestyle.jl |  2 +-
 4 files changed, 45 insertions(+), 45 deletions(-)

diff --git a/src/factorizations.jl b/src/factorizations.jl
index 1ab3963..5f7903e 100644
--- a/src/factorizations.jl
+++ b/src/factorizations.jl
@@ -29,7 +29,7 @@ for f in (
     @eval begin
         function $f(
                 style::MatricizeStyle, A,
-                perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}};
+                perm_codomain, perm_domain;
                 kwargs...
             )
             ndims(A) == length(perm_codomain) + length(perm_domain) ||
@@ -64,7 +64,7 @@ for f in (
     @eval begin
         function $f(
                 style::MatricizeStyle, A,
-                perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}};
+                perm_codomain, perm_domain;
                 kwargs...
             )
             A_mat = matricizeperm(style, A, perm_codomain, perm_domain)
@@ -102,7 +102,7 @@ for f in (
         end
         function $f(
                 A,
-                perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}};
+                perm_codomain, perm_domain;
                 kwargs...
             )
             return $f(MatricizeStyle(A), A, perm_codomain, perm_domain; kwargs...)
@@ -143,7 +143,7 @@ end
 
 """
     TensorAlgebra.tr(A, labels_A, labels_codomain, labels_domain)
-    TensorAlgebra.tr(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}})
+    TensorAlgebra.tr(A, perm_codomain, perm_domain)
     TensorAlgebra.tr(A, ndims_codomain::Val)
 
 Trace of a generic N-dimensional array `A` interpreted as a linear map from its domain to its
@@ -173,7 +173,7 @@ end
 function tr(A, ndims_codomain::Val)
     return tr(MatricizeStyle(A), A, ndims_codomain)
 end
-function tr(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}})
+function tr(A, perm_codomain, perm_domain)
     return LinearAlgebra.tr(matricizeperm(A, perm_codomain, perm_domain))
 end
 function tr(A, labels_A, labels_codomain, labels_domain)
@@ -184,7 +184,7 @@ end
 
 """
     qr_compact(A, labels_A, labels_codomain, labels_domain; kwargs...) -> Q, R
-    qr_compact(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> Q, R
+    qr_compact(A, perm_codomain, perm_domain; kwargs...) -> Q, R
     qr_compact(A, ndims_codomain::Val; kwargs...) -> Q, R
 
 Compute the compact QR decomposition of a generic N-dimensional array, by interpreting it
@@ -202,7 +202,7 @@ qr_compact
 
 """
     qr_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> Q, R
-    qr_full(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> Q, R
+    qr_full(A, perm_codomain, perm_domain; kwargs...) -> Q, R
     qr_full(A, ndims_codomain::Val; kwargs...) -> Q, R
 
 Compute the full QR decomposition of a generic N-dimensional array, by interpreting it as
@@ -220,7 +220,7 @@ qr_full
 
 """
     lq_compact(A, labels_A, labels_codomain, labels_domain; kwargs...) -> L, Q
-    lq_compact(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> L, Q
+    lq_compact(A, perm_codomain, perm_domain; kwargs...) -> L, Q
     lq_compact(A, ndims_codomain::Val; kwargs...) -> L, Q
 
 Compute the compact LQ decomposition of a generic N-dimensional array, by interpreting it
@@ -238,7 +238,7 @@ lq_compact
 
 """
     lq_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> L, Q
-    lq_full(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> L, Q
+    lq_full(A, perm_codomain, perm_domain; kwargs...) -> L, Q
     lq_full(A, ndims_codomain::Val; kwargs...) -> L, Q
 
 Compute the full LQ decomposition of a generic N-dimensional array, by interpreting it as
@@ -256,7 +256,7 @@ lq_full
 
 """
     left_polar(A, labels_A, labels_codomain, labels_domain; kwargs...) -> W, P
-    left_polar(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> W, P
+    left_polar(A, perm_codomain, perm_domain; kwargs...) -> W, P
     left_polar(A, ndims_codomain::Val; kwargs...) -> W, P
 
 Compute the left polar decomposition of a generic N-dimensional array, by interpreting it as
@@ -273,7 +273,7 @@ left_polar
 
 """
     right_polar(A, labels_A, labels_codomain, labels_domain; kwargs...) -> P, W
-    right_polar(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> P, W
+    right_polar(A, perm_codomain, perm_domain; kwargs...) -> P, W
     right_polar(A, ndims_codomain::Val; kwargs...) -> P, W
 
 Compute the right polar decomposition of a generic N-dimensional array, by interpreting it as
@@ -290,7 +290,7 @@ right_polar
 
 """
     left_orth(A, labels_A, labels_codomain, labels_domain; kwargs...) -> V, C
-    left_orth(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> V, C
+    left_orth(A, perm_codomain, perm_domain; kwargs...) -> V, C
     left_orth(A, ndims_codomain::Val; kwargs...) -> V, C
 
 Compute the left orthogonal decomposition of a generic N-dimensional array, by interpreting it as
@@ -307,7 +307,7 @@ left_orth
 
 """
     right_orth(A, labels_A, labels_codomain, labels_domain; kwargs...) -> C, V
-    right_orth(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> C, V
+    right_orth(A, perm_codomain, perm_domain; kwargs...) -> C, V
     right_orth(A, ndims_codomain::Val; kwargs...) -> C, V
 
 Compute the right orthogonal decomposition of a generic N-dimensional array, by interpreting it as
@@ -383,7 +383,7 @@ end
 
 """
     svd_compact(A, labels_A, labels_codomain, labels_domain; kwargs...) -> U, S, Vᴴ
-    svd_compact(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> U, S, Vᴴ
+    svd_compact(A, perm_codomain, perm_domain; kwargs...) -> U, S, Vᴴ
     svd_compact(A, ndims_codomain::Val; kwargs...) -> U, S, Vᴴ
 
 Compute the compact (thin) SVD of a generic N-dimensional array, by interpreting it as a
@@ -396,7 +396,7 @@ svd_compact
 
 """
     svd_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> U, S, Vᴴ
-    svd_full(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> U, S, Vᴴ
+    svd_full(A, perm_codomain, perm_domain; kwargs...) -> U, S, Vᴴ
     svd_full(A, ndims_codomain::Val; kwargs...) -> U, S, Vᴴ
 
 Compute the full (thick) SVD of a generic N-dimensional array, by interpreting it as a
@@ -409,7 +409,7 @@ svd_full
 
 """
     svd_trunc(A, labels_A, labels_codomain, labels_domain; trunc, kwargs...) -> U, S, Vᴴ, ϵ
-    svd_trunc(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; trunc, kwargs...) -> U, S, Vᴴ, ϵ
+    svd_trunc(A, perm_codomain, perm_domain; trunc, kwargs...) -> U, S, Vᴴ, ϵ
     svd_trunc(A, ndims_codomain::Val; trunc, kwargs...) -> U, S, Vᴴ, ϵ
 
 Compute the truncated SVD of a generic N-dimensional array, by interpreting it as a linear
@@ -446,7 +446,7 @@ svd_trunc
 
 """
     svd_vals(A, labels_A, labels_codomain, labels_domain) -> S
-    svd_vals(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}) -> S
+    svd_vals(A, perm_codomain, perm_domain) -> S
     svd_vals(A, ndims_codomain::Val) -> S
 
 Compute the singular values of a generic N-dimensional array, by interpreting it as a
@@ -459,7 +459,7 @@ svd_vals
 
 """
     eigh_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> D, V
-    eigh_full(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> D, V
+    eigh_full(A, perm_codomain, perm_domain; kwargs...) -> D, V
     eigh_full(A, ndims_codomain::Val; kwargs...) -> D, V
 
 Compute the eigenvalue decomposition of a generic N-dimensional array interpreted as a
@@ -472,7 +472,7 @@ eigh_full
 
 """
     eig_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> D, V
-    eig_full(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> D, V
+    eig_full(A, perm_codomain, perm_domain; kwargs...) -> D, V
     eig_full(A, ndims_codomain::Val; kwargs...) -> D, V
 
 Compute the eigenvalue decomposition of a generic N-dimensional array interpreted as a
@@ -486,7 +486,7 @@ eig_full
 
 """
     eigh_trunc(A, labels_A, labels_codomain, labels_domain; trunc, kwargs...) -> D, V
-    eigh_trunc(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; trunc, kwargs...) -> D, V
+    eigh_trunc(A, perm_codomain, perm_domain; trunc, kwargs...) -> D, V
     eigh_trunc(A, ndims_codomain::Val; trunc, kwargs...) -> D, V
 
 Truncated Hermitian eigenvalue decomposition, like [`eigh_full`](@ref) but keeping only the
@@ -498,7 +498,7 @@ eigh_trunc
 
 """
     eig_trunc(A, labels_A, labels_codomain, labels_domain; trunc, kwargs...) -> D, V
-    eig_trunc(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; trunc, kwargs...) -> D, V
+    eig_trunc(A, perm_codomain, perm_domain; trunc, kwargs...) -> D, V
     eig_trunc(A, ndims_codomain::Val; trunc, kwargs...) -> D, V
 
 Truncated general eigenvalue decomposition, like [`eig_full`](@ref) but keeping only the
@@ -510,7 +510,7 @@ eig_trunc
 
 """
     eigh_vals(A, labels_A, labels_codomain, labels_domain; kwargs...) -> D
-    eigh_vals(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> D
+    eigh_vals(A, perm_codomain, perm_domain; kwargs...) -> D
     eigh_vals(A, ndims_codomain::Val; kwargs...) -> D
 
 Compute the eigenvalues of a generic N-dimensional array interpreted as a Hermitian linear
@@ -522,7 +522,7 @@ eigh_vals
 
 """
     eig_vals(A, labels_A, labels_codomain, labels_domain; kwargs...) -> D
-    eig_vals(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> D
+    eig_vals(A, perm_codomain, perm_domain; kwargs...) -> D
     eig_vals(A, ndims_codomain::Val; kwargs...) -> D
 
 Compute the eigenvalues of a generic N-dimensional array interpreted as a general
@@ -535,7 +535,7 @@ eig_vals
 
 """
     left_null(A, labels_A, labels_codomain, labels_domain; kwargs...) -> N
-    left_null(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> N
+    left_null(A, perm_codomain, perm_domain; kwargs...) -> N
     left_null(A, ndims_codomain::Val; kwargs...) -> N
 
 Compute the left nullspace of a generic N-dimensional array, by interpreting it as
@@ -572,7 +572,7 @@ end
 
 """
     right_null(A, labels_A, labels_codomain, labels_domain; kwargs...) -> Nᴴ
-    right_null(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> Nᴴ
+    right_null(A, perm_codomain, perm_domain; kwargs...) -> Nᴴ
     right_null(A, ndims_codomain::Val::Val; kwargs...) -> Nᴴ
 
 Compute the right nullspace of a generic N-dimensional array, by interpreting it as
@@ -609,7 +609,7 @@ end
 
 """
     gram_eigh_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> X
-    gram_eigh_full(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> X
+    gram_eigh_full(A, perm_codomain, perm_domain; kwargs...) -> X
     gram_eigh_full(A, ndims_codomain::Val; kwargs...) -> X
 
 Gram factorization of a generic N-dimensional array, interpreting it as a
@@ -665,7 +665,7 @@ end
 
 """
     gram_eigh_full_with_pinv(A, labels_A, labels_codomain, labels_domain; kwargs...) -> X, Y
-    gram_eigh_full_with_pinv(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> X, Y
+    gram_eigh_full_with_pinv(A, perm_codomain, perm_domain; kwargs...) -> X, Y
     gram_eigh_full_with_pinv(A, ndims_codomain::Val; kwargs...) -> X, Y
 
 Like [`gram_eigh_full`](@ref), but additionally returns `Y ≈ pinv(X)` such
@@ -727,7 +727,7 @@ end
 
 """
     sqrth_safe(A, labels_A, labels_codomain, labels_domain; kwargs...) -> P
-    sqrth_safe(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> P
+    sqrth_safe(A, perm_codomain, perm_domain; kwargs...) -> P
     sqrth_safe(A, ndims_codomain::Val; kwargs...) -> P
 
 Square root of a generic N-dimensional array, interpreting it as a
@@ -750,7 +750,7 @@ sqrth_safe
 
 """
     invsqrth_safe(A, labels_A, labels_codomain, labels_domain; kwargs...) -> P
-    invsqrth_safe(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> P
+    invsqrth_safe(A, perm_codomain, perm_domain; kwargs...) -> P
     invsqrth_safe(A, ndims_codomain::Val; kwargs...) -> P
 
 Pseudo-inverse square root of a generic N-dimensional array, interpreting
@@ -784,7 +784,7 @@ end
 
 """
     project_hermitian(A, labels_A, labels_codomain, labels_domain; kwargs...) -> H
-    project_hermitian(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> H
+    project_hermitian(A, perm_codomain, perm_domain; kwargs...) -> H
     project_hermitian(A, ndims_codomain::Val; kwargs...) -> H
 
 Hermitian part `(M + M') / 2` of a generic N-dimensional array, interpreting
@@ -804,7 +804,7 @@ end
 
 """
     sqrth_invsqrth_safe(A, labels_A, labels_codomain, labels_domain; kwargs...) -> P, Pinv
-    sqrth_invsqrth_safe(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...) -> P, Pinv
+    sqrth_invsqrth_safe(A, perm_codomain, perm_domain; kwargs...) -> P, Pinv
     sqrth_invsqrth_safe(A, ndims_codomain::Val; kwargs...) -> P, Pinv
 
 Square root and pseudo-inverse square root of a generic N-dimensional
@@ -833,7 +833,7 @@ end
 
 """
     TensorAlgebra.one(A, labels_A, labels_codomain, labels_domain) -> Id
-    TensorAlgebra.one(A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}) -> Id
+    TensorAlgebra.one(A, perm_codomain, perm_domain) -> Id
     TensorAlgebra.one(A, ndims_codomain::Val) -> Id
 
 Construct the identity operator tensor whose shape mirrors `A`, interpreted as a
@@ -901,14 +901,14 @@ end
 # `TensorMap` backend fills the identity through TensorKit rather than MatrixAlgebraKit).
 function one(
         style::MatricizeStyle, A,
-        perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}};
+        perm_codomain, perm_domain;
         kwargs...
     )
     A_perm = bipermutedims(A, perm_codomain, perm_domain)
     return one!!(style, A_perm, Val(length(perm_codomain)); kwargs...)
 end
 function one(
-        A, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}; kwargs...
+        A, perm_codomain, perm_domain; kwargs...
     )
     return one(MatricizeStyle(A), A, perm_codomain, perm_domain; kwargs...)
 end
diff --git a/src/matricize.jl b/src/matricize.jl
index 4b88f0c..ac146d9 100644
--- a/src/matricize.jl
+++ b/src/matricize.jl
@@ -107,7 +107,7 @@ end
 # guaranteed to be a copy.
 function matricizecopy(
         style::MatricizeStyle, a,
-        perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
+        perm_codomain, perm_domain
     )
     a_perm = bipermutedims(a, perm_codomain, perm_domain)
     return matricize(style, a_perm, Val(length(perm_codomain)))
@@ -115,7 +115,7 @@ end
 
 function matricizeperm(
         a,
-        perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
+        perm_codomain, perm_domain
     )
     return matricizeperm(MatricizeStyle(a), a, perm_codomain, perm_domain)
 end
@@ -124,7 +124,7 @@ end
 # `matricizeopperm`.
 function matricizeperm(
         style::MatricizeStyle, a,
-        perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
+        perm_codomain, perm_domain
     )
     return matricizeopperm(style, identity, a, perm_codomain, perm_domain)
 end
@@ -142,7 +142,7 @@ copy, depending on the matricize style and array type. The caller should treat t
 as read-only.
 """
 function matricizeopperm(
-        op, a, perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
+        op, a, perm_codomain, perm_domain
     )
     return matricizeopperm(MatricizeStyle(a), op, a, perm_codomain, perm_domain)
 end
@@ -155,7 +155,7 @@ isidentityperm(perm::Tuple{Vararg{Int}}) = perm == ntuple(identity, length(perm)
 # alias `a` and must be treated as read-only, matching the docstring.
 function matricizeopperm(
         style::MatricizeStyle, op, a,
-        perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
+        perm_codomain, perm_domain
     )
     check_biperm(a, perm_codomain, perm_domain)
     op === identity && isidentityperm((perm_codomain..., perm_domain...)) &&
@@ -175,7 +175,7 @@ end
 ismatricizeview(::MatricizeStyle, a, ndims_codomain::Val) = false
 function ismatricizeview(
         style::MatricizeStyle, a,
-        perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
+        perm_codomain, perm_domain
     )
     isidentityperm((perm_codomain..., perm_domain...)) || return false
     return ismatricizeview(style, a, Val(length(perm_codomain)))
@@ -213,13 +213,13 @@ end
 # the same forward bipermutation to both.
 function unmatricize!(
         a_dest, m,
-        perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
+        perm_codomain, perm_domain
     )
     return unmatricize!(MatricizeStyle(m), a_dest, m, perm_codomain, perm_domain)
 end
 function unmatricize!(
         style::MatricizeStyle, a_dest, m,
-        perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
+        perm_codomain, perm_domain
     )
     biperm_src = BiTuple(perm_codomain, perm_domain)
     ndims(a_dest) == length(biperm_src) ||
diff --git a/src/matrixfunctions.jl b/src/matrixfunctions.jl
index ebbda11..1747244 100644
--- a/src/matrixfunctions.jl
+++ b/src/matrixfunctions.jl
@@ -50,7 +50,7 @@ for f in MATRIX_FUNCTIONS
 
         function $f(
                 style::MatricizeStyle, a,
-                perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}};
+                perm_codomain, perm_domain;
                 kwargs...
             )
             a_mat = matricizeperm(style, a, perm_codomain, perm_domain)
@@ -63,7 +63,7 @@ for f in MATRIX_FUNCTIONS
         end
         function $f(
                 a,
-                perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}};
+                perm_codomain, perm_domain;
                 kwargs...
             )
             return $f(MatricizeStyle(a), a, perm_codomain, perm_domain; kwargs...)
diff --git a/test/test_matricizestyle.jl b/test/test_matricizestyle.jl
index 038ac12..8d13c40 100644
--- a/test/test_matricizestyle.jl
+++ b/test/test_matricizestyle.jl
@@ -18,7 +18,7 @@ module MatricizeStyleTestUtils
     end
     function TA.unmatricize!(
             ::MyArrayMatricize, a_dest::MyArray, m,
-            perm_codomain::Tuple{Vararg{Int}}, perm_domain::Tuple{Vararg{Int}}
+            perm_codomain, perm_domain
         )
         TA.unmatricize!(
             TA.ReshapeMatricize(), a_dest.parent, m, perm_codomain, perm_domain

From f5d7f6ad4eefd6a4463b7f8eb67299676d7b608a Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 19:36:35 -0400
Subject: [PATCH 08/17] Ask about identity bipermutations with both halves

Both callers had a bipermutation in hand and were splatting it to reach
`isidentityperm`, the same shape problem `isbiperm` fixed. Also drops the
Ellipsis spellings from the matricize tests, which the removed normalizing
tier supported.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/bituple.jl      |  8 ++++++++
 src/matricize.jl    |  6 ++----
 test/test_basics.jl | 10 ++++------
 3 files changed, 14 insertions(+), 10 deletions(-)

diff --git a/src/bituple.jl b/src/bituple.jl
index 0501193..f898ee7 100644
--- a/src/bituple.jl
+++ b/src/bituple.jl
@@ -75,3 +75,11 @@ bipartition(t::Tuple, bt::BiTuple) = bipartition(t, bt.t1, bt.t2)
 # bipartitioned permutation. The two halves only make sense jointly, so this takes them as a pair
 # rather than leaving every caller to splat and call `isperm`.
 isbiperm(perm_codomain, perm_domain) = isperm((perm_codomain..., perm_domain...))
+
+# Whether `perm_codomain` and `perm_domain` are the identity bipermutation, i.e. leave every
+# dimension where it is. Takes the halves for the same reason `isbiperm` does: both call sites had
+# a bipermutation in hand and were splatting it to ask.
+function isidentitybiperm(perm_codomain, perm_domain)
+    perm = (perm_codomain..., perm_domain...)
+    return perm == ntuple(identity, length(perm))
+end
diff --git a/src/matricize.jl b/src/matricize.jl
index ac146d9..70b2beb 100644
--- a/src/matricize.jl
+++ b/src/matricize.jl
@@ -146,8 +146,6 @@ function matricizeopperm(
     )
     return matricizeopperm(MatricizeStyle(a), op, a, perm_codomain, perm_domain)
 end
-# Whether `perm` is the identity permutation `(1, …, n)`.
-isidentityperm(perm::Tuple{Vararg{Int}}) = perm == ntuple(identity, length(perm))
 
 # The identity bipermutation is a no-op permute, so `matricize` runs directly on `a` (a view
 # for dense, a gather without the extra permute copy for graded); the fast path requires
@@ -158,7 +156,7 @@ function matricizeopperm(
         perm_codomain, perm_domain
     )
     check_biperm(a, perm_codomain, perm_domain)
-    op === identity && isidentityperm((perm_codomain..., perm_domain...)) &&
+    op === identity && isidentitybiperm(perm_codomain, perm_domain) &&
         return matricize(style, a, Val(length(perm_codomain)))
     a_perm_op = permutedimsop(op, a, perm_codomain, perm_domain)
     return matricize(style, a_perm_op, Val(length(perm_codomain)))
@@ -177,7 +175,7 @@ function ismatricizeview(
         style::MatricizeStyle, a,
         perm_codomain, perm_domain
     )
-    isidentityperm((perm_codomain..., perm_domain...)) || return false
+    isidentitybiperm(perm_codomain, perm_domain) || return false
     return ismatricizeview(style, a, Val(length(perm_codomain)))
 end
 
diff --git a/test/test_basics.jl b/test/test_basics.jl
index 218ce0f..8f45811 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -62,17 +62,15 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         a_fused = matricizeperm(a, (3, 1, 2), (4,))
         @test eltype(a_fused) === elt
         @test a_fused ≈ reshape(permutedims(a, (3, 1, 2, 4)), (24, 5))
-        a_fused = matricizeperm(a, (..,), (3, 1))
+        a_fused = matricizeperm(a, (2, 4), (3, 1))
         @test eltype(a_fused) === elt
         @test a_fused ≈ reshape(permutedims(a, (2, 4, 3, 1)), (15, 8))
-        a_fused = matricizeperm(a, (3, 1), (..,))
-        @test eltype(a_fused) === elt
-        @test a_fused ≈ reshape(permutedims(a, (3, 1, 2, 4)), (8, 15))
 
-        a_fused = matricizeperm(a, (), (..,))
+        # Degenerate splits: everything in the domain, then everything in the codomain.
+        a_fused = matricizeperm(a, (), (1, 2, 3, 4))
         @test eltype(a_fused) === elt
         @test a_fused ≈ reshape(a, (1, 120))
-        a_fused = matricizeperm(a, (..,), ())
+        a_fused = matricizeperm(a, (1, 2, 3, 4), ())
         @test eltype(a_fused) === elt
         @test a_fused ≈ reshape(a, (120, 1))
 

From 77cc3c7a7f2fed42757b7deb3d5081b0ce162220 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 20:50:58 -0400
Subject: [PATCH 09/17] Give matricize one primitive per operation

A style now implements four bipermutation hooks, `allocate_output`,
`matricizeop!`, `matricizeopview` and `is_output_view`, and the copy and
maybe-alias forms are derived. Allocation is a hook because only the style
knows its fused axes, and because it is what makes the copy path terminate.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/TensorAlgebra.jl               |   2 +-
 src/contract/contract_matricize.jl |  16 ++-
 src/diagonal.jl                    |   6 +-
 src/factorizations.jl              |  32 +++--
 src/matricize.jl                   | 205 +++++++++++++++--------------
 src/matrixfunctions.jl             |   4 +-
 6 files changed, 147 insertions(+), 118 deletions(-)

diff --git a/src/TensorAlgebra.jl b/src/TensorAlgebra.jl
index 2853d1f..0ac1bc9 100644
--- a/src/TensorAlgebra.jl
+++ b/src/TensorAlgebra.jl
@@ -9,7 +9,7 @@ export contract, contract!, dual, eig_full, eig_trunc, eig_vals, eigh_full, eigh
 if VERSION >= v"1.11.0-DEV.469"
     eval(
         Meta.parse(
-            "public biperm, bipartition, cat_similar, concatenate, concatenate!, ContractAlgorithm, contractopadd!, data, datatype, directsum, flattenlinear, label_type, matricizeopperm, permutedims, permutedims!, scalar, similar_map, TensorOperationsAlgorithm, to_range, tr, tryflattenlinear, ungrade, zero!, scale!, permuteddims, PermutedDims"
+            "public biperm, bipartition, cat_similar, concatenate, concatenate!, ContractAlgorithm, contractopadd!, data, datatype, directsum, flattenlinear, is_output_view, label_type, matricize, matricizeop, matricizeop!, matricizeopcopy, matricizeopview, permutedims, permutedims!, scalar, similar_map, TensorOperationsAlgorithm, to_range, tr, tryflattenlinear, ungrade, zero!, scale!, permuteddims, PermutedDims"
         )
     )
 end
diff --git a/src/contract/contract_matricize.jl b/src/contract/contract_matricize.jl
index 1593d73..0b5fc95 100644
--- a/src/contract/contract_matricize.jl
+++ b/src/contract/contract_matricize.jl
@@ -16,16 +16,20 @@ function contractopadd!(
         a1, biperm1_codomain, biperm1_domain,
         a2, biperm2_codomain, biperm2_domain
     )
-    a1_mat = matricizeopperm(
+    a1_mat = matricizeop(
         algorithm.left_matricize_style, op1, a1, biperm1_codomain, biperm1_domain
     )
-    a2_mat = matricizeopperm(
+    a2_mat = matricizeop(
         algorithm.right_matricize_style, op2, a2, biperm2_codomain, biperm2_domain
     )
     output_style = algorithm.output_matricize_style
-    if ismatricizeview(output_style, a_dest, invperm_codomain, invperm_domain)
+    if is_output_view(
+            matricizeop, output_style, identity, a_dest, invperm_codomain, invperm_domain
+        )
         # The matricization shares `a_dest`'s memory, so the matmul is the whole operation.
-        a_dest_mat = matricizeview(output_style, a_dest, Val(length(invperm_codomain)))
+        a_dest_mat = matricizeopview(
+            output_style, identity, a_dest, invperm_codomain, invperm_domain
+        )
         mul!(a_dest_mat, a1_mat, a2_mat, α, β)
     elseif iszero(β)
         # `β` is a strong zero, so `a_dest`'s current data is irrelevant: let the matmul
@@ -38,7 +42,9 @@ function contractopadd!(
     else
         # `a_dest`'s data contributes through `β`, so gather it, multiply into the gathered
         # copy, and scatter back.
-        a_dest_mat = matricizecopy(output_style, a_dest, invperm_codomain, invperm_domain)
+        a_dest_mat = matricizeopcopy(
+            output_style, identity, a_dest, invperm_codomain, invperm_domain
+        )
         mul!(a_dest_mat, a1_mat, a2_mat, α, β)
         unmatricize!(output_style, a_dest, a_dest_mat, invperm_codomain, invperm_domain)
     end
diff --git a/src/diagonal.jl b/src/diagonal.jl
index 7a46589..0488eff 100644
--- a/src/diagonal.jl
+++ b/src/diagonal.jl
@@ -44,7 +44,11 @@ end
 # A `Diagonal` is already a matrix; the `(1 codomain, 1 domain)` matricization is the identity
 # reshape, so the memory-sharing matricization is `a` itself (keeping it a `Diagonal` for the
 # `Diagonal`-specialized consumers downstream).
-matricizeview(::ReshapeMatricize, a::Diagonal, ::Val{1}) = a
+function matricizeopview(
+        ::ReshapeMatricize, op, a::Diagonal, perm_codomain::Tuple{Int}, perm_domain::Tuple{Int}
+    )
+    return a
+end
 # A `{1,1}` unmatricize (one codomain axis, one domain axis) is the endomorphism identity: the
 # result stays `Diagonal`, so return `m` directly. The generic `check_input(unmatricize, ...)`
 # validates the axis lengths against `m`'s size.
diff --git a/src/factorizations.jl b/src/factorizations.jl
index 5f7903e..ad1a1ca 100644
--- a/src/factorizations.jl
+++ b/src/factorizations.jl
@@ -16,7 +16,7 @@ using MatrixAlgebraKit: MatrixAlgebraKit
 # Owned tier: the matrix-level entries mutate their input, so the perm form materializes an
 # owned matricization following MatrixAlgebraKit's `f(A) = f!(copy_input(f, A))` convention —
 # a memory-sharing matricization is materialized through `MatrixAlgebraKit.copy_input`, while
-# the `matricizecopy` gather is owned by contract and is donated directly (with
+# the `matricizeopcopy` gather is owned by contract and is donated directly (with
 # `copy_input` still applied when the eltype must change) — and the wrapper calls the mutating
 # entry unconditionally.
 for f in (
@@ -34,11 +34,16 @@ for f in (
             )
             ndims(A) == length(perm_codomain) + length(perm_domain) ||
                 throw(ArgumentError("Invalid bipermutation"))
-            A_mat = if ismatricizeview(style, A, perm_codomain, perm_domain)
-                A_shared = matricizeview(style, A, Val(length(perm_codomain)))
+            A_mat =
+            if is_output_view(
+                    matricizeop, style, identity, A, perm_codomain, perm_domain
+                )
+                A_shared =
+                    matricizeopview(style, identity, A, perm_codomain, perm_domain)
                 MatrixAlgebraKit.copy_input(MatrixAlgebraKit.$f, A_shared)
             else
-                A_gather = matricizecopy(style, A, perm_codomain, perm_domain)
+                A_gather =
+                    matricizeopcopy(style, identity, A, perm_codomain, perm_domain)
                 if eltype(A_gather) === float(eltype(A_gather))
                     A_gather
                 else
@@ -56,7 +61,7 @@ for f in (
 end
 
 # Read-only tier: the matrix-level entries never mutate their input (they copy internally), so
-# the perm form consumes the maybe-alias `matricizeperm` matricization directly.
+# the perm form consumes the maybe-alias `matricize` matricization directly.
 for f in (
         :gram_eigh_full, :gram_eigh_full_with_pinv,
         :sqrth_safe, :invsqrth_safe, :sqrth_invsqrth_safe,
@@ -67,7 +72,7 @@ for f in (
                 perm_codomain, perm_domain;
                 kwargs...
             )
-            A_mat = matricizeperm(style, A, perm_codomain, perm_domain)
+            A_mat = matricize(style, A, perm_codomain, perm_domain)
             F = MatrixAlgebra.$f(A_mat; kwargs...)
             axes_codomain, axes_domain = bipartition_axes(
                 map(i -> axes(A, i), (perm_codomain..., perm_domain...)),
@@ -174,7 +179,7 @@ function tr(A, ndims_codomain::Val)
     return tr(MatricizeStyle(A), A, ndims_codomain)
 end
 function tr(A, perm_codomain, perm_domain)
-    return LinearAlgebra.tr(matricizeperm(A, perm_codomain, perm_domain))
+    return LinearAlgebra.tr(matricize(A, perm_codomain, perm_domain))
 end
 function tr(A, labels_A, labels_codomain, labels_domain)
     perm_codomain, perm_domain =
@@ -858,14 +863,14 @@ julia> A = randn(2, 3, 2, 3);
 
 julia> Id = TensorAlgebra.one(A, (:a, :b, :c, :d), (:a, :b), (:c, :d));
 
-julia> matricize(Id, Val(2)) ≈ I
+julia> matricize(Id, (1, 2), (3, 4)) ≈ I
 true
 ```
 """
 function one end
 
 function one!!(style::MatricizeStyle, A, ndims_codomain::Val; kwargs...)
-    A_mat = matricize(style, A, ndims_codomain)
+    A_mat = matricize(style, A, trivialbiperm(A, ndims_codomain)...)
     MatrixAlgebraKit.one!(A_mat)
     axes_codomain, axes_domain = bipartition_axes(axes(A), ndims_codomain)
     return unmatricize(style, A_mat, axes_codomain, axes_domain)
@@ -878,11 +883,14 @@ end
 # matricization directly when the style declares one at this split, and otherwise fills a
 # gathered matrix and scatters it back with `unmatricize!`.
 function one!(style::MatricizeStyle, A, ndims_codomain::Val; kwargs...)
-    if ismatricizeview(style, A, ndims_codomain)
-        MatrixAlgebraKit.one!(matricizeview(style, A, ndims_codomain))
+    perm_codomain, perm_domain = trivialbiperm(A, ndims_codomain)
+    if is_output_view(matricizeop, style, identity, A, perm_codomain, perm_domain)
+        MatrixAlgebraKit.one!(
+            matricizeopview(style, identity, A, perm_codomain, perm_domain)
+        )
         return A
     end
-    A_mat = matricizecopy(style, A, ndims_codomain)
+    A_mat = matricizeopcopy(style, identity, A, perm_codomain, perm_domain)
     MatrixAlgebraKit.one!(A_mat)
     return unmatricize!(style, A, A_mat, ndims_codomain)
 end
diff --git a/src/matricize.jl b/src/matricize.jl
index 70b2beb..aa6b3fe 100644
--- a/src/matricize.jl
+++ b/src/matricize.jl
@@ -70,116 +70,111 @@ function bipermutedims!(
 end
 
 # =====================================  matricize  ========================================
-# Copy convention: `bipermutedims`/`permutedims` always copy (Base `permutedims` semantics). At
-# the trivial (`Val`) split the sharing story is exact: `matricizeview` shares `a`'s memory,
-# `matricizecopy` returns fresh storage the caller owns, and `matricize` aliases `a` iff
-# `ismatricizeview` — so a consumer that writes into a destination checks the trait and writes
-# through `matricizeview`, and a consumer that mutates an input either owns a `matricizecopy`
-# result by contract or materializes an owned matrix with `MatrixAlgebraKit.copy_input` (see the
-# owned tier in `factorizations.jl`). `matricizeperm`/`matricizeopperm` keep maybe-alias
-# semantics (the result may view or copy; treat it as read-only) until the planned op/perm-form
-# trait lands, and `matricizeview` deliberately has no perm form pending that op/perm-layer
-# design.
+# A style implements four hooks, all taking the operation, the array and the bipermutation:
+#
+#   `allocate_output(matricizeop, style, op, a, pc, pd)`  the matrix destination
+#   `matricizeop!(dest, style, op, a, pc, pd)`            write the matricization into it
+#   `matricizeopview(style, op, a, pc, pd)`               partial: the aliasing form
+#   `is_output_view(matricizeop, style, op, a, pc, pd)`   whether the aliasing form applies
+#
+# Everything else is derived. `matricizeopcopy` allocates and writes, so it always returns fresh
+# storage the caller owns. `matricizeop` returns the view where the style declares one and the copy
+# otherwise, i.e. it has maybe-alias semantics and its result must be treated as read-only.
+# `matricize` is `matricizeop` at `identity`.
+#
+# Allocation is a hook rather than generic machinery because computing a matricized destination
+# needs the fused axes, which only the style knows: TensorAlgebra deliberately has no generic
+# axis-fusion interface. It is also what makes the copy path terminate, since `matricizeop!` is a
+# distinct function from the router rather than a re-entry into it.
 
-# `matricize` at the trivial split routes on the style's sharing declaration. Styles implement
-# the three leaves (`ismatricizeview`, `matricizeview`, `matricizecopy`) rather than overloading
-# `matricize` itself. This assumes the permutation was already performed.
-function matricize(style::MatricizeStyle, a, ndims_codomain::Val)
-    ismatricizeview(style, a, ndims_codomain) &&
-        return matricizeview(style, a, ndims_codomain)
-    return matricizecopy(style, a, ndims_codomain)
+"""
+    matricizeop(op, a, perm_codomain, perm_domain)
+
+Matricize `a` across the bipermutation with the element-wise operation `op` folded in, i.e. a
+matrix representing `op.(permutedims(a, (perm_codomain..., perm_domain...)))` with the codomain
+fused to rows and the domain to columns.
+
+Has "maybe alias" semantics: the result may share `a`'s memory or be fresh storage, depending on
+the style and the array type. Treat it as read-only. Use `matricizeopcopy` for a matrix the caller
+owns, and `matricizeopview` (partial) for one guaranteed to alias.
+"""
+function matricizeop(op, a, perm_codomain, perm_domain)
+    return matricizeop(MatricizeStyle(a), op, a, perm_codomain, perm_domain)
 end
-function matricize(a, ndims_codomain::Val)
-    return matricize(MatricizeStyle(a), a, ndims_codomain)
+function matricizeop(style::MatricizeStyle, op, a, perm_codomain, perm_domain)
+    check_biperm(a, perm_codomain, perm_domain)
+    is_output_view(matricizeop, style, op, a, perm_codomain, perm_domain) &&
+        return matricizeopview(style, op, a, perm_codomain, perm_domain)
+    return matricizeopcopy(style, op, a, perm_codomain, perm_domain)
 end
 
-# Partial: defined only where `ismatricizeview` is `true`, and always returns a matricization
-# sharing `a`'s memory (the `StridedView` partial-constructor pattern).
-function matricizeview(style::MatricizeStyle, a, ndims_codomain::Val)
-    return throw(MethodError(matricizeview, (style, a, ndims_codomain)))
+"""
+    matricize(a, perm_codomain, perm_domain)
+
+`matricizeop` at `identity`. Has the same maybe-alias semantics.
+"""
+function matricize(a, perm_codomain, perm_domain)
+    return matricizeop(identity, a, perm_codomain, perm_domain)
 end
-# Total: always returns a matricization in fresh storage the caller owns.
-function matricizecopy(style::MatricizeStyle, a, ndims_codomain::Val)
-    return throw(MethodError(matricizecopy, (style, a, ndims_codomain)))
+function matricize(style::MatricizeStyle, a, perm_codomain, perm_domain)
+    return matricizeop(style, identity, a, perm_codomain, perm_domain)
 end
 
-# `bipermutedims` always copies and `matricize` might return a view, so the result is
-# guaranteed to be a copy.
-function matricizecopy(
-        style::MatricizeStyle, a,
-        perm_codomain, perm_domain
-    )
-    a_perm = bipermutedims(a, perm_codomain, perm_domain)
-    return matricize(style, a_perm, Val(length(perm_codomain)))
+# Total: always fresh storage the caller owns.
+function matricizeopcopy(op, a, perm_codomain, perm_domain)
+    return matricizeopcopy(MatricizeStyle(a), op, a, perm_codomain, perm_domain)
+end
+function matricizeopcopy(style::MatricizeStyle, op, a, perm_codomain, perm_domain)
+    check_biperm(a, perm_codomain, perm_domain)
+    dest = allocate_output(matricizeop, style, op, a, perm_codomain, perm_domain)
+    return matricizeop!(dest, style, op, a, perm_codomain, perm_domain)
 end
 
-function matricizeperm(
-        a,
-        perm_codomain, perm_domain
+# Partial: defined only where `is_output_view` is `true`, and always returns a matricization
+# sharing `a`'s memory (the `StridedView` partial-constructor pattern).
+function matricizeopview(style::MatricizeStyle, op, a, perm_codomain, perm_domain)
+    return throw(
+        MethodError(matricizeopview, (style, op, a, perm_codomain, perm_domain))
     )
-    return matricizeperm(MatricizeStyle(a), a, perm_codomain, perm_domain)
 end
-# Thin wrapper around `matricizeopperm` with identity op — the actual matricization logic
-# (and the matricize-style overload point for folding ops into matricization) lives in
-# `matricizeopperm`.
-function matricizeperm(
-        style::MatricizeStyle, a,
-        perm_codomain, perm_domain
+
+# Required of every style: write the matricization of `a` into `dest`.
+function matricizeop!(dest, style::MatricizeStyle, op, a, perm_codomain, perm_domain)
+    return throw(
+        MethodError(matricizeop!, (dest, style, op, a, perm_codomain, perm_domain))
     )
-    return matricizeopperm(style, identity, a, perm_codomain, perm_domain)
 end
 
-# ==================================  matricizeopperm  =====================================
-
-"""
-    matricizeopperm(op, a, perm_codomain, perm_domain)
-
-Matricize `a` with element-wise operation `op` folded in. Returns a matrix representing
-`op.(matricizeperm(a, perm_codomain, perm_domain))`.
-
-Has "maybe alias" semantics: the result may be a view/wrapper aliasing `a` or a fresh
-copy, depending on the matricize style and array type. The caller should treat the result
-as read-only.
-"""
-function matricizeopperm(
-        op, a, perm_codomain, perm_domain
+# Required of every style: the matrix destination `matricizeop!` writes into.
+function allocate_output(
+        ::typeof(matricizeop), style::MatricizeStyle, op, a, perm_codomain, perm_domain
+    )
+    return throw(
+        MethodError(
+            allocate_output,
+            (matricizeop, style, op, a, perm_codomain, perm_domain)
+        )
     )
-    return matricizeopperm(MatricizeStyle(a), op, a, perm_codomain, perm_domain)
 end
 
-# The identity bipermutation is a no-op permute, so `matricize` runs directly on `a` (a view
-# for dense, a gather without the extra permute copy for graded); the fast path requires
-# `op === identity`, since a plain view cannot carry a fused `op` like `conj`. The result may
-# alias `a` and must be treated as read-only, matching the docstring.
-function matricizeopperm(
-        style::MatricizeStyle, op, a,
-        perm_codomain, perm_domain
-    )
-    check_biperm(a, perm_codomain, perm_domain)
-    op === identity && isidentitybiperm(perm_codomain, perm_domain) &&
-        return matricize(style, a, Val(length(perm_codomain)))
-    a_perm_op = permutedimsop(op, a, perm_codomain, perm_domain)
-    return matricize(style, a_perm_op, Val(length(perm_codomain)))
+# The trivial bipermutation for a rank-`N` array split after `ndims_codomain` dimensions. The
+# `Val` entry points that remain build it to reach the bipermutation hooks.
+function trivialbiperm(a, ndims_codomain::Val{K}) where {K}
+    return ntuple(identity, Val(K)), ntuple(i -> K + i, Val(ndims(a) - K))
 end
 
-# ==================================  ismatricizeview  =====================================
-# `true` iff `matricize(style, a, ndims_codomain)` shares `a`'s memory (writes to it are writes
-# to `a`) — the `isstrided`/`StridedView` pattern (also TensorKit's `has_shared_permute` and
-# TensorOperations' `isblasdestination`). Styles overload the `Val` (trivial split) form to
-# declare which splits share; the bipermutation form delegates to it at the identity and is
-# `false` (fail-safe) everywhere else. A general `ismatricizeview(style, op, a, perm_codomain,
-# perm_domain)` form (op and bipermutation view-sets) is planned; these are its
-# `op === identity` special cases.
-ismatricizeview(::MatricizeStyle, a, ndims_codomain::Val) = false
-function ismatricizeview(
-        style::MatricizeStyle, a,
-        perm_codomain, perm_domain
+# ==================================  is_output_view  ======================================
+# `true` iff `matricizeop(style, op, a, perm_codomain, perm_domain)` shares `a`'s memory, so that
+# writes to the result are writes to `a`. The `isstrided`/`StridedView` pattern, and the same
+# question TensorKit asks with `has_shared_permute` and TensorOperations with `isblasdestination`.
+# Keyed on the operation like the other function-keyed hooks (`check_input`, `allocate_output`,
+# `output_axes`), so the predicate's arguments are exactly the call's arguments.
+function is_output_view(
+        ::typeof(matricizeop), ::MatricizeStyle, op, a, perm_codomain, perm_domain
     )
-    isidentitybiperm(perm_codomain, perm_domain) || return false
-    return ismatricizeview(style, a, Val(length(perm_codomain)))
+    return false
 end
 
-# ====================================  unmatricize  =======================================
 # Split form: `axes_codomain` and `axes_domain` are the destination axes for the codomain and
 # domain groups, given codomain-facing (un-dualized), the same convention as `similar_map`. A
 # matricize style stores the domain axes dualized, so its overload re-dualizes them with `conj`
@@ -207,7 +202,7 @@ end
 # The bipermutation maps the destination's dimension order to the matrix's: `axes(a_dest)` grouped
 # by it gives the legs in `m`'s order, and the result is permuted back by its inverse. It is not
 # intrinsically an inverse permutation — the matricized-contraction destination path happens to
-# derive it as `invperm(biperm_dest)`, while a `matricizeperm`/`unmatricize!` round trip passes
+# derive it as `invperm(biperm_dest)`, while a `matricize`/`unmatricize!` round trip passes
 # the same forward bipermutation to both.
 function unmatricize!(
         a_dest, m,
@@ -250,16 +245,32 @@ end
 # Defaults to ReshapeMatricize, a simple reshape
 struct ReshapeMatricize <: MatricizeStyle end
 MatricizeStyle(::Type{<:AbstractArray}) = ReshapeMatricize()
-# A dense reshape matricization is a lazy wrapper at any split, so it always shares memory.
-ismatricizeview(::ReshapeMatricize, a, ndims_codomain::Val) = true
-function matricizeview(::ReshapeMatricize, a, ndims_codomain::Val)
-    unval(ndims_codomain) ≤ ndims(a) ||
-        throw(ArgumentError("Codomain length exceeds number of dimensions."))
-    size_codomain, size_domain = bipartition(size(a), ndims_codomain)
+# A dense reshape shares memory only when the data is already in codomain-then-domain order and
+# no operation has to be folded in: a reshape can neither reorder nor carry a `conj`.
+function is_output_view(
+        ::typeof(matricizeop), ::ReshapeMatricize, op, a, perm_codomain, perm_domain
+    )
+    return op === identity && isidentitybiperm(perm_codomain, perm_domain)
+end
+function matricizeopview(::ReshapeMatricize, op, a, perm_codomain, perm_domain)
+    size_codomain, size_domain = bipartition(size(a), Val(length(perm_codomain)))
     return reshape(a, (prod(size_codomain), prod(size_domain)))
 end
-function matricizecopy(style::ReshapeMatricize, a, ndims_codomain::Val)
-    return copy(matricizeview(style, a, ndims_codomain))
+function allocate_output(
+        ::typeof(matricizeop), ::ReshapeMatricize, op, a, perm_codomain, perm_domain
+    )
+    T = Base.promote_op(op, eltype(a))
+    size_codomain = map(i -> size(a, i), perm_codomain)
+    size_domain = map(i -> size(a, i), perm_domain)
+    return similar(a, T, (prod(size_codomain), prod(size_domain)))
+end
+# The destination is a dense matrix, so reshaping it to the permuted tensor shape is a view and
+# the permuted-add writes straight through it.
+function matricizeop!(dest, ::ReshapeMatricize, op, a, perm_codomain, perm_domain)
+    perm = (perm_codomain..., perm_domain...)
+    dest_tensor = reshape(dest, map(i -> size(a, i), perm))
+    bipermutedimsopadd!(dest_tensor, op, a, perm_codomain, perm_domain, true, false)
+    return dest
 end
 # The matricized input's rows must be the fused codomain and its columns the fused domain.
 # `reshape` alone only checks the total element count, so a wrong split with the right total
diff --git a/src/matrixfunctions.jl b/src/matrixfunctions.jl
index 1747244..4b5f524 100644
--- a/src/matrixfunctions.jl
+++ b/src/matrixfunctions.jl
@@ -32,7 +32,7 @@ const MATRIX_FUNCTIONS = [
 ]
 
 # The matrix functions never mutate their input (they allocate their own outputs), so the
-# permuted forms consume the maybe-alias `matricizeperm` matricization read-only, skipping
+# permuted forms consume the maybe-alias `matricize` matricization read-only, skipping
 # the eager `bipermutedims` copy at the identity bipermutation.
 for f in MATRIX_FUNCTIONS
     @eval begin
@@ -53,7 +53,7 @@ for f in MATRIX_FUNCTIONS
                 perm_codomain, perm_domain;
                 kwargs...
             )
-            a_mat = matricizeperm(style, a, perm_codomain, perm_domain)
+            a_mat = matricize(style, a, perm_codomain, perm_domain)
             axes_codomain, axes_domain = bipartition_axes(
                 map(i -> axes(a, i), (perm_codomain..., perm_domain...)),
                 Val(length(perm_codomain))

From b9a884c399bd1ee4ba9deff7bbe59a5394e457fc Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 20:54:17 -0400
Subject: [PATCH 10/17] Say identity, not trivial, for the do-nothing
 bipermutation

`identitybiperm` now matches `isidentitybiperm` and sits beside it. The
symmetry sense of `trivial` is a different concept, so the two never apply
to the same object.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/bituple.jl        | 7 +++++++
 src/factorizations.jl | 4 ++--
 src/matricize.jl      | 6 ------
 3 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/src/bituple.jl b/src/bituple.jl
index f898ee7..60c78d5 100644
--- a/src/bituple.jl
+++ b/src/bituple.jl
@@ -83,3 +83,10 @@ function isidentitybiperm(perm_codomain, perm_domain)
     perm = (perm_codomain..., perm_domain...)
     return perm == ntuple(identity, length(perm))
 end
+
+# The identity bipermutation for a rank-`N` array split after `ndims_codomain` dimensions, i.e.
+# the one `isidentitybiperm` accepts. Transitional: only the `Val` entry points that have yet to be
+# removed build it, to reach the bipermutation hooks.
+function identitybiperm(a, ndims_codomain::Val{K}) where {K}
+    return ntuple(identity, Val(K)), ntuple(i -> K + i, Val(ndims(a) - K))
+end
diff --git a/src/factorizations.jl b/src/factorizations.jl
index ad1a1ca..67dd1e8 100644
--- a/src/factorizations.jl
+++ b/src/factorizations.jl
@@ -870,7 +870,7 @@ true
 function one end
 
 function one!!(style::MatricizeStyle, A, ndims_codomain::Val; kwargs...)
-    A_mat = matricize(style, A, trivialbiperm(A, ndims_codomain)...)
+    A_mat = matricize(style, A, identitybiperm(A, ndims_codomain)...)
     MatrixAlgebraKit.one!(A_mat)
     axes_codomain, axes_domain = bipartition_axes(axes(A), ndims_codomain)
     return unmatricize(style, A_mat, axes_codomain, axes_domain)
@@ -883,7 +883,7 @@ end
 # matricization directly when the style declares one at this split, and otherwise fills a
 # gathered matrix and scatters it back with `unmatricize!`.
 function one!(style::MatricizeStyle, A, ndims_codomain::Val; kwargs...)
-    perm_codomain, perm_domain = trivialbiperm(A, ndims_codomain)
+    perm_codomain, perm_domain = identitybiperm(A, ndims_codomain)
     if is_output_view(matricizeop, style, identity, A, perm_codomain, perm_domain)
         MatrixAlgebraKit.one!(
             matricizeopview(style, identity, A, perm_codomain, perm_domain)
diff --git a/src/matricize.jl b/src/matricize.jl
index aa6b3fe..3bc4acf 100644
--- a/src/matricize.jl
+++ b/src/matricize.jl
@@ -157,12 +157,6 @@ function allocate_output(
     )
 end
 
-# The trivial bipermutation for a rank-`N` array split after `ndims_codomain` dimensions. The
-# `Val` entry points that remain build it to reach the bipermutation hooks.
-function trivialbiperm(a, ndims_codomain::Val{K}) where {K}
-    return ntuple(identity, Val(K)), ntuple(i -> K + i, Val(ndims(a) - K))
-end
-
 # ==================================  is_output_view  ======================================
 # `true` iff `matricizeop(style, op, a, perm_codomain, perm_domain)` shares `a`'s memory, so that
 # writes to the result are writes to `a`. The `isstrided`/`StridedView` pattern, and the same

From facedf4f85eb69ba9e7f15c1d4592ac131dd0305 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 21:33:07 -0400
Subject: [PATCH 11/17] Move the tests onto the matricize hooks

Also migrates five Val-forwarding calls in the factorization wrappers that
a literal-only search had missed, and the second custom style in the
factorization tests.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 ext/TensorAlgebraTensorKitExt.jl | 57 ++++++++++-----------
 src/factorizations.jl            | 10 ++--
 test/test_basics.jl              | 50 +++++++++----------
 test/test_diagonal.jl            |  2 +-
 test/test_exports.jl             |  5 +-
 test/test_factorizations.jl      | 69 +++++++++++++++++++-------
 test/test_matricize.jl           | 85 ++++++++++++++++++--------------
 test/test_matricizestyle.jl      | 28 ++++++++---
 test/test_tensorkitext.jl        |  2 +-
 9 files changed, 181 insertions(+), 127 deletions(-)

diff --git a/ext/TensorAlgebraTensorKitExt.jl b/ext/TensorAlgebraTensorKitExt.jl
index aa770db..58fa709 100644
--- a/ext/TensorAlgebraTensorKitExt.jl
+++ b/ext/TensorAlgebraTensorKitExt.jl
@@ -243,44 +243,41 @@ end
 struct TensorKitMatricize <: TensorAlgebra.MatricizeStyle end
 TensorAlgebra.MatricizeStyle(::Type{<:AbstractTensorMap}) = TensorKitMatricize()
 
-# `permute` at the tensor's own codomain/domain split is trivial and returns `t` itself, so
-# the matching split is the one memory-sharing matricization (TensorKit's own
-# `has_shared_permute` notion); any other split regroups into a fresh `TensorMap`.
-function TensorAlgebra.ismatricizeview(
-        ::TensorKitMatricize, ::AbstractTensorMap{<:Any, <:Any, K}, ::Val{K}
-    ) where {K}
-    return true
+# `permute` at the tensor's own codomain/domain split is trivial and returns `t` itself, so that
+# split is the one memory-sharing matricization (TensorKit's own `has_shared_permute` notion). Any
+# other split, or a folded `conj`, regroups into a fresh `TensorMap`.
+function TensorAlgebra.is_output_view(
+        ::typeof(TensorAlgebra.matricizeop), ::TensorKitMatricize, op,
+        t::AbstractTensorMap, perm_codomain, perm_domain
+    )
+    return op === identity &&
+        TensorAlgebra.isidentitybiperm(perm_codomain, perm_domain) &&
+        length(perm_codomain) == numout(t)
 end
-TensorAlgebra.ismatricizeview(::TensorKitMatricize, ::AbstractTensorMap, ::Val) = false
-function TensorAlgebra.matricizeview(
-        ::TensorKitMatricize, t::AbstractTensorMap{<:Any, <:Any, K}, ::Val{K}
-    ) where {K}
+function TensorAlgebra.matricizeopview(
+        ::TensorKitMatricize, op, t::AbstractTensorMap, perm_codomain, perm_domain
+    )
     return t
 end
-function TensorAlgebra.matricizecopy(
-        ::TensorKitMatricize, t::AbstractTensorMap, ndims_codomain::Val{K}
-    ) where {K}
-    N = numind(t)
-    return permute(
-        t,
-        (ntuple(identity, Val(K)), ntuple(i -> K + i, Val(N - K)));
-        copy = true
+# A `TensorMap`'s matricization is a regrouping of its indices, so the destination is a `TensorMap`
+# over the regrouped space and the write is the ordinary permuted-add. `bipermutedimsopadd!` above
+# routes that through `tensoradd!`, which realizes the permutation, the `op === conj` conjugation
+# and the scaling in one call, so no separate handling of `op` is needed here.
+function TensorAlgebra.allocate_output(
+        ::typeof(TensorAlgebra.matricizeop), ::TensorKitMatricize, op,
+        t::AbstractTensorMap, perm_codomain, perm_domain
     )
+    return similar(t, permute(space(t), (perm_codomain, perm_domain)))
 end
-
-# The identity fill on the regrouped map is TensorKit's own `one!` (MatrixAlgebraKit's
-# `one!` speaks `AbstractMatrix` only).
-function TensorAlgebra.one!!(
-        style::TensorKitMatricize, A::AbstractTensorMap, ndims_codomain::Val; kwargs...
+function TensorAlgebra.matricizeop!(
+        dest::AbstractTensorMap, ::TensorKitMatricize, op,
+        t::AbstractTensorMap, perm_codomain, perm_domain
+    )
+    return TensorAlgebra.bipermutedimsopadd!(
+        dest, op, t, perm_codomain, perm_domain, true, false
     )
-    return TensorKit.one!(TensorAlgebra.matricize(style, A, ndims_codomain))
 end
 
-# `unmatricize` reconstructs the codomain/domain axes from the matrix `m`. A `TensorMap` already
-# is the linear map its space describes, so the only valid request is the one whose codomain/domain
-# split matches `m`'s own space, and `unmatricize` returns `m` unchanged. The domain axes arrive
-# codomain-facing (un-dualized), which is exactly TensorKit's domain convention, so they build the
-# domain `ProductSpace` directly.
 function TensorAlgebra.unmatricize(
         ::TensorKitMatricize, m::AbstractTensorMap, axes_codomain, axes_domain
     )
diff --git a/src/factorizations.jl b/src/factorizations.jl
index 67dd1e8..83bb998 100644
--- a/src/factorizations.jl
+++ b/src/factorizations.jl
@@ -173,7 +173,7 @@ true
 ```
 """
 function tr(style::MatricizeStyle, A, ndims_codomain::Val)
-    return LinearAlgebra.tr(matricize(style, A, ndims_codomain))
+    return LinearAlgebra.tr(matricize(style, A, identitybiperm(A, ndims_codomain)...))
 end
 function tr(A, ndims_codomain::Val)
     return tr(MatricizeStyle(A), A, ndims_codomain)
@@ -559,7 +559,7 @@ The output satisfies `N' * A ≈ 0` and `N' * N ≈ I`.
 left_null
 
 function left_null!!(style::MatricizeStyle, A, ndims_codomain::Val; kwargs...)
-    A_mat = matricize(style, A, ndims_codomain)
+    A_mat = matricize(style, A, identitybiperm(A, ndims_codomain)...)
     N = MatrixAlgebraKit.left_null!(A_mat; kwargs...)
     axes_codomain = first(bipartition(axes(A), ndims_codomain))
     return unmatricize(style, N, axes_codomain, (conj(axes(N, ndims(N))),))
@@ -596,7 +596,7 @@ The output satisfies `A * Nᴴ' ≈ 0` and `Nᴴ * Nᴴ' ≈ I`.
 right_null
 
 function right_null!!(style::MatricizeStyle, A, ndims_codomain::Val; kwargs...)
-    A_mat = matricize(style, A, ndims_codomain)
+    A_mat = matricize(style, A, identitybiperm(A, ndims_codomain)...)
     Nᴴ = MatrixAlgebraKit.right_null!(A_mat; kwargs...)
     _, axes_domain = bipartition_axes(axes(A), ndims_codomain)
     return unmatricize(style, Nᴴ, (axes(Nᴴ, 1),), axes_domain)
@@ -652,7 +652,7 @@ gram_eigh_full
 function gram_eigh_full!!(
         style::MatricizeStyle, A, ndims_codomain::Val; kwargs...
     )
-    A_mat = matricize(style, A, ndims_codomain)
+    A_mat = matricize(style, A, identitybiperm(A, ndims_codomain)...)
     X = MatrixAlgebra.gram_eigh_full!!(A_mat; kwargs...)
     axes_codomain = first(bipartition(axes(A), ndims_codomain))
     return unmatricize(style, X, axes_codomain, (conj(axes(X, ndims(X))),))
@@ -711,7 +711,7 @@ gram_eigh_full_with_pinv
 function gram_eigh_full_with_pinv!!(
         style::MatricizeStyle, A, ndims_codomain::Val; kwargs...
     )
-    A_mat = matricize(style, A, ndims_codomain)
+    A_mat = matricize(style, A, identitybiperm(A, ndims_codomain)...)
     X, Y = MatrixAlgebra.gram_eigh_full_with_pinv!!(A_mat; kwargs...)
     axes_codomain = first(bipartition(axes(A), ndims_codomain))
     return unmatricize(style, X, axes_codomain, (conj(axes(X, ndims(X))),)),
diff --git a/test/test_basics.jl b/test/test_basics.jl
index 8f45811..b0a119c 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -1,7 +1,7 @@
 import TensorAlgebra
 using StableRNGs: StableRNG
 using TensorAlgebra: BiTuple, ContractAlgorithm, bipermutedims, bipermutedims!, contract,
-    contract!, contractadd!, length_codomain, length_domain, matricizeperm, unmatricize,
+    contract!, contractadd!, length_codomain, length_domain, matricize, unmatricize,
     unmatricize!
 using TensorOperations: TensorOperations
 using Test: @test, @test_broken, @test_throws, @testset
@@ -53,68 +53,68 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
     @testset "matricize (eltype=$elt)" for elt in elts
         a = randn(elt, 2, 3, 4, 5)
 
-        a_fused = matricizeperm(a, (1, 2), (3, 4))
+        a_fused = matricize(a, (1, 2), (3, 4))
         @test eltype(a_fused) === elt
         @test a_fused ≈ reshape(a, 6, 20)
-        a_fused = matricizeperm(a, (3, 1), (2, 4))
+        a_fused = matricize(a, (3, 1), (2, 4))
         @test eltype(a_fused) === elt
         @test a_fused ≈ reshape(permutedims(a, (3, 1, 2, 4)), (8, 15))
-        a_fused = matricizeperm(a, (3, 1, 2), (4,))
+        a_fused = matricize(a, (3, 1, 2), (4,))
         @test eltype(a_fused) === elt
         @test a_fused ≈ reshape(permutedims(a, (3, 1, 2, 4)), (24, 5))
-        a_fused = matricizeperm(a, (2, 4), (3, 1))
+        a_fused = matricize(a, (2, 4), (3, 1))
         @test eltype(a_fused) === elt
         @test a_fused ≈ reshape(permutedims(a, (2, 4, 3, 1)), (15, 8))
 
         # Degenerate splits: everything in the domain, then everything in the codomain.
-        a_fused = matricizeperm(a, (), (1, 2, 3, 4))
+        a_fused = matricize(a, (), (1, 2, 3, 4))
         @test eltype(a_fused) === elt
         @test a_fused ≈ reshape(a, (1, 120))
-        a_fused = matricizeperm(a, (1, 2, 3, 4), ())
+        a_fused = matricize(a, (1, 2, 3, 4), ())
         @test eltype(a_fused) === elt
         @test a_fused ≈ reshape(a, (120, 1))
 
-        @test_throws MethodError matricizeperm(a, (1, 2), (3,), (4,))
-        @test_throws MethodError matricizeperm(a, (1, 2, 3, 4))
-        @test_throws ArgumentError matricizeperm(a, (1, 2), (3,))
+        @test_throws MethodError matricize(a, (1, 2), (3,), (4,))
+        @test_throws MethodError matricize(a, (1, 2, 3, 4))
+        @test_throws ArgumentError matricize(a, (1, 2), (3,))
 
         v = ones(elt, 2)
-        a_fused = matricizeperm(v, (1,), ())
+        a_fused = matricize(v, (1,), ())
         @test eltype(a_fused) === elt
         @test a_fused ≈ ones(elt, 2, 1)
-        a_fused = matricizeperm(v, (), (1,))
+        a_fused = matricize(v, (), (1,))
         @test eltype(a_fused) === elt
         @test a_fused ≈ ones(elt, 1, 2)
 
-        a_fused = matricizeperm(ones(elt), (), ())
+        a_fused = matricize(ones(elt), (), ())
         @test eltype(a_fused) === elt
         @test a_fused ≈ ones(elt, 1, 1)
     end
 
-    @testset "matricizeopperm (eltype=$elt)" for elt in elts
+    @testset "matricizeop (eltype=$elt)" for elt in elts
         rng = StableRNG(123)
         a = randn(rng, elt, 2, 3, 4)
 
         # identity op: should match matricize exactly
-        m = TensorAlgebra.matricizeopperm(identity, a, (1,), (2, 3))
-        m_ref = matricizeperm(a, (1,), (2, 3))
+        m = TensorAlgebra.matricizeop(identity, a, (1,), (2, 3))
+        m_ref = matricize(a, (1,), (2, 3))
         @test m ≈ m_ref
 
-        m = TensorAlgebra.matricizeopperm(identity, a, (3, 1), (2,))
-        m_ref = matricizeperm(a, (3, 1), (2,))
+        m = TensorAlgebra.matricizeop(identity, a, (3, 1), (2,))
+        m_ref = matricize(a, (3, 1), (2,))
         @test m ≈ m_ref
 
-        m = TensorAlgebra.matricizeopperm(identity, a, (2, 3), (1,))
-        m_ref = matricizeperm(a, (2, 3), (1,))
+        m = TensorAlgebra.matricizeop(identity, a, (2, 3), (1,))
+        m_ref = matricize(a, (2, 3), (1,))
         @test m ≈ m_ref
 
         # conj op
-        m = TensorAlgebra.matricizeopperm(conj, a, (1,), (2, 3))
-        m_ref = conj.(matricizeperm(a, (1,), (2, 3)))
+        m = TensorAlgebra.matricizeop(conj, a, (1,), (2, 3))
+        m_ref = conj.(matricize(a, (1,), (2, 3)))
         @test m ≈ m_ref
 
-        m = TensorAlgebra.matricizeopperm(conj, a, (3, 1), (2,))
-        m_ref = conj.(matricizeperm(a, (3, 1), (2,)))
+        m = TensorAlgebra.matricizeop(conj, a, (3, 1), (2,))
+        m_ref = conj.(matricize(a, (3, 1), (2,)))
         @test m ≈ m_ref
     end
 
@@ -136,7 +136,7 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         unmatricize!(a, m, (1, 2), (3, 4))
         @test a ≈ a0
 
-        m1 = matricizeperm(a0, perm_codomain, perm_domain)
+        m1 = matricize(a0, perm_codomain, perm_domain)
         a = similar(a0)
         unmatricize!(a, m1, perm_codomain, perm_domain)
         @test a ≈ a0
diff --git a/test/test_diagonal.jl b/test/test_diagonal.jl
index c492918..45fcefe 100644
--- a/test/test_diagonal.jl
+++ b/test/test_diagonal.jl
@@ -41,7 +41,7 @@ using Test: @test, @test_throws, @testset
     end
 
     @testset "matricize(1, 1) is the identity reshape" begin
-        m = TensorAlgebra.matricize(TensorAlgebra.ReshapeMatricize(), d, Val(1))
+        m = TensorAlgebra.matricize(TensorAlgebra.ReshapeMatricize(), d, (1,), (2,))
         @test m === d
     end
 
diff --git a/test/test_exports.jl b/test/test_exports.jl
index 5bab6ad..0ff6b33 100644
--- a/test/test_exports.jl
+++ b/test/test_exports.jl
@@ -43,8 +43,9 @@ using Test: @test, @testset
                 :biperm, :bipartition, :cat_similar,
                 :concatenate, :concatenate!, :ContractAlgorithm, :contractopadd!, :data,
                 :datatype, :directsum,
-                :flattenlinear, :label_type,
-                :matricizeopperm, :permutedims, :permutedims!, :scalar, :similar_map,
+                :flattenlinear, :is_output_view, :label_type,
+                :matricize, :matricizeop, :matricizeop!, :matricizeopcopy, :matricizeopview,
+                :permutedims, :permutedims!, :scalar, :similar_map,
                 :TensorOperationsAlgorithm,
                 :to_range, :tr, :tryflattenlinear, :ungrade, :zero!, :scale!,
                 :permuteddims, :PermutedDims,
diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl
index f893445..d69a007 100644
--- a/test/test_factorizations.jl
+++ b/test/test_factorizations.jl
@@ -7,6 +7,9 @@ using TensorAlgebra: TensorAlgebra, contract, eig_full, eig_vals, eigh_full, eig
 using Test: @test, @testset
 using TestExtras: @constinferred
 
+# Matricize without permuting: the identity bipermutation for a split after `k` dimensions.
+splitperms(a, k) = (ntuple(identity, k), ntuple(i -> k + i, ndims(a) - k))
+
 elts = (Float64, ComplexF64)
 
 # QR Decomposition
@@ -373,7 +376,7 @@ end
     @test size(Id) == size(A)
     @test eltype(Id) === T
 
-    @test TensorAlgebra.matricize(Id, Val(2)) ≈ I
+    @test TensorAlgebra.matricize(Id, splitperms(Id, 2)...) ≈ I
 
     # `Val`, perm, and label entries agree.
     @test TensorAlgebra.one(A, Val(2)) ≈ Id
@@ -386,7 +389,7 @@ end
     B = randn(T, 2, 2, 2, 2)
     labels_B = (:a, :c, :b, :d)
     Id_perm = TensorAlgebra.one(B, labels_B, labels_cod, labels_dom)
-    @test TensorAlgebra.matricize(Id_perm, Val(2)) ≈ I
+    @test TensorAlgebra.matricize(Id_perm, splitperms(Id_perm, 2)...) ≈ I
     # Perm- and biperm-tuple forms agree with the label form.
     @test TensorAlgebra.one(B, (1, 3), (2, 4)) ≈ Id_perm
 
@@ -394,12 +397,12 @@ end
     C = randn(T, 2, 3, 2, 3)
     Cret = @constinferred TensorAlgebra.one!(C, Val(2))
     @test Cret === C
-    @test TensorAlgebra.matricize(C, Val(2)) ≈ I
+    @test TensorAlgebra.matricize(C, splitperms(C, 2)...) ≈ I
     @test C ≈ TensorAlgebra.one(A, Val(2))
 
     # `unmatricize!` scatters a fused matrix back into an existing array.
     D = randn(T, 2, 3, 2, 3)
-    Dmat = TensorAlgebra.matricize(D, Val(2))
+    Dmat = TensorAlgebra.matricize(D, splitperms(D, 2)...)
     E = similar(D)
     Eret = TensorAlgebra.unmatricize!(E, Dmat, Val(2))
     @test Eret === E
@@ -430,23 +433,25 @@ end
         (((1, 2), (3,)), ((3, 1), (2,)), ((3,), (1, 2)), ((2,), (3, 1)))
         k = length(perm_codomain)
         A_perm = TensorAlgebra.bipermutedims(A, perm_codomain, perm_domain)
-        A_mat = TensorAlgebra.matricize(A_perm, Val(k))
+        A_mat = TensorAlgebra.matricize(A_perm, splitperms(A_perm, k)...)
         for f in (qr_compact, lq_compact, left_orth, right_orth)
             X, Y = f(A, perm_codomain, perm_domain)
-            @test TensorAlgebra.matricize(X, Val(k)) *
-                TensorAlgebra.matricize(Y, Val(1)) ≈ A_mat
+            @test TensorAlgebra.matricize(X, splitperms(X, k)...) *
+                TensorAlgebra.matricize(Y, splitperms(Y, 1)...) ≈ A_mat
         end
         for f in (svd_compact, svd_trunc)
             U, S, Vᴴ = f(A, perm_codomain, perm_domain)
-            U_mat = TensorAlgebra.matricize(U, Val(k))
-            @test U_mat * S * TensorAlgebra.matricize(Vᴴ, Val(1)) ≈ A_mat
+            U_mat = TensorAlgebra.matricize(U, splitperms(U, k)...)
+            @test U_mat * S * TensorAlgebra.matricize(Vᴴ, splitperms(Vᴴ, 1)...) ≈ A_mat
             @test U_mat' * U_mat ≈ I
         end
         @test svd_vals(A, perm_codomain, perm_domain) ≈ LinearAlgebra.svdvals(A_mat)
-        N = TensorAlgebra.matricize(left_null(A, perm_codomain, perm_domain), Val(k))
+        N_tensor = left_null(A, perm_codomain, perm_domain)
+        N = TensorAlgebra.matricize(N_tensor, splitperms(N_tensor, k)...)
         @test norm(N' * A_mat) ≈ 0 atol = 1.0e-13
         @test N' * N ≈ I
-        Nᴴ = TensorAlgebra.matricize(right_null(A, perm_codomain, perm_domain), Val(1))
+        Nᴴ_tensor = right_null(A, perm_codomain, perm_domain)
+        Nᴴ = TensorAlgebra.matricize(Nᴴ_tensor, splitperms(Nᴴ_tensor, 1)...)
         @test norm(A_mat * Nᴴ') ≈ 0 atol = 1.0e-13
         @test Nᴴ * Nᴴ' ≈ I
         @test A == Acopy
@@ -456,9 +461,9 @@ end
     for (perm_codomain, perm_domain) in
         (((1, 2), (3, 4)), ((3, 4), (1, 2)), ((2, 3), (4, 1)))
         B_perm = TensorAlgebra.bipermutedims(B, perm_codomain, perm_domain)
-        B_mat = Matrix(TensorAlgebra.matricize(B_perm, Val(2)))
+        B_mat = Matrix(TensorAlgebra.matricize(B_perm, splitperms(B_perm, 2)...))
         D, V = eig_full(B, perm_codomain, perm_domain)
-        V_mat = TensorAlgebra.matricize(V, Val(2))
+        V_mat = TensorAlgebra.matricize(V, splitperms(V, 2)...)
         @test B_mat * V_mat ≈ V_mat * D
         sortvals(v) = sort(v; by = x -> (real(x), imag(x)))
         @test sortvals(eig_vals(B, perm_codomain, perm_domain)) ≈
@@ -481,11 +486,37 @@ module FactorizationMatricizeTestUtils
     end
     struct AliasingMatricize <: TA.MatricizeStyle end
     TA.MatricizeStyle(::Type{<:AliasingArray}) = AliasingMatricize()
-    function TA.matricize(::AliasingMatricize, a::AliasingArray, ndims_codomain::Val)
-        return TA.matricize(TA.ReshapeMatricize(), a.parent, ndims_codomain)
+    # Delegate every hook to the dense style on the unwrapped parent, so the matricization
+    # aliases exactly where a plain `Array`'s would.
+    unwrap(a::AliasingArray) = a.parent
+    unwrap(a::AbstractArray) = a
+    function TA.is_output_view(
+            ::typeof(TA.matricizeop), ::AliasingMatricize, op, a, perm_codomain, perm_domain
+        )
+        return TA.is_output_view(
+            TA.matricizeop, TA.ReshapeMatricize(), op, unwrap(a), perm_codomain, perm_domain
+        )
+    end
+    function TA.matricizeopview(
+            ::AliasingMatricize, op, a, perm_codomain, perm_domain
+        )
+        return TA.matricizeopview(
+            TA.ReshapeMatricize(), op, unwrap(a), perm_codomain, perm_domain
+        )
     end
-    function TA.matricize(::AliasingMatricize, a::AbstractArray, ndims_codomain::Val)
-        return TA.matricize(TA.ReshapeMatricize(), a, ndims_codomain)
+    function TA.allocate_output(
+            ::typeof(TA.matricizeop), ::AliasingMatricize, op, a, perm_codomain, perm_domain
+        )
+        return TA.allocate_output(
+            TA.matricizeop, TA.ReshapeMatricize(), op, unwrap(a), perm_codomain, perm_domain
+        )
+    end
+    function TA.matricizeop!(
+            dest, ::AliasingMatricize, op, a, perm_codomain, perm_domain
+        )
+        return TA.matricizeop!(
+            dest, TA.ReshapeMatricize(), op, unwrap(a), perm_codomain, perm_domain
+        )
     end
     function TA.unmatricize(::AliasingMatricize, m, axes_codomain, axes_domain)
         return AliasingArray(
@@ -512,9 +543,9 @@ end
     parent_copy = copy(parent)
     Q, R = qr_compact(A, Val(1))
     @test parent == parent_copy
-    Q_mat = TensorAlgebra.matricize(Q, Val(1))
+    Q_mat = TensorAlgebra.matricize(Q, splitperms(Q, 1)...)
     @test eltype(Q_mat) === Float64
-    @test Q_mat * TensorAlgebra.matricize(R, Val(1)) ≈ reshape(parent, 2, 12)
+    @test Q_mat * TensorAlgebra.matricize(R, splitperms(R, 1)...) ≈ reshape(parent, 2, 12)
     @test svd_vals(A, (1,), (2, 3)) ≈ LinearAlgebra.svdvals(reshape(float.(parent), 2, 12))
     @test parent == parent_copy
 end
diff --git a/test/test_matricize.jl b/test/test_matricize.jl
index d16eba7..0582b4b 100644
--- a/test/test_matricize.jl
+++ b/test/test_matricize.jl
@@ -1,6 +1,6 @@
 using StableRNGs: StableRNG
-using TensorAlgebra: TensorAlgebra, ReshapeMatricize, ismatricizeview, matricize,
-    matricizecopy, matricizeopperm, matricizeperm, matricizeview
+using TensorAlgebra: TensorAlgebra, ReshapeMatricize, is_output_view, matricize,
+    matricizeop, matricizeop!, matricizeopcopy, matricizeopview
 using Test: @test, @test_throws, @testset
 
 # A non-`ReshapeMatricize` style, to check the always-safe generic fallback.
@@ -14,86 +14,95 @@ function matricize_ref(a, perm_codomain, perm_domain)
     return reshape(a_perm, (nrow, ncol))
 end
 
-@testset "maybe-view matricizeperm (eltype=$elt)" for elt in (Float64, ComplexF64)
+@testset "maybe-view matricize (eltype=$elt)" for elt in (Float64, ComplexF64)
     a = randn(StableRNG(123), elt, 2, 3, 4)
 
     # Identity bipermutation: correct values and a view aliasing `a`.
-    m = matricizeperm(a, (1,), (2, 3))
+    m = matricize(a, (1,), (2, 3))
     @test m ≈ matricize_ref(a, (1,), (2, 3))
     @test Base.mightalias(m, a)
 
     # Every other bipermutation is a fresh permuted copy in matricized layout (no lazy
     # wrappers), including the codomain/domain swap.
     for (pc, pd) in (((2, 3), (1,)), ((3, 1), (2,)))
-        m = matricizeperm(a, pc, pd)
+        m = matricize(a, pc, pd)
         @test m ≈ matricize_ref(a, pc, pd)
         @test m isa Matrix
         @test !Base.mightalias(m, a)
     end
-    @test_throws ArgumentError matricizeperm(a, (1,), (2,))
+    @test_throws ArgumentError matricize(a, (1,), (2,))
 
     # `conj` cannot ride a view, so it copies even on the identity bipermutation.
-    m = matricizeopperm(conj, a, (1,), (2, 3))
+    m = matricizeop(conj, a, (1,), (2, 3))
     @test m ≈ conj.(matricize_ref(a, (1,), (2, 3)))
     @test !Base.mightalias(m, a)
-    m = matricizeopperm(conj, a, (2, 3), (1,))
+    m = matricizeop(conj, a, (2, 3), (1,))
     @test m ≈ conj.(matricize_ref(a, (2, 3), (1,)))
     @test !Base.mightalias(m, a)
 end
 
-@testset "ismatricizeview" begin
+@testset "is_output_view" begin
     a = randn(StableRNG(321), 2, 3, 4)
     style = ReshapeMatricize()
 
-    # A dense reshape matricization shares memory at every trivial split.
-    @test ismatricizeview(style, a, Val(1))
-    @test ismatricizeview(style, a, (1,), (2, 3))
+    # A dense reshape shares memory at the identity bipermutation.
+    @test is_output_view(matricizeop, style, identity, a, (1,), (2, 3))
+    @test is_output_view(matricizeop, style, identity, a, (), (1, 2, 3))
+    @test is_output_view(matricizeop, style, identity, a, (1, 2, 3), ())
 
-    # The bipermutation form declares sharing only at the identity: a swap or interleaving
-    # bipermutation routes through the consumers' gather branches.
-    @test !ismatricizeview(style, a, (2, 3), (1,))
-    @test !ismatricizeview(style, a, (3, 1), (2,))
+    # Not at a swap or an interleaving, which route through the gather branch instead.
+    @test !is_output_view(matricizeop, style, identity, a, (2, 3), (1,))
+    @test !is_output_view(matricizeop, style, identity, a, (3, 1), (2,))
+
+    # And never with an operation folded in, since a reshape cannot carry a `conj`.
+    @test !is_output_view(matricizeop, style, conj, a, (1,), (2, 3))
 
     # A generic style declares nothing (fail-safe default).
-    @test !ismatricizeview(DummyMatricize(), a, Val(1))
-    @test !ismatricizeview(DummyMatricize(), a, (1,), (2, 3))
+    @test !is_output_view(matricizeop, DummyMatricize(), identity, a, (1,), (2, 3))
 
     # Writes to the shared matricization are writes to `a`.
-    m = matricizeview(style, a, Val(1))
+    m = matricizeopview(style, identity, a, (1,), (2, 3))
     @test m == matricize_ref(a, (1,), (2, 3))
     m[1, 1] = 42
     @test a[1, 1, 1] == 42
 end
 
-@testset "ismatricizeview coherence" begin
+@testset "is_output_view coherence" begin
     rng = StableRNG(11)
     a = randn(rng, 2, 3, 4)
     style = ReshapeMatricize()
 
-    # A declared share means `matricizeview` (and so `matricize`) aliases `a`, while
-    # `matricizecopy` never does.
+    # A declared share means `matricizeopview` (and so `matricize`) aliases `a`, while
+    # `matricizeopcopy` never does.
     for K in 0:3
-        if ismatricizeview(style, a, Val(K))
-            m = matricizeview(style, a, Val(K))
-            @test Base.mightalias(m, a)
-            @test matricize(style, a, Val(K)) == m
-        end
-        @test !Base.mightalias(matricizecopy(style, a, Val(K)), a)
-
-        # The perm form of the copy leaf is owned too, and at the trivial bipermutation it
-        # matches the `Val` form.
         pc = ntuple(identity, K)
         pd = ntuple(i -> K + i, 3 - K)
-        m_perm = matricizecopy(style, a, pc, pd)
-        @test !Base.mightalias(m_perm, a)
-        @test m_perm == matricizecopy(style, a, Val(K))
+        if is_output_view(matricizeop, style, identity, a, pc, pd)
+            m = matricizeopview(style, identity, a, pc, pd)
+            @test Base.mightalias(m, a)
+            @test matricize(style, a, pc, pd) == m
+        end
+        m_copy = matricizeopcopy(style, identity, a, pc, pd)
+        @test !Base.mightalias(m_copy, a)
+        @test m_copy == matricize_ref(a, pc, pd)
     end
     for (pc, pd) in (((2, 3), (1,)), ((3, 1), (2,)))
-        m = matricizecopy(style, a, pc, pd)
+        m = matricizeopcopy(style, identity, a, pc, pd)
         @test m ≈ matricize_ref(a, pc, pd)
         @test !Base.mightalias(m, a)
     end
-    @test_throws ArgumentError matricizecopy(style, a, (1,), (2,))
+    @test_throws ArgumentError matricizeopcopy(style, identity, a, (1,), (2,))
+
+    # The allocation and write hooks compose into the copy form.
+    for (pc, pd) in (((1,), (2, 3)), ((3, 1), (2,)))
+        for op in (identity, conj)
+            dest = TensorAlgebra.allocate_output(matricizeop, style, op, a, pc, pd)
+            @test size(dest) == size(matricize_ref(a, pc, pd))
+            matricizeop!(dest, style, op, a, pc, pd)
+            @test dest ≈ op.(matricize_ref(a, pc, pd))
+            @test dest ≈ matricizeopcopy(style, op, a, pc, pd)
+        end
+    end
 
     # Both destination branches of a consumer (`contractadd!`) behave: the shared-view route
     # for the identity destination bipermutation and the gather/scatter route otherwise.
@@ -123,14 +132,14 @@ end
 
     # Identity-bipermutation view tracks an in-place update of `a`.
     a = randn(rng, 2, 3, 4)
-    m = matricizeperm(a, (1,), (2, 3))
+    m = matricize(a, (1,), (2, 3))
     a .= randn(rng, 2, 3, 4)
     @test m ≈ matricize_ref(a, (1,), (2, 3))
 
     # Permuted copies are independent of later updates to `a`.
     for (pc, pd) in (((2, 3), (1,)), ((3, 1), (2,)))
         a = randn(rng, 2, 3, 4)
-        m = matricizeperm(a, pc, pd)
+        m = matricize(a, pc, pd)
         snapshot = copy(m)
         a .= a .+ 1
         @test m == snapshot
diff --git a/test/test_matricizestyle.jl b/test/test_matricizestyle.jl
index 8d13c40..3c623e2 100644
--- a/test/test_matricizestyle.jl
+++ b/test/test_matricizestyle.jl
@@ -9,12 +9,28 @@ module MatricizeStyleTestUtils
     end
     struct MyArrayMatricize <: TA.MatricizeStyle end
     TA.MatricizeStyle(::Type{<:MyArray}) = MyArrayMatricize()
-    # Minimal fold/unfold leaves so a round-trip (`one!`) can run through the custom style:
-    # both dispatch on `MyArrayMatricize`, so an unfold whose style was re-derived from the
-    # plain fused matrix instead of threaded through would miss them and error.
-    TA.ismatricizeview(::MyArrayMatricize, a, ::Val) = false
-    function TA.matricizecopy(::MyArrayMatricize, a::MyArray, ndims_codomain::Val)
-        return TA.matricizecopy(TA.ReshapeMatricize(), a.parent, ndims_codomain)
+    # Minimal hooks so a round trip (`one!`) can run through the custom style. All of them
+    # dispatch on `MyArrayMatricize`, so a path whose style was re-derived from the plain fused
+    # matrix instead of threaded through would miss them and error.
+    function TA.is_output_view(
+            ::typeof(TA.matricizeop), ::MyArrayMatricize, op, a, perm_codomain, perm_domain
+        )
+        return false
+    end
+    function TA.allocate_output(
+            ::typeof(TA.matricizeop), ::MyArrayMatricize, op, a::MyArray,
+            perm_codomain, perm_domain
+        )
+        return TA.allocate_output(
+            TA.matricizeop, TA.ReshapeMatricize(), op, a.parent, perm_codomain, perm_domain
+        )
+    end
+    function TA.matricizeop!(
+            dest, ::MyArrayMatricize, op, a::MyArray, perm_codomain, perm_domain
+        )
+        return TA.matricizeop!(
+            dest, TA.ReshapeMatricize(), op, a.parent, perm_codomain, perm_domain
+        )
     end
     function TA.unmatricize!(
             ::MyArrayMatricize, a_dest::MyArray, m,
diff --git a/test/test_tensorkitext.jl b/test/test_tensorkitext.jl
index 557b498..30b36a2 100644
--- a/test/test_tensorkitext.jl
+++ b/test/test_tensorkitext.jl
@@ -60,7 +60,7 @@ using Test: @test, @test_throws, @testset
         # `unmatricize` takes the domain axes codomain-facing (un-dualized), so pass `B`, `C1`
         # directly rather than the dualized `space(t, 3)`, `space(t, 4)`.
         axes_domain = (B, C1)
-        m = matricize(t, Val(2))
+        m = matricize(t, (1, 2), (3, 4))
         @test space(m) == space(t)
         back = unmatricize(m, axes_codomain, axes_domain)
         @test back ≈ t

From 21eb4618f4b6a5098a654be09e969920fcedecc7 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 22:35:37 -0400
Subject: [PATCH 12/17] Offer the split-only spelling as convenience again

Matricizing without permuting is worth a short spelling. Removing `Val` as
a dispatch tier was what mattered: a style implements the bipermutation
hooks and never these, so the copy path cannot recurse through the router.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/bituple.jl   |  4 ++--
 src/matricize.jl | 16 ++++++++++++++++
 2 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/src/bituple.jl b/src/bituple.jl
index 60c78d5..d9344a5 100644
--- a/src/bituple.jl
+++ b/src/bituple.jl
@@ -85,8 +85,8 @@ function isidentitybiperm(perm_codomain, perm_domain)
 end
 
 # The identity bipermutation for a rank-`N` array split after `ndims_codomain` dimensions, i.e.
-# the one `isidentitybiperm` accepts. Transitional: only the `Val` entry points that have yet to be
-# removed build it, to reach the bipermutation hooks.
+# the one `isidentitybiperm` accepts. The split-only `Val` conveniences build it to reach the
+# bipermutation forms.
 function identitybiperm(a, ndims_codomain::Val{K}) where {K}
     return ntuple(identity, Val(K)), ntuple(i -> K + i, Val(ndims(a) - K))
 end
diff --git a/src/matricize.jl b/src/matricize.jl
index 3bc4acf..23cff0f 100644
--- a/src/matricize.jl
+++ b/src/matricize.jl
@@ -120,6 +120,22 @@ function matricize(style::MatricizeStyle, a, perm_codomain, perm_domain)
     return matricizeop(style, identity, a, perm_codomain, perm_domain)
 end
 
+# Split-only convenience: matricize after `ndims_codomain` dimensions without permuting. Sugar over
+# the bipermutation forms, not a dispatch tier. A style implements the hooks above and never these,
+# which is what keeps the copy path from recursing back through the router.
+function matricize(a, ndims_codomain::Val)
+    return matricize(a, identitybiperm(a, ndims_codomain)...)
+end
+function matricize(style::MatricizeStyle, a, ndims_codomain::Val)
+    return matricize(style, a, identitybiperm(a, ndims_codomain)...)
+end
+function matricizeop(op, a, ndims_codomain::Val)
+    return matricizeop(op, a, identitybiperm(a, ndims_codomain)...)
+end
+function matricizeop(style::MatricizeStyle, op, a, ndims_codomain::Val)
+    return matricizeop(style, op, a, identitybiperm(a, ndims_codomain)...)
+end
+
 # Total: always fresh storage the caller owns.
 function matricizeopcopy(op, a, perm_codomain, perm_domain)
     return matricizeopcopy(MatricizeStyle(a), op, a, perm_codomain, perm_domain)

From 9e41545562efe649166e93858ac13e244dce57d7 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 22:48:42 -0400
Subject: [PATCH 13/17] Remove gram_eigh_full, gram_eigh_full_with_pinv, and
 sqrth_invsqrth_safe

The gram factorizations were only used by higher-level network code, so they belong in the package that needs them. `sqrth_invsqrth_safe` only saves an eigendecomposition over calling the two separately.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/MatrixAlgebra.jl        | 121 +--------------------------
 src/TensorAlgebra.jl        |   4 +-
 src/factorizations.jl       | 158 ++----------------------------------
 test/test_exports.jl        |   6 --
 test/test_factorizations.jl |  58 +------------
 test/test_matrixalgebra.jl  |  36 --------
 6 files changed, 11 insertions(+), 372 deletions(-)

diff --git a/src/MatrixAlgebra.jl b/src/MatrixAlgebra.jl
index 82aeece..05e5aba 100644
--- a/src/MatrixAlgebra.jl
+++ b/src/MatrixAlgebra.jl
@@ -1,15 +1,12 @@
 module MatrixAlgebra
 
-export gram_eigh_full,
-    gram_eigh_full_with_pinv,
-    invsqrt_diag_safe,
+export invsqrt_diag_safe,
     invsqrth_safe,
     pow_diag_safe,
     pow_diag_safe!,
     powh_safe,
     sqrt_diag_safe,
-    sqrth_safe,
-    sqrth_invsqrth_safe
+    sqrth_safe
 
 using LinearAlgebra: LinearAlgebra, Diagonal, isdiag, norm
 using MatrixAlgebraKit: MatrixAlgebraKit as MAK
@@ -165,120 +162,6 @@ $(_clamp_kwargs_doc("M"))
 """
 invsqrth_safe(M; kwargs...) = powh_safe(M, -1 // 2; kwargs...)
 
-"""
-    sqrth_invsqrth_safe(M; alg=nothing, atol=0, rtol=eps(real(eltype(M)))^(3//4)) -> M^(1//2), M^(-1//2)
-
-Square root and pseudo-inverse square root of a Hermitian positive
-semi-definite matrix, from a single eigendecomposition. Equivalent
-to `(sqrth_safe(M; ...), invsqrth_safe(M; ...))` but with the
-eigendecomposition computed once. Eigenvalues below tolerance are clamped
-to zero in both factors (Moore-Penrose convention for the inverse).
-
-The input must be Hermitian (as for `MatrixAlgebraKit.eigh_full`): project
-with `MatrixAlgebraKit.project_hermitian` first if it is Hermitian only up
-to numerical noise.
-
-## Keyword arguments
-
-  - `alg`: forwarded to `MatrixAlgebraKit.eigh_full`.
-
-$(_clamp_kwargs_doc("M"))
-"""
-function sqrth_invsqrth_safe(M; alg = nothing, kwargs...)
-    if isdiag(M)
-        return pow_diag_safe(M, 1 // 2; kwargs...), pow_diag_safe(M, -1 // 2; kwargs...)
-    end
-    D, V = MAK.eigh_full(M; alg)
-    return V * pow_diag_safe(D, 1 // 2; kwargs...) * V',
-        V * pow_diag_safe(D, -1 // 2; kwargs...) * V'
-end
-
-for (gram, gram_with_pinv, eigh_full) in (
-        (:gram_eigh_full, :gram_eigh_full_with_pinv, :eigh_full),
-        (:gram_eigh_full!!, :gram_eigh_full_with_pinv!!, :eigh_full!),
-    )
-    @eval begin
-        function $gram(A::AbstractMatrix; alg = nothing, kwargs...)
-            D, V = MAK.$eigh_full(A; alg)
-            return V * sqrth_safe(D; kwargs...)
-        end
-        function $gram_with_pinv(A::AbstractMatrix; alg = nothing, kwargs...)
-            D, V = MAK.$eigh_full(A; alg)
-            return V * sqrth_safe(D; kwargs...), invsqrth_safe(D; kwargs...) * V'
-        end
-    end
-end
-
-"""
-    gram_eigh_full(A::AbstractMatrix; alg=nothing, atol=0, rtol=eps(real(eltype(A)))^(3//4)) -> X
-
-Gram factorization of a Hermitian positive semi-definite matrix via its
-eigendecomposition (balanced eigh): returns `X = V * sqrth_safe(D; atol, rtol)`
-such that `A ≈ X * X'`, where `A = V * D * V'`. The square-root of `D` is
-absorbed symmetrically into the two factors of the eigendecomposition.
-Eigenvalues below `tol` (see [`pow_diag_safe`](@ref)) are clamped to zero.
-
-## Keyword arguments
-
-  - `alg`: forwarded to `MatrixAlgebraKit.eigh_full`.
-
-$(_clamp_kwargs_doc("A"))
-
-# Examples
-
-```jldoctest
-julia> using TensorAlgebra.MatrixAlgebra: gram_eigh_full
-
-julia> B = [1.0 0.5; 0.5 2.0];
-
-julia> A = B' * B;
-
-julia> X = gram_eigh_full(A);
-
-julia> X * X' ≈ A
-true
-```
-
-See also [`gram_eigh_full_with_pinv`](@ref).
-"""
-gram_eigh_full
-
-"""
-    gram_eigh_full_with_pinv(A::AbstractMatrix; alg=nothing, atol=0, rtol=eps(real(eltype(A)))^(3//4)) -> X, Y
-
-Like [`gram_eigh_full`](@ref), but additionally returns
-`Y = invsqrth_safe(D; atol, rtol) * V' ≈ pinv(X)`, a left inverse of `X`
-on the rank subspace: `Y * X ≈ I`. Eigenvalues below `tol` are clamped to
-zero in both factors.
-
-## Keyword arguments
-
-  - `alg`: forwarded to `MatrixAlgebraKit.eigh_full`.
-
-$(_clamp_kwargs_doc("A"))
-
-# Examples
-
-```jldoctest
-julia> using LinearAlgebra: I
-
-julia> using TensorAlgebra.MatrixAlgebra: gram_eigh_full_with_pinv
-
-julia> B = [1.0 0.5; 0.5 2.0];
-
-julia> A = B' * B;
-
-julia> X, Y = gram_eigh_full_with_pinv(A);
-
-julia> X * X' ≈ A
-true
-
-julia> Y * X ≈ I
-true
-```
-"""
-gram_eigh_full_with_pinv
-
 using MatrixAlgebraKit: MatrixAlgebraKit, TruncationStrategy
 
 struct TruncationDegenerate{Strategy <: TruncationStrategy, T <: Real} <: TruncationStrategy
diff --git a/src/TensorAlgebra.jl b/src/TensorAlgebra.jl
index 0ac1bc9..6ae6d18 100644
--- a/src/TensorAlgebra.jl
+++ b/src/TensorAlgebra.jl
@@ -1,9 +1,9 @@
 module TensorAlgebra
 
 export contract, contract!, dual, eig_full, eig_trunc, eig_vals, eigh_full, eigh_trunc,
-    eigh_vals, gram_eigh_full, gram_eigh_full_with_pinv, invsqrth_safe, isdual, left_null,
+    eigh_vals, invsqrth_safe, isdual, left_null,
     left_orth, left_polar, lq_compact, lq_full, project_hermitian, qr_compact,
-    qr_full, right_null, right_orth, right_polar, sqrth_invsqrth_safe, sqrth_safe,
+    qr_full, right_null, right_orth, right_polar, sqrth_safe,
     svd_compact, svd_full, svd_trunc, svd_vals
 
 if VERSION >= v"1.11.0-DEV.469"
diff --git a/src/factorizations.jl b/src/factorizations.jl
index 83bb998..87fbad5 100644
--- a/src/factorizations.jl
+++ b/src/factorizations.jl
@@ -63,8 +63,7 @@ end
 # Read-only tier: the matrix-level entries never mutate their input (they copy internally), so
 # the perm form consumes the maybe-alias `matricize` matricization directly.
 for f in (
-        :gram_eigh_full, :gram_eigh_full_with_pinv,
-        :sqrth_safe, :invsqrth_safe, :sqrth_invsqrth_safe,
+        :sqrth_safe, :invsqrth_safe,
     )
     @eval begin
         function $f(
@@ -90,8 +89,8 @@ for f in (
         :left_polar, :right_polar, :left_orth, :right_orth,
         :svd_compact, :svd_full, :svd_trunc, :svd_vals,
         :eigh_full, :eig_full, :eigh_trunc, :eig_trunc, :eigh_vals, :eig_vals,
-        :left_null, :right_null, :gram_eigh_full, :gram_eigh_full_with_pinv,
-        :sqrth_safe, :invsqrth_safe, :sqrth_invsqrth_safe, :project_hermitian,
+        :left_null, :right_null,
+        :sqrth_safe, :invsqrth_safe, :project_hermitian,
     )
     @eval begin
         function $f(style::MatricizeStyle, A, ndims_codomain::Val{K}; kwargs...) where {K}
@@ -612,124 +611,6 @@ function unmatricize_factors(
     return unmatricize(style, Nᴴ, (axes(Nᴴ, 1),), axes_domain)
 end
 
-"""
-    gram_eigh_full(A, labels_A, labels_codomain, labels_domain; kwargs...) -> X
-    gram_eigh_full(A, perm_codomain, perm_domain; kwargs...) -> X
-    gram_eigh_full(A, ndims_codomain::Val; kwargs...) -> X
-
-Gram factorization of a generic N-dimensional array, interpreting it as a
-Hermitian positive semi-definite linear map from the domain to the codomain
-dimensions. Returns `X` such that `A ≈ X * X'` (contracted on the rank leg),
-i.e. the codomain axes of `X` match the codomain axes of `A` and `X` has a
-single trailing rank axis.
-
-## Keyword arguments
-
-  - `alg`: forwarded to `MatrixAlgebraKit.eigh_full`.
-
-$(MatrixAlgebra._clamp_kwargs_doc("A"))
-
-# Examples
-
-```jldoctest
-julia> using TensorAlgebra: contract, gram_eigh_full
-
-julia> B = randn(3, 2, 2);
-
-julia> A = contract((:a, :b, :c, :d), conj(B), (:r, :a, :b), B, (:r, :c, :d));
-
-julia> X = gram_eigh_full(A, (:a, :b, :c, :d), (:a, :b), (:c, :d));
-
-julia> A ≈ contract((:a, :b, :c, :d), X, (:a, :b, :r), conj(X), (:c, :d, :r))
-true
-```
-
-See also [`gram_eigh_full_with_pinv`](@ref) and
-[`MatrixAlgebra.gram_eigh_full`](@ref).
-"""
-gram_eigh_full
-
-function gram_eigh_full!!(
-        style::MatricizeStyle, A, ndims_codomain::Val; kwargs...
-    )
-    A_mat = matricize(style, A, identitybiperm(A, ndims_codomain)...)
-    X = MatrixAlgebra.gram_eigh_full!!(A_mat; kwargs...)
-    axes_codomain = first(bipartition(axes(A), ndims_codomain))
-    return unmatricize(style, X, axes_codomain, (conj(axes(X, ndims(X))),))
-end
-function gram_eigh_full!!(A, ndims_codomain::Val; kwargs...)
-    return gram_eigh_full!!(MatricizeStyle(A), A, ndims_codomain; kwargs...)
-end
-
-function unmatricize_factors(
-        ::typeof(gram_eigh_full), style::MatricizeStyle, X,
-        axes_codomain, axes_domain
-    )
-    return unmatricize(style, X, axes_codomain, (conj(axes(X, ndims(X))),))
-end
-
-"""
-    gram_eigh_full_with_pinv(A, labels_A, labels_codomain, labels_domain; kwargs...) -> X, Y
-    gram_eigh_full_with_pinv(A, perm_codomain, perm_domain; kwargs...) -> X, Y
-    gram_eigh_full_with_pinv(A, ndims_codomain::Val; kwargs...) -> X, Y
-
-Like [`gram_eigh_full`](@ref), but additionally returns `Y ≈ pinv(X)` such
-that `Y * X ≈ I` on the rank subspace (a left inverse). The codomain axes
-of `X` match the codomain axes of `A`; `Y` has a leading rank axis followed
-by the codomain axes.
-
-## Keyword arguments
-
-  - `alg`: forwarded to `MatrixAlgebraKit.eigh_full`.
-
-$(MatrixAlgebra._clamp_kwargs_doc("A"))
-
-# Examples
-
-```jldoctest
-julia> using LinearAlgebra: I
-
-julia> using TensorAlgebra: contract, gram_eigh_full_with_pinv
-
-julia> B = randn(8, 2, 2);
-
-julia> A = contract((:a, :b, :c, :d), conj(B), (:r, :a, :b), B, (:r, :c, :d));
-
-julia> X, Y = gram_eigh_full_with_pinv(A, (:a, :b, :c, :d), (:a, :b), (:c, :d));
-
-julia> A ≈ contract((:a, :b, :c, :d), X, (:a, :b, :r), conj(X), (:c, :d, :r))
-true
-
-julia> contract((:r, :s), Y, (:r, :a, :b), X, (:a, :b, :s)) ≈ I
-true
-```
-
-See also [`MatrixAlgebra.gram_eigh_full_with_pinv`](@ref).
-"""
-gram_eigh_full_with_pinv
-
-function gram_eigh_full_with_pinv!!(
-        style::MatricizeStyle, A, ndims_codomain::Val; kwargs...
-    )
-    A_mat = matricize(style, A, identitybiperm(A, ndims_codomain)...)
-    X, Y = MatrixAlgebra.gram_eigh_full_with_pinv!!(A_mat; kwargs...)
-    axes_codomain = first(bipartition(axes(A), ndims_codomain))
-    return unmatricize(style, X, axes_codomain, (conj(axes(X, ndims(X))),)),
-        unmatricize(style, Y, (axes(Y, 1),), axes_codomain)
-end
-function gram_eigh_full_with_pinv!!(A, ndims_codomain::Val; kwargs...)
-    return gram_eigh_full_with_pinv!!(MatricizeStyle(A), A, ndims_codomain; kwargs...)
-end
-
-function unmatricize_factors(
-        ::typeof(gram_eigh_full_with_pinv), style::MatricizeStyle, F,
-        axes_codomain, axes_domain
-    )
-    X, Y = F
-    return unmatricize(style, X, axes_codomain, (conj(axes(X, ndims(X))),)),
-        unmatricize(style, Y, (axes(Y, 1),), axes_codomain)
-end
-
 """
     sqrth_safe(A, labels_A, labels_codomain, labels_domain; kwargs...) -> P
     sqrth_safe(A, perm_codomain, perm_domain; kwargs...) -> P
@@ -748,7 +629,7 @@ up to numerical noise.
 
 $(MatrixAlgebra._clamp_kwargs_doc("A"))
 
-See also [`invsqrth_safe`](@ref), [`sqrth_invsqrth_safe`](@ref), and
+See also [`invsqrth_safe`](@ref) and
 [`MatrixAlgebra.sqrth_safe`](@ref).
 """
 sqrth_safe
@@ -771,7 +652,7 @@ first if it is Hermitian only up to numerical noise.
 
 $(MatrixAlgebra._clamp_kwargs_doc("A"))
 
-See also [`sqrth_safe`](@ref), [`sqrth_invsqrth_safe`](@ref), and
+See also [`sqrth_safe`](@ref) and
 [`MatrixAlgebra.invsqrth_safe`](@ref).
 """
 invsqrth_safe
@@ -807,35 +688,6 @@ function unmatricize_factors(
     return unmatricize(style, H_mat, axes_codomain, axes_domain)
 end
 
-"""
-    sqrth_invsqrth_safe(A, labels_A, labels_codomain, labels_domain; kwargs...) -> P, Pinv
-    sqrth_invsqrth_safe(A, perm_codomain, perm_domain; kwargs...) -> P, Pinv
-    sqrth_invsqrth_safe(A, ndims_codomain::Val; kwargs...) -> P, Pinv
-
-Square root and pseudo-inverse square root of a generic N-dimensional
-array (see [`sqrth_safe`](@ref) and [`invsqrth_safe`](@ref)), from a
-single eigendecomposition. Both results carry the same codomain and
-domain axes as `A`.
-
-## Keyword arguments
-
-  - `alg`: forwarded to `MatrixAlgebraKit.eigh_full`.
-
-$(MatrixAlgebra._clamp_kwargs_doc("A"))
-
-See also [`MatrixAlgebra.sqrth_invsqrth_safe`](@ref).
-"""
-sqrth_invsqrth_safe
-
-function unmatricize_factors(
-        ::typeof(sqrth_invsqrth_safe), style::MatricizeStyle, F,
-        axes_codomain, axes_domain
-    )
-    P_mat, Pinv_mat = F
-    return unmatricize(style, P_mat, axes_codomain, axes_domain),
-        unmatricize(style, Pinv_mat, axes_codomain, axes_domain)
-end
-
 """
     TensorAlgebra.one(A, labels_A, labels_codomain, labels_domain) -> Id
     TensorAlgebra.one(A, perm_codomain, perm_domain) -> Id
diff --git a/test/test_exports.jl b/test/test_exports.jl
index 0ff6b33..5cb1aea 100644
--- a/test/test_exports.jl
+++ b/test/test_exports.jl
@@ -13,8 +13,6 @@ using Test: @test, @testset
         :eigh_full,
         :eigh_trunc,
         :eigh_vals,
-        :gram_eigh_full,
-        :gram_eigh_full_with_pinv,
         :invsqrth_safe,
         :isdual,
         :left_null,
@@ -28,7 +26,6 @@ using Test: @test, @testset
         :right_null,
         :right_orth,
         :right_polar,
-        :sqrth_invsqrth_safe,
         :sqrth_safe,
         :svd_compact,
         :svd_full,
@@ -56,15 +53,12 @@ using Test: @test, @testset
 
     exports = [
         :MatrixAlgebra,
-        :gram_eigh_full,
-        :gram_eigh_full_with_pinv,
         :invsqrt_diag_safe,
         :invsqrth_safe,
         :pow_diag_safe,
         :pow_diag_safe!,
         :powh_safe,
         :sqrt_diag_safe,
-        :sqrth_invsqrth_safe,
         :sqrth_safe,
     ]
     @test issetequal(names(TensorAlgebra.MatrixAlgebra), exports)
diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl
index d69a007..d9f62c6 100644
--- a/test/test_factorizations.jl
+++ b/test/test_factorizations.jl
@@ -1,9 +1,8 @@
 using LinearAlgebra: LinearAlgebra, Diagonal, I, diag, norm
 using MatrixAlgebraKit: truncrank
 using TensorAlgebra: TensorAlgebra, contract, eig_full, eig_vals, eigh_full, eigh_vals,
-    gram_eigh_full, gram_eigh_full_with_pinv, left_null, left_orth, left_polar, lq_compact,
-    lq_full, qr_compact, qr_full, right_null, right_orth, right_polar, svd_compact,
-    svd_full, svd_trunc, svd_vals
+    left_null, left_orth, left_polar, lq_compact, lq_full, qr_compact, qr_full, right_null,
+    right_orth, right_polar, svd_compact, svd_full, svd_trunc, svd_vals
 using Test: @test, @testset
 using TestExtras: @constinferred
 
@@ -306,59 +305,6 @@ end
     @test A ≈ contract(labels_A, P, (labels_P..., :w), W, (:w, labels_W...))
 end
 
-# Gram factorization
-# ------------------
-# Build a Hermitian positive semi-definite tensor A[a,b,c,d] with codomain
-# (a, b) and domain (c, d): pick a random B[k, a, b] (k = aux), then form
-# A = B' * B over k. By construction A ≈ X' * X for X[r, a, b] with rank r
-# bounded by k (rank leg first, following the Cholesky `A = U' * U`
-# convention).
-@testset "Full-rank gram_eigh_full ($T)" for T in elts
-    B = randn(T, 6, 2, 3) # k = 6, codomain = (a, b) of size 2*3 = 6 -> full rank
-    A = contract((:a, :b, :c, :d), conj(B), (:k, :a, :b), B, (:k, :c, :d))
-    labels_A = (:a, :b, :c, :d)
-    labels_X = (:a, :b)
-    labels_Y = (:c, :d)
-
-    Acopy = copy(A)
-    X = @constinferred gram_eigh_full(A, labels_A, labels_X, labels_Y)
-    @test A == Acopy # should not have altered initial array
-    A′ = contract(labels_A, X, (:a, :b, :r), conj(X), (:c, :d, :r))
-    @test A ≈ A′
-    @test size(X, ndims(X)) == size(A, 1) * size(A, 2)
-
-    # `Val`, perm, and label entries agree.
-    @test gram_eigh_full(A, Val(2)) ≈ X
-    @test gram_eigh_full(A, (1, 2), (3, 4)) ≈ X
-
-    # `with_pinv` variant: Y is a left inverse of X (Y * X ≈ I on the
-    # rank subspace).
-    X2, Y2 = @constinferred gram_eigh_full_with_pinv(A, labels_A, labels_X, labels_Y)
-    @test A ≈ contract(labels_A, X2, (:a, :b, :r), conj(X2), (:c, :d, :r))
-    YX = contract((:r, :s), Y2, (:r, :a, :b), X2, (:a, :b, :s))
-    @test YX ≈ I
-end
-
-@testset "Rank-deficient gram_eigh_full ($T)" for T in elts
-    B = randn(T, 4, 2, 3) # k = 4 < codomain dim 6, so A is rank-4
-    A = contract((:a, :b, :c, :d), conj(B), (:k, :a, :b), B, (:k, :c, :d))
-
-    # Recovery of A is independent of the `rtol` cutoff because all
-    # nonzero eigenvalues sit far above any reasonable threshold.
-    X = gram_eigh_full(A, Val(2); rtol = 1.0e-10)
-    @test A ≈ contract(
-        (:a, :b, :c, :d), X, (:a, :b, :r), conj(X), (:c, :d, :r)
-    )
-
-    # Moore–Penrose-like identity: X * Y * X ≈ X when Y is pinv(X). With
-    # cod-first X and rank-first Y, contract Y[r, a, b] * X[a, b, s] → P[r, s]
-    # (projector onto the rank subspace), then X * P → X.
-    X2, Y2 = gram_eigh_full_with_pinv(A, Val(2); rtol = 1.0e-10)
-    P = contract((:r, :s), Y2, (:r, :a, :b), X2, (:a, :b, :s))
-    XP = contract((:c, :d, :r), X2, (:c, :d, :s), P, (:s, :r))
-    @test XP ≈ X2
-end
-
 # one (identity tensor)
 # ---------------------
 # An identity tensor matricized along its codomain/domain partition is the
diff --git a/test/test_matrixalgebra.jl b/test/test_matrixalgebra.jl
index e279ea4..bda1a27 100644
--- a/test/test_matrixalgebra.jl
+++ b/test/test_matrixalgebra.jl
@@ -147,42 +147,6 @@ elts = (Float32, Float64, ComplexF32, ComplexF64)
         @test norm(ũ * s̃ * ṽ) ≈ 0
     end
 
-    @testset "gram_eigh_full" begin
-        n = 5
-        # Full-rank Hermitian PSD. Use a tall random factor so `B' * B`
-        # is comfortably full rank even at Float32 precision (a square
-        # random `B` can produce a `B' * B` whose smallest eigenvalue
-        # falls below the default rtol clamp on some seeds).
-        rng = StableRNG(123)
-        B = randn(rng, elt, 2n, n)
-        A = B' * B
-        X = MatrixAlgebra.gram_eigh_full(A)
-        @test X * X' ≈ A
-        @test size(X) == (n, n)
-
-        X2, Y2 = MatrixAlgebra.gram_eigh_full_with_pinv(A)
-        @test X2 * X2' ≈ A
-        @test Y2 * X2 ≈ I(n)
-
-        # `!!` variant accepts a destroyable copy.
-        Xb = MatrixAlgebra.gram_eigh_full!!(copy(A))
-        @test Xb * Xb' ≈ A
-
-        # Rank deficient: A is n×n of rank k < n. Recovery of A still holds;
-        # X * Y is the projector onto the rank-k codomain subspace
-        # (idempotent, rank k), and X * P ≈ X (Moore–Penrose).
-        k = 3
-        Brd = randn(rng, elt, k, n)
-        Ard = Brd' * Brd
-        Xrd, Yrd = MatrixAlgebra.gram_eigh_full_with_pinv(
-            Ard; rtol = sqrt(eps(real(elt)))
-        )
-        @test Xrd * Xrd' ≈ Ard
-        P = Xrd * Yrd
-        @test P * P ≈ P
-        @test P * Xrd ≈ Xrd
-    end
-
     @testset "powh_safe / sqrth_safe / invsqrth_safe" begin
         n = 4
         rng = StableRNG(123)

From ec471e83f2224da38c5d73feab8cdffd5ce6c7dd Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 15 Sep 2026 22:56:31 -0400
Subject: [PATCH 14/17] Let a style overload the total matricize copy

A graded array already gets owned storage out of `permutedimsop`, whose stored matrix is the answer, so decomposing the copy into an allocation plus an in-place write would copy that storage a second time.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/matricize.jl | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/matricize.jl b/src/matricize.jl
index 23cff0f..830872e 100644
--- a/src/matricize.jl
+++ b/src/matricize.jl
@@ -86,6 +86,12 @@ end
 # needs the fused axes, which only the style knows: TensorAlgebra deliberately has no generic
 # axis-fusion interface. It is also what makes the copy path terminate, since `matricizeop!` is a
 # distinct function from the router rather than a re-entry into it.
+#
+# `matricizeopcopy` is itself an overload point for a style whose owned matricization already falls
+# out of an allocating operation it has (for a graded array, permuting into fresh storage whose
+# stored matrix is the answer). Such a style overloads the copy instead of
+# `allocate_output`/`matricizeop!`, which would copy that storage a second time, and then owes
+# only `matricizeopview` and `is_output_view`.
 
 """
     matricizeop(op, a, perm_codomain, perm_domain)

From 4a830205e125fb9ffbf94f2bf3f9dedbc1e52abc Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 16 Sep 2026 00:17:59 -0400
Subject: [PATCH 15/17] Give the contract interface one name per argument form

The labels forms keep the plain names and the bipermutation forms take a `perm` marker, freeing `contract` to become variadic over operands. `allocate_contract_output` is gone in favor of overloading `allocate_output`, and a generic `select_algorithm` sits above the per-operation resolvers.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 ext/TensorAlgebraTensorKitExt.jl              |   2 +-
 .../TensorAlgebraTensorOperationsExt.jl       |   8 +-
 src/TensorAlgebra.jl                          |   6 +-
 src/algorithm.jl                              |  39 +++++
 src/contract/allocate_output.jl               |  18 +--
 src/contract/contract.jl                      | 134 ++++++++++++++----
 src/contract/contract_matricize.jl            |   2 +-
 src/diagonal.jl                               |  34 +++--
 src/factorizations.jl                         |   6 +-
 test/test_basics.jl                           |  41 +++---
 test/test_exports.jl                          |  10 +-
 test/test_factorizations.jl                   |  82 ++++++-----
 test/test_matricize.jl                        |   2 +-
 test/test_mooncakeext.jl                      |   6 +-
 14 files changed, 279 insertions(+), 111 deletions(-)
 create mode 100644 src/algorithm.jl

diff --git a/ext/TensorAlgebraTensorKitExt.jl b/ext/TensorAlgebraTensorKitExt.jl
index 58fa709..90bbb92 100644
--- a/ext/TensorAlgebraTensorKitExt.jl
+++ b/ext/TensorAlgebraTensorKitExt.jl
@@ -293,7 +293,7 @@ end
 # already implements through its TensorOperations interface. Route the generic `contract`
 # there: `zero!` clears the `similar_map`-allocated destination, and the default algorithm
 # hands the in-place contraction to the TensorOperations backend (see the TensorOperations
-# extension's `contractopadd!`).
+# extension's `contractpermopadd!`).
 TensorAlgebra.zero!(t::AbstractTensorMap) = VectorInterface.zerovector!(t)
 
 # A `TensorMap` is not an `AbstractArray`, so the generic in-place `TensorAlgebra` operations
diff --git a/ext/TensorAlgebraTensorOperationsExt/TensorAlgebraTensorOperationsExt.jl b/ext/TensorAlgebraTensorOperationsExt/TensorAlgebraTensorOperationsExt.jl
index 2b62f80..79c84f3 100644
--- a/ext/TensorAlgebraTensorOperationsExt/TensorAlgebraTensorOperationsExt.jl
+++ b/ext/TensorAlgebraTensorOperationsExt/TensorAlgebraTensorOperationsExt.jl
@@ -22,7 +22,7 @@ end
 # ----------------------------------------------------------------
 
 # not in-place
-function TA.contract(
+function TA.contractperm(
         algorithm::TensorOperationsAlgorithm,
         perm_dest_codomain, perm_dest_domain,
         a1::AbstractArray, perm1_codomain, perm1_domain,
@@ -39,7 +39,7 @@ function TA.contract(
     )
 end
 
-function TA.contract(
+function TA.contractalign(
         algorithm::TensorOperationsAlgorithm,
         labels_dest,
         a1::AbstractArray, labels1,
@@ -56,7 +56,7 @@ function TA.contract(
 end
 
 # in-place
-function TA.contractopadd!(
+function TA.contractpermopadd!(
         algorithm::TensorOperationsAlgorithm,
         a_dest, perm_dest_codomain, perm_dest_domain,
         op1, a1, perm1_codomain, perm1_domain,
@@ -89,7 +89,7 @@ function TO.tensorcontract!(
     )
     op1 = conj1 ? conj : identity
     op2 = conj2 ? conj : identity
-    return TA.contractopadd!(
+    return TA.contractpermopadd!(
         backend,
         a_dest, permblocks_dest...,
         op1, a1, permblocks1...,
diff --git a/src/TensorAlgebra.jl b/src/TensorAlgebra.jl
index 6ae6d18..b3b04fd 100644
--- a/src/TensorAlgebra.jl
+++ b/src/TensorAlgebra.jl
@@ -1,6 +1,7 @@
 module TensorAlgebra
 
-export contract, contract!, dual, eig_full, eig_trunc, eig_vals, eigh_full, eigh_trunc,
+export contract, contract!, contractalign, dual, eig_full, eig_trunc, eig_vals, eigh_full,
+    eigh_trunc,
     eigh_vals, invsqrth_safe, isdual, left_null,
     left_orth, left_polar, lq_compact, lq_full, project_hermitian, qr_compact,
     qr_full, right_null, right_orth, right_polar, sqrth_safe,
@@ -9,7 +10,7 @@ export contract, contract!, dual, eig_full, eig_trunc, eig_vals, eigh_full, eigh
 if VERSION >= v"1.11.0-DEV.469"
     eval(
         Meta.parse(
-            "public biperm, bipartition, cat_similar, concatenate, concatenate!, ContractAlgorithm, contractopadd!, data, datatype, directsum, flattenlinear, is_output_view, label_type, matricize, matricizeop, matricizeop!, matricizeopcopy, matricizeopview, permutedims, permutedims!, scalar, similar_map, TensorOperationsAlgorithm, to_range, tr, tryflattenlinear, ungrade, zero!, scale!, permuteddims, PermutedDims"
+            "public allocate_output, biperm, bipartition, cat_similar, check_input, concatenate, concatenate!, ContractAlgorithm, contractopadd!, contractperm, contractperm!, contractpermadd!, contractpermopadd!, data, datatype, directsum, flattenlinear, is_output_view, label_type, matricize, matricizeop, matricizeop!, matricizeopcopy, matricizeopview, output_axes, select_algorithm, default_algorithm, permutedims, permutedims!, scalar, similar_map, TensorOperationsAlgorithm, to_range, tr, tryflattenlinear, ungrade, zero!, scale!, permuteddims, PermutedDims"
         )
     )
 end
@@ -25,6 +26,7 @@ include("concatenate.jl")
 include("directsum.jl")
 include("dual.jl")
 include("to_range.jl")
+include("algorithm.jl")
 include("contract/contractalgorithm.jl")
 include("contract/contract.jl")
 include("contract/contract_labels.jl")
diff --git a/src/algorithm.jl b/src/algorithm.jl
new file mode 100644
index 0000000..28775bb
--- /dev/null
+++ b/src/algorithm.jl
@@ -0,0 +1,39 @@
+"""
+    TensorAlgebra.select_algorithm(f, args...; alg = nothing, kwargs...)
+
+Resolve the algorithm operation `f` should run with on `args`. An `alg` of `nothing` defers to
+[`TensorAlgebra.default_algorithm`](@ref); anything else is validated and passed through.
+
+This is the forward-facing layer the per-operation resolvers sit under, so a caller that does not
+care which operation it is dispatching writes `select_algorithm(f, ...)` and a backend registers
+its choice on `default_algorithm(f, ...)`.
+"""
+function select_algorithm(f, args...; alg = nothing, kwargs...)
+    isnothing(alg) && return default_algorithm(f, args...; kwargs...)
+    return select_algorithm_specified(f, alg, args...; kwargs...)
+end
+
+"""
+    TensorAlgebra.default_algorithm(f, args...)
+    TensorAlgebra.default_algorithm(f, argtypes::Type...)
+
+The algorithm operation `f` runs with on `args` when the caller names none. The types form is the
+registration point for a storage type; the values form defaults to it.
+
+Each operation bridges to its own resolver, so `default_algorithm(contract, A1, A2)` is
+[`TensorAlgebra.default_contract_algorithm`](@ref).
+"""
+function default_algorithm(f, args...; kwargs...)
+    return default_algorithm(f, map(typeof, args)...; kwargs...)
+end
+function default_algorithm(f, argtypes::Type...; kwargs...)
+    return throw(MethodError(default_algorithm, (f, argtypes...)))
+end
+
+# `alg` named something. A resolved algorithm object passes through; anything else is a caller
+# error, reported against the operation rather than as a `MethodError` from inside the resolver.
+function select_algorithm_specified(f, alg, args...; kwargs...)
+    return throw(
+        ArgumentError("`$alg` is not an algorithm for `$f`")
+    )
+end
diff --git a/src/contract/allocate_output.jl b/src/contract/allocate_output.jl
index 13b392c..5fb3210 100644
--- a/src/contract/allocate_output.jl
+++ b/src/contract/allocate_output.jl
@@ -68,13 +68,16 @@ function output_axes(
         axes_uncontracted, perm_dest_codomain, perm_dest_domain
     )
     # The operand axes are stored/dualized, so un-dualize the domain axes into the codomain-facing
-    # construction convention shared by `allocate_contract_output`, `similar_map`, and `unmatricize`
-    # (a no-op on dense axes).
+    # construction convention shared by `similar_map` and `unmatricize` (a no-op on dense axes).
     return axes_codomain_dest, conj.(axes_domain_dest)
 end
 
 # TODO: Use `ArrayLayouts`-like `MulAdd` object,
 # i.e. `ContractAdd`?
+# The destination `contract` writes into. A structured operand type overloads this directly, deriving
+# the axes and element type from `output_axes` and `Base.promote_op` as below; the permutations are
+# part of the signature because the contraction pattern is not recoverable from the destination leg
+# counts alone.
 function allocate_output(
         ::typeof(contract),
         perm_dest_codomain, perm_dest_domain,
@@ -97,14 +100,5 @@ function allocate_output(
         a2, perm2_codomain, perm2_domain
     )
     T = Base.promote_op(matprod, eltype(a1), eltype(a2))
-    return allocate_contract_output(a1, a2, T, axes_codomain_dest, axes_domain_dest)
-end
-
-# Allocate the output container for `contract`: the operand types, the output element type and
-# axes (domain codomain-facing), and the output's codomain/domain leg counts (the axes tuple
-# lengths) select the container type. Internal to TensorAlgebra, not a public extension point:
-# the leg counts identify the contraction pattern only for matrix-shaped operands (see the
-# `Diagonal` method in `diagonal.jl`), so external structured types should not overload it.
-function allocate_contract_output(a1, a2, T, axes_codomain::Tuple, axes_domain::Tuple)
-    return zero!(similar_map(a1, T, axes_codomain, axes_domain))
+    return zero!(similar_map(a1, T, axes_codomain_dest, axes_domain_dest))
 end
diff --git a/src/contract/contract.jl b/src/contract/contract.jl
index f31f14e..a3d47e1 100644
--- a/src/contract/contract.jl
+++ b/src/contract/contract.jl
@@ -2,15 +2,65 @@
 # TODO: Add `scaledcontract(a1, labels1, a2, labels2, α) = α * contract(a1, labels1, a2, labels2)`.
 
 # contract (labels)
+"""
+    contract(a1, labels1, a2, labels2, ...; alg = nothing) -> a_dest, labels_dest
+
+Contract the arrays over the labels they share, returning the result along with the labels of its
+dimensions. A label appearing on two operands is summed over, one appearing on a single operand
+survives, and `labels_dest` reports the surviving labels in the order the result carries them.
+
+Operands past the second are contracted one pair at a time from left to right, so the call
+expresses the contraction order rather than requesting an optimized one.
+
+```jldoctest
+julia> using TensorAlgebra: contract
+
+julia> a, b = randn(2, 3), randn(3, 4);
+
+julia> ab, labels = contract(a, (:i, :j), b, (:j, :k));
+
+julia> labels
+(:i, :k)
+
+julia> ab ≈ a * b
+true
+```
+
+See also [`contractalign`](@ref) to name the output dimensions and get back just the array.
+"""
 function contract(a1, labels1, a2, labels2; kwargs...)
     # Optionally convert the labels to a representation cheaper to run the bookkeeping on (see
     # `label_type`). `encode_contraction_labels`/`decode_contraction_labels` are no-ops unless the label type opts in.
     l1, l2 = encode_contraction_labels(labels1, labels2)
     l_dest = contract_labels(l1, l2)
-    a_dest = contract(l_dest, a1, l1, a2, l2; kwargs...)
+    a_dest = contractalign(l_dest, a1, l1, a2, l2; kwargs...)
     return a_dest, decode_contraction_labels(l_dest, labels1, labels2)
 end
-function contract(
+function contract(a1, labels1, a2, labels2, a3, labels3, rest...; kwargs...)
+    check_alternating_labels(contract, rest)
+    a12, labels12 = contract(a1, labels1, a2, labels2; kwargs...)
+    return contract(a12, labels12, a3, labels3, rest...; kwargs...)
+end
+
+"""
+    contractalign(labels_dest, a1, labels1, a2, labels2, ...; alg = nothing) -> a_dest
+
+Contract the arrays over the labels they share into a result whose dimensions carry
+`labels_dest`, which must be the surviving labels in some order.
+
+This is [`contract`](@ref) with the output specified, so it returns the array on its own. The
+name matches `ITensorBase.align`: arrange the result's dimensions to match the labels given.
+
+```jldoctest
+julia> using TensorAlgebra: contractalign
+
+julia> a, b = randn(2, 3), randn(3, 4);
+
+julia> contractalign((:k, :i), a, (:i, :j), b, (:j, :k)) ≈ permutedims(a * b, (2, 1))
+true
+```
+"""
+function contractalign(
         labels_dest, a1, labels1, a2, labels2; kwargs...
     )
     t1 = ntuple(i -> labels1[i], Val(ndims(a1)))
@@ -18,7 +68,7 @@ function contract(
     contracted1 = map(in(t2), t1)
     # Cross into a `Val(K)` method (a function-barrier on the contracted count) so the
     # bipartitioned permutations and the contraction below them are type-stable.
-    return _contract(
+    return _contractalign(
         Val(count(contracted1)),
         labels_dest,
         a1,
@@ -29,17 +79,38 @@ function contract(
         kwargs...
     )
 end
-function _contract(
+# Only the last pair lands on the requested output labels; the ones before it infer their own.
+function contractalign(
+        labels_dest, a1, labels1, a2, labels2, a3, labels3, rest...; kwargs...
+    )
+    check_alternating_labels(contractalign, rest)
+    a12, labels12 = contract(a1, labels1, a2, labels2; kwargs...)
+    return contractalign(labels_dest, a12, labels12, a3, labels3, rest...; kwargs...)
+end
+function _contractalign(
         ::Val{K}, labels_dest, a1, labels1, a2, labels2,
         contracted1; kwargs...
     ) where {K}
     biperm_dest, biperm1, biperm2 =
         biperms(contract, Val(K), labels_dest, labels1, labels2, contracted1)
-    return contract(biperm_dest..., a1, biperm1..., a2, biperm2...; kwargs...)
+    return contractperm(biperm_dest..., a1, biperm1..., a2, biperm2...; kwargs...)
 end
 
-# contract (bipartitioned permutations)
-function contract(
+# The variadic forms take arrays and labels in alternating positions, so a trailing group with an
+# odd length is a miscount at the call site rather than something to diagnose further down.
+function check_alternating_labels(f, rest::Tuple)
+    iseven(length(rest)) || throw(
+        ArgumentError(
+            "`$f` takes each array followed by its labels, so the trailing arguments must come in pairs"
+        )
+    )
+    return nothing
+end
+
+# contractperm (bipartitioned permutations)
+# `perm` marks the whole biperm ladder: every rung has a labels-form sibling under the plain name,
+# and once `contract` is variadic over operands the two can no longer be told apart by arity.
+function contractperm(
         a1, perm1_codomain, perm1_domain,
         a2, perm2_codomain, perm2_domain;
         kwargs...
@@ -48,14 +119,14 @@ function contract(
     Ndest = Val(length(perm1_codomain) + length(perm2_domain))
     perm_dest_codomain, perm_dest_domain =
         bipartition(ntuple(identity, Ndest), Ndest_codomain)
-    return contract(
+    return contractperm(
         perm_dest_codomain, perm_dest_domain,
         a1, perm1_codomain, perm1_domain,
         a2, perm2_codomain, perm2_domain;
         kwargs...
     )
 end
-function contract(
+function contractperm(
         perm_dest_codomain, perm_dest_domain,
         a1, perm1_codomain, perm1_domain,
         a2, perm2_codomain, perm2_domain;
@@ -67,7 +138,7 @@ function contract(
         a1, perm1_codomain, perm1_domain,
         a2, perm2_codomain, perm2_domain
     )
-    return contract!(
+    return contractperm!(
         a_dest, perm_dest_codomain, perm_dest_domain,
         a1, perm1_codomain, perm1_domain,
         a2, perm2_codomain, perm2_domain;
@@ -86,13 +157,13 @@ function contract!(
         a_dest, labels_dest, a1, labels1, a2, labels2, true, false; kwargs...
     )
 end
-function contract!(
+function contractperm!(
         a_dest, perm_dest_codomain, perm_dest_domain,
         a1, perm1_codomain, perm1_domain,
         a2, perm2_codomain, perm2_domain;
         kwargs...
     )
-    return contractadd!(
+    return contractpermadd!(
         a_dest, perm_dest_codomain, perm_dest_domain,
         a1, perm1_codomain, perm1_domain,
         a2, perm2_codomain, perm2_domain,
@@ -112,15 +183,15 @@ function contractadd!(
         a_dest, labels_dest, identity, a1, labels1, identity, a2, labels2, α, β; kwargs...
     )
 end
-# contractadd! (bipartitioned permutations)
-function contractadd!(
+# contractpermadd! (bipartitioned permutations)
+function contractpermadd!(
         a_dest, perm_dest_codomain, perm_dest_domain,
         a1, perm1_codomain, perm1_domain,
         a2, perm2_codomain, perm2_domain,
         α::Number, β::Number;
         kwargs...
     )
-    return contractopadd!(
+    return contractpermopadd!(
         a_dest, perm_dest_codomain, perm_dest_domain,
         identity, a1, perm1_codomain, perm1_domain,
         identity, a2, perm2_codomain, perm2_domain,
@@ -153,17 +224,17 @@ function _contractopadd!(
     ) where {K}
     biperm_dest, biperm1, biperm2 =
         biperms(contract, Val(K), labels_dest, labels1, labels2, contracted1)
-    return contractopadd!(
+    return contractpermopadd!(
         a_dest, biperm_dest..., op1, a1, biperm1..., op2, a2, biperm2..., α, β; kwargs...
     )
 end
-# contractopadd! (bipartitioned permutations, algorithm selection)
-function contractopadd!(
+# contractpermopadd! (bipartitioned permutations, algorithm selection)
+function contractpermopadd!(
         a_dest, perm_dest_codomain, perm_dest_domain,
         op1, a1, perm1_codomain, perm1_domain,
         op2, a2, perm2_codomain, perm2_domain,
         α::Number, β::Number;
-        alg = DefaultContractAlgorithm(), kwargs...
+        alg = nothing, kwargs...
     )
     check_input(
         contract!,
@@ -171,8 +242,8 @@ function contractopadd!(
         a1, perm1_codomain, perm1_domain,
         a2, perm2_codomain, perm2_domain
     )
-    algorithm = select_contract_algorithm(alg, a1, a2; kwargs...)
-    return contractopadd!(
+    algorithm = select_algorithm(contract, a1, a2; alg, kwargs...)
+    return contractpermopadd!(
         algorithm,
         a_dest, perm_dest_codomain, perm_dest_domain,
         op1, a1, perm1_codomain, perm1_domain,
@@ -180,9 +251,9 @@ function contractopadd!(
         α, β
     )
 end
-# contractopadd! (dispatched on the algorithm, bipartitioned permutations)
+# contractpermopadd! (dispatched on the algorithm, bipartitioned permutations)
 # Required interface if not using matricized contraction
-function contractopadd!(
+function contractpermopadd!(
         algorithm::ContractAlgorithm,
         a_dest, perm_dest_codomain, perm_dest_domain,
         op1, a1, perm1_codomain, perm1_domain,
@@ -191,7 +262,7 @@ function contractopadd!(
     )
     return throw(
         MethodError(
-            contractopadd!,
+            contractpermopadd!,
             (
                 algorithm,
                 a_dest, perm_dest_codomain, perm_dest_domain,
@@ -202,3 +273,18 @@ function contractopadd!(
         )
     )
 end
+
+# Bridges from the operation-generic algorithm layer in `algorithm.jl` down to the contraction
+# resolvers. They live here rather than beside those resolvers because dispatching on
+# `::typeof(contract)` needs `contract` to exist, and `contractalgorithm.jl` is included first for
+# the algorithm types this file's signatures use.
+function default_algorithm(::typeof(contract), A1::Type, A2::Type; kwargs...)
+    algorithm = default_contract_algorithm(A1, A2)
+    reject_algorithm_kwargs(algorithm; kwargs...)
+    return algorithm
+end
+function select_algorithm_specified(
+        ::typeof(contract), alg::ContractAlgorithm, a1, a2; kwargs...
+    )
+    return select_contract_algorithm(alg, a1, a2; kwargs...)
+end
diff --git a/src/contract/contract_matricize.jl b/src/contract/contract_matricize.jl
index 0b5fc95..409d032 100644
--- a/src/contract/contract_matricize.jl
+++ b/src/contract/contract_matricize.jl
@@ -1,6 +1,6 @@
 using LinearAlgebra: mul!
 
-function contractopadd!(
+function contractpermopadd!(
         algorithm::Matricize,
         a_dest::AbstractArray, biperm_dest_codomain, biperm_dest_domain,
         op1, a1::AbstractArray, biperm1_codomain, biperm1_domain,
diff --git a/src/diagonal.jl b/src/diagonal.jl
index 0488eff..d4befc6 100644
--- a/src/diagonal.jl
+++ b/src/diagonal.jl
@@ -69,13 +69,31 @@ function unmatricize(
     return unmatricize(style, copyto!(similar(m, axes(m)), m), axes_codomain, axes_domain)
 end
 
-# Contracting two `Diagonal`s over a single leg is the matmul/endomorphism pattern
-# `Diagonal * Diagonal = Diagonal` (all transpose variants `[i,j]*[j,k]`, `[i,j]*[k,j]`, ...),
-# whose `{1,1}` output stays `Diagonal`, so allocate one. Every other output shape (rank-4 outer
-# product, scalar full contraction) is not representable as a `Diagonal` and falls back to the
-# generic dense allocation, matching `Diagonal`/dense mixing.
-function allocate_contract_output(
-        a1::Diagonal, a2::Diagonal, T, axes_codomain::Tuple{Any}, axes_domain::Tuple{Any}
+function allocate_output(
+        ::typeof(contract),
+        perm_dest_codomain, perm_dest_domain,
+        a1::Diagonal, perm1_codomain, perm1_domain,
+        a2::Diagonal, perm2_codomain, perm2_domain
+    )
+    check_input(
+        contract, a1, perm1_codomain, perm1_domain, a2, perm2_codomain, perm2_domain
+    )
+    axes_codomain_dest, axes_domain_dest = output_axes(
+        contract,
+        perm_dest_codomain, perm_dest_domain,
+        a1, perm1_codomain, perm1_domain,
+        a2, perm2_codomain, perm2_domain
     )
-    return Diagonal(zero!(similar(a1.diag, T, (only(axes_codomain),))))
+    T = Base.promote_op(matprod, eltype(a1), eltype(a2))
+    # Contracting two `Diagonal`s over a single leg, leaving one free leg on each, is the
+    # matmul/endomorphism pattern `Diagonal * Diagonal = Diagonal` (all transpose variants
+    # `[i,j]*[j,k]`, `[i,j]*[k,j]`, ...), whose `{1,1}` output stays `Diagonal`. Every other
+    # pattern (rank-4 outer product, scalar full contraction) is not representable as a
+    # `Diagonal` and takes the generic dense allocation, matching `Diagonal`/dense mixing.
+    is_matmul =
+        length(perm1_codomain) == 1 && length(perm1_domain) == 1 &&
+        length(perm2_domain) == 1
+    is_matmul ||
+        return zero!(similar_map(a1, T, axes_codomain_dest, axes_domain_dest))
+    return Diagonal(zero!(similar(a1.diag, T, (only(axes_codomain_dest),))))
 end
diff --git a/src/factorizations.jl b/src/factorizations.jl
index 87fbad5..6344071 100644
--- a/src/factorizations.jl
+++ b/src/factorizations.jl
@@ -429,15 +429,15 @@ truncation error `ϵ`, the 2-norm of the discarded singular values.
 # Examples
 
 ```jldoctest
-julia> using TensorAlgebra: svd_trunc, contract
+julia> using TensorAlgebra: svd_trunc, contractalign
 
 julia> A = randn(4, 4);
 
 julia> U, S, Vᴴ, ϵ = svd_trunc(A, (:i, :j), (:i,), (:j,));
 
-julia> SV = contract((:u, :j), S, (:u, :v), Vᴴ, (:v, :j));
+julia> SV = contractalign((:u, :j), S, (:u, :v), Vᴴ, (:v, :j));
 
-julia> contract((:i, :j), U, (:i, :u), SV, (:u, :j)) ≈ A
+julia> contractalign((:i, :j), U, (:i, :u), SV, (:u, :j)) ≈ A
 true
 
 julia> isapprox(ϵ, 0; atol = 1e-10)
diff --git a/test/test_basics.jl b/test/test_basics.jl
index b0a119c..2ffabf8 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -1,8 +1,8 @@
 import TensorAlgebra
 using StableRNGs: StableRNG
 using TensorAlgebra: BiTuple, ContractAlgorithm, bipermutedims, bipermutedims!, contract,
-    contract!, contractadd!, length_codomain, length_domain, matricize, unmatricize,
-    unmatricize!
+    contract!, contractadd!, contractalign, length_codomain, length_domain, matricize,
+    unmatricize, unmatricize!
 using TensorOperations: TensorOperations
 using Test: @test, @test_broken, @test_throws, @testset
 
@@ -172,15 +172,22 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         a2 = randn(3, 4)
         # A keyword no algorithm can consume must name itself, not surface as a `MethodError`
         # from inside the resolver.
-        @test_throws ArgumentError contract((1, 3), a1, (1, 2), a2, (2, 3); nonsense = 1)
-        @test_throws ArgumentError contract(
+        @test_throws ArgumentError contractalign(
+            (1, 3),
+            a1,
+            (1, 2),
+            a2,
+            (2, 3);
+            nonsense = 1
+        )
+        @test_throws ArgumentError contractalign(
             (1, 3), a1, (1, 2), a2, (2, 3); alg = TensorAlgebra.Matricize(), nonsense = 1
         )
         # A non-algorithm passed as `alg` says so rather than erroring with "Not implemented".
         @test_throws ArgumentError TensorAlgebra.select_contract_algorithm(:nope, a1, a2)
         # The supported spellings still work.
-        @test contract((1, 3), a1, (1, 2), a2, (2, 3)) ≈ a1 * a2
-        @test contract(
+        @test contractalign((1, 3), a1, (1, 2), a2, (2, 3)) ≈ a1 * a2
+        @test contractalign(
             (1, 3), a1, (1, 2), a2, (2, 3); alg = TensorAlgebra.Matricize()
         ) ≈ a1 * a2
     end
@@ -201,8 +208,8 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         a_dest = ones(elt_dest, (1, 1))
         @test_throws ArgumentError contract(a1, (1, 2, 4), a2, (2, 3))
         @test_throws ArgumentError contract(a1, (1, 2), a2, (2, 3, 4))
-        @test_throws ArgumentError contract((1, 3, 4), a1, (1, 2), a2, (2, 3))
-        @test_throws ArgumentError contract((1, 3), a1, (1, 2), a2, (2, 4))
+        @test_throws ArgumentError contractalign((1, 3, 4), a1, (1, 2), a2, (2, 3))
+        @test_throws ArgumentError contractalign((1, 3), a1, (1, 2), a2, (2, 4))
         @test_throws ArgumentError contract!(a_dest, (1, 3, 4), a1, (1, 2), a2, (2, 3))
 
         dims = (2, 3, 4, 5, 6, 7, 8, 9, 10)
@@ -241,14 +248,14 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
             @test a_dest ≈ a_dest_tensoroperations
 
             # Specify destination labels
-            a_dest = contract(labels_dest, a1, labels1, a2, labels2)
-            a_dest_tensoroperations = contract(
+            a_dest = contractalign(labels_dest, a1, labels1, a2, labels2)
+            a_dest_tensoroperations = contractalign(
                 labels_dest, a1, labels1, a2, labels2; alg = alg_tensoroperations
             )
             @test a_dest ≈ a_dest_tensoroperations
 
-            a_dest = contract(labels_dest′, a1, labels1, a2, labels2)
-            a_dest_tensoroperations = contract(
+            a_dest = contractalign(labels_dest′, a1, labels1, a2, labels2)
+            a_dest_tensoroperations = contractalign(
                 labels_dest′, a1, labels1, a2, labels2; alg = alg_tensoroperations
             )
             @test a_dest ≈ a_dest_tensoroperations
@@ -286,7 +293,7 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         @test labels_dest == [L(1), L(4)]
 
         # Specifying the destination labels still works for opted-in types.
-        a_dest = contract([L(1), L(4)], a1, (L(1), L(2), L(3)), a2, (L(2), L(3), L(4)))
+        a_dest = contractalign([L(1), L(4)], a1, (L(1), L(2), L(3)), a2, (L(2), L(3), L(4)))
         @test a_dest ≈ a_ref
 
         # Empty labels (e.g. a scalar operand) are handled.
@@ -309,7 +316,7 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         @test eltype(a_dest) === elt_dest
         @test a_dest ≈ reshape(vec(a1) * transpose(vec(a2)), (size(a1)..., size(a2)...))
 
-        a_dest = contract(("i", "k", "j", "l"), a1, ("i", "j"), a2, ("k", "l"))
+        a_dest = contractalign(("i", "k", "j", "l"), a1, ("i", "j"), a2, ("k", "l"))
         @test eltype(a_dest) === elt_dest
         @test a_dest ≈ permutedims(
             reshape(vec(a1) * transpose(vec(a2)), (size(a1)..., size(a2)...)), (1, 3, 2, 4)
@@ -456,17 +463,17 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         size_dest_example = (3, 5, 2, 4)
 
         # Array-scalar contraction.
-        a_dest = contract(labels_dest_example, a, labels_a, s, ())
+        a_dest = contractalign(labels_dest_example, a, labels_a, s, ())
         @test size(a_dest) == size_dest_example
         @test a_dest ≈ permutedims(a, (2, 4, 1, 3)) * s[]
 
         # Scalar-array contraction.
-        a_dest = contract(labels_dest_example, s, (), a, labels_a)
+        a_dest = contractalign(labels_dest_example, s, (), a, labels_a)
         @test size(a_dest) == size_dest_example
         @test a_dest ≈ permutedims(a, (2, 4, 1, 3)) * s[]
 
         # Scalar-scalar contraction.
-        a_dest = contract((), s, (), t, ())
+        a_dest = contractalign((), s, (), t, ())
         @test size(a_dest) == ()
         @test a_dest[] ≈ s[] * t[]
 
diff --git a/test/test_exports.jl b/test/test_exports.jl
index 5cb1aea..89f67df 100644
--- a/test/test_exports.jl
+++ b/test/test_exports.jl
@@ -6,6 +6,7 @@ using Test: @test, @testset
         :TensorAlgebra,
         :contract,
         :contract!,
+        :contractalign,
         :dual,
         :eig_full,
         :eig_trunc,
@@ -37,12 +38,15 @@ using Test: @test, @testset
         append!(
             exports,
             [
-                :biperm, :bipartition, :cat_similar,
-                :concatenate, :concatenate!, :ContractAlgorithm, :contractopadd!, :data,
+                :allocate_output, :biperm, :bipartition, :cat_similar, :check_input,
+                :concatenate, :concatenate!, :ContractAlgorithm, :contractopadd!,
+                :contractperm, :contractperm!, :contractpermadd!, :contractpermopadd!,
+                :data,
                 :datatype, :directsum,
                 :flattenlinear, :is_output_view, :label_type,
                 :matricize, :matricizeop, :matricizeop!, :matricizeopcopy, :matricizeopview,
-                :permutedims, :permutedims!, :scalar, :similar_map,
+                :default_algorithm, :output_axes, :permutedims, :select_algorithm,
+                :permutedims!, :scalar, :similar_map,
                 :TensorOperationsAlgorithm,
                 :to_range, :tr, :tryflattenlinear, :ungrade, :zero!, :scale!,
                 :permuteddims, :PermutedDims,
diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl
index d9f62c6..c2daeaf 100644
--- a/test/test_factorizations.jl
+++ b/test/test_factorizations.jl
@@ -1,8 +1,8 @@
 using LinearAlgebra: LinearAlgebra, Diagonal, I, diag, norm
 using MatrixAlgebraKit: truncrank
-using TensorAlgebra: TensorAlgebra, contract, eig_full, eig_vals, eigh_full, eigh_vals,
-    left_null, left_orth, left_polar, lq_compact, lq_full, qr_compact, qr_full, right_null,
-    right_orth, right_polar, svd_compact, svd_full, svd_trunc, svd_vals
+using TensorAlgebra: TensorAlgebra, contract, contractalign, eig_full, eig_vals, eigh_full,
+    eigh_vals, left_null, left_orth, left_polar, lq_compact, lq_full, qr_compact, qr_full,
+    right_null, right_orth, right_polar, svd_compact, svd_full, svd_trunc, svd_vals
 using Test: @test, @testset
 using TestExtras: @constinferred
 
@@ -22,15 +22,15 @@ elts = (Float64, ComplexF64)
     Acopy = copy(A)
     Q, R = @constinferred qr_full(A, labels_A, labels_Q, labels_R)
     @test A == Acopy # should not have altered initial array
-    A′ = contract(labels_A, Q, (labels_Q..., :q), R, (:q, labels_R...))
+    A′ = contractalign(labels_A, Q, (labels_Q..., :q), R, (:q, labels_R...))
     @test A ≈ A′
     @test size(Q, 1) * size(Q, 2) == size(Q, 3) # Q is unitary
 
     Q, R = qr_full(A, (2, 1), (4, 3))
-    @test A ≈ contract(labels_A, Q, (labels_Q..., :q), R, (:q, labels_R...))
+    @test A ≈ contractalign(labels_A, Q, (labels_Q..., :q), R, (:q, labels_R...))
 
     Q, R = qr_full(A, Val(2))
-    @test A ≈ contract((:a, :b, :c, :d), Q, (:a, :b, :q), R, (:q, :c, :d))
+    @test A ≈ contractalign((:a, :b, :c, :d), Q, (:a, :b, :q), R, (:q, :c, :d))
 end
 
 @testset "Compact QR ($T)" for T in elts
@@ -42,7 +42,7 @@ end
     Acopy = copy(A)
     Q, R = @constinferred qr_compact(A, labels_A, labels_Q, labels_R)
     @test A == Acopy # should not have altered initial array
-    A′ = contract(labels_A, Q, (labels_Q..., :q), R, (:q, labels_R...))
+    A′ = contractalign(labels_A, Q, (labels_Q..., :q), R, (:q, labels_R...))
     @test A ≈ A′
     @test size(Q, 3) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4))
 end
@@ -58,12 +58,12 @@ end
     Acopy = copy(A)
     L, Q = @constinferred lq_full(A, labels_A, labels_L, labels_Q)
     @test A == Acopy # should not have altered initial array
-    A′ = contract(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...))
+    A′ = contractalign(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...))
     @test A ≈ A′
     @test size(Q, 1) == size(Q, 2) * size(Q, 3) # Q is unitary
 
     L, Q = lq_full(A, (2, 1), (4, 3))
-    @test A ≈ contract(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...))
+    @test A ≈ contractalign(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...))
 end
 
 @testset "Compact LQ ($T)" for T in elts
@@ -75,7 +75,7 @@ end
     Acopy = copy(A)
     L, Q = @constinferred lq_compact(A, labels_A, labels_L, labels_Q)
     @test A == Acopy # should not have altered initial array
-    A′ = contract(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...))
+    A′ = contractalign(labels_A, L, (labels_L..., :q), Q, (:q, labels_Q...))
     @test A ≈ A′
     @test size(Q, 1) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4)) # Q is unitary
 end
@@ -95,8 +95,8 @@ end
     # `D` is returned bare (the spectrum over the internal bond), which is a `Diagonal`.
     @test D isa Diagonal
 
-    AV = contract((:a, :b, :D), A, labels_A, V, (labels_V′..., :D))
-    VD = contract((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D))
+    AV = contractalign((:a, :b, :D), A, labels_A, V, (labels_V′..., :D))
+    VD = contractalign((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D))
     @test AV ≈ VD
 
     Dvals = eig_vals(A, labels_A, labels_V, labels_V′)
@@ -118,8 +118,8 @@ end
     @test eltype(V) == eltype(A)
     @test D isa Diagonal
 
-    AV = contract((:a, :b, :D), A, labels_A, V, (labels_V′..., :D))
-    VD = contract((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D))
+    AV = contractalign((:a, :b, :D), A, labels_A, V, (labels_V′..., :D))
+    VD = contractalign((:a, :b, :D), V, (labels_V..., :D′), D, (:D′, :D))
     @test AV ≈ VD
 
     Dvals = eigh_vals(A, labels_A, labels_V, labels_V′)
@@ -139,26 +139,26 @@ end
     U, S, Vᴴ = @constinferred svd_full(A, labels_A, labels_U, labels_Vᴴ)
     @test A == Acopy # should not have altered initial array
     US, labels_US = contract(U, (labels_U..., :u), S, (:u, :v))
-    A′ = contract(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...))
+    A′ = contractalign(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...))
     @test A ≈ A′
     @test size(U, 1) * size(U, 2) == size(U, 3) # U is unitary
     @test size(Vᴴ, 1) == size(Vᴴ, 2) * size(Vᴴ, 3) # V is unitary
 
     U, S, Vᴴ = svd_full(A, (2, 1), (4, 3))
     US, labels_US = contract(U, (labels_U..., :u), S, (:u, :v))
-    @test A ≈ contract(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...))
+    @test A ≈ contractalign(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...))
 
     U, S, Vᴴ = @constinferred svd_full(A, labels_A, labels_A, ())
     @test A == Acopy # should not have altered initial array
     US, labels_US = contract(U, (labels_A..., :u), S, (:u, :v))
-    A′ = contract(labels_A, US, labels_US, Vᴴ, (:v,))
+    A′ = contractalign(labels_A, US, labels_US, Vᴴ, (:v,))
     @test A ≈ A′
     @test size(Vᴴ, 1) == 1
 
     U, S, Vᴴ = @constinferred svd_full(A, labels_A, (), labels_A)
     @test A == Acopy # should not have altered initial array
     US, labels_US = contract(U, (:u,), S, (:u, :v))
-    A′ = contract(labels_A, US, labels_US, Vᴴ, (:v, labels_A...))
+    A′ = contractalign(labels_A, US, labels_US, Vᴴ, (:v, labels_A...))
     @test A ≈ A′
     @test size(U, 2) == 1
 end
@@ -174,7 +174,7 @@ end
     @test A == Acopy # should not have altered initial array
     @test S isa Diagonal
     US, labels_US = contract(U, (labels_U..., :u), S, (:u, :v))
-    A′ = contract(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...))
+    A′ = contractalign(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...))
     @test A ≈ A′
     k = min(size(S)...)
     @test size(U, 3) == k == size(Vᴴ, 1)
@@ -185,14 +185,14 @@ end
     U, S, Vᴴ = @constinferred svd_compact(A, labels_A, labels_A, ())
     @test A == Acopy # should not have altered initial array
     US, labels_US = contract(U, (labels_A..., :u), S, (:u, :v))
-    A′ = contract(labels_A, US, labels_US, Vᴴ, (:v,))
+    A′ = contractalign(labels_A, US, labels_US, Vᴴ, (:v,))
     @test A ≈ A′
     @test size(U, ndims(U)) == 1 == size(Vᴴ, 1)
 
     U, S, Vᴴ = @constinferred svd_compact(A, labels_A, (), labels_A)
     @test A == Acopy # should not have altered initial array
     US, labels_US = contract(U, (:u,), S, (:u, :v))
-    A′ = contract(labels_A, US, labels_US, Vᴴ, (:v, labels_A...))
+    A′ = contractalign(labels_A, US, labels_US, Vᴴ, (:v, labels_A...))
     @test A ≈ A′
     @test size(U, 1) == 1 == size(Vᴴ, 1)
 end
@@ -212,7 +212,7 @@ end
 
     @test A == Acopy # should not have altered initial array
     US, labels_US = contract(U, (labels_U..., :u), S, (:u, :v))
-    A′ = contract(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...))
+    A′ = contractalign(labels_A, US, labels_US, Vᴴ, (:v, labels_Vᴴ...))
     @test norm(A - A′) ≈ S_untrunc[end]
     @test size(S, 1) == size(S_untrunc, 1) - 1
     # `ϵ` is the 2-norm of the discarded singular values (here the single dropped value).
@@ -229,18 +229,36 @@ end
     N = @constinferred left_null(A, labels_A, labels_codomain, labels_domain)
     @test A == Acopy # should not have altered initial array
     # N^ba_n' * A^ba_dc = 0
-    NA = contract((:n, labels_domain...), conj(N), (labels_codomain..., :n), A, labels_A)
+    NA = contractalign(
+        (:n, labels_domain...),
+        conj(N),
+        (labels_codomain..., :n),
+        A,
+        labels_A
+    )
     @test norm(NA) ≈ 0 atol = 1.0e-14
     NN =
-        contract((:n, :n′), conj(N), (labels_codomain..., :n), N, (labels_codomain..., :n′))
+        contractalign(
+        (:n, :n′),
+        conj(N),
+        (labels_codomain..., :n),
+        N,
+        (labels_codomain..., :n′)
+    )
     @test NN ≈ LinearAlgebra.I
 
     Nᴴ = @constinferred right_null(A, labels_A, labels_codomain, labels_domain)
     @test A == Acopy # should not have altered initial array
     # A^ba_dc * N^dc_n' = 0
-    AN = contract((labels_codomain..., :n), A, labels_A, conj(Nᴴ), (:n, labels_domain...))
+    AN = contractalign(
+        (labels_codomain..., :n),
+        A,
+        labels_A,
+        conj(Nᴴ),
+        (:n, labels_domain...)
+    )
     @test norm(AN) ≈ 0 atol = 1.0e-14
-    NN = contract((:n, :n′), Nᴴ, (:n, labels_domain...), Nᴴ, (:n′, labels_domain...))
+    NN = contractalign((:n, :n′), Nᴴ, (:n, labels_domain...), Nᴴ, (:n′, labels_domain...))
 end
 
 @testset "Left polar ($T)" for T in elts
@@ -252,7 +270,7 @@ end
     Acopy = copy(A)
     W, P = left_polar(A, labels_A, labels_W, labels_P)
     @test A == Acopy # should not have altered initial array
-    A′ = contract(labels_A, W, (labels_W..., :w), P, (:w, labels_P...))
+    A′ = contractalign(labels_A, W, (labels_W..., :w), P, (:w, labels_P...))
     @test A ≈ A′
     @test size(W, 3) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4))
 end
@@ -266,7 +284,7 @@ end
     Acopy = copy(A)
     P, W = right_polar(A, labels_A, labels_P, labels_W)
     @test A == Acopy # should not have altered initial array
-    A′ = contract(labels_A, P, (labels_P..., :w), W, (:w, labels_W...))
+    A′ = contractalign(labels_A, P, (labels_P..., :w), W, (:w, labels_W...))
     @test A ≈ A′
     @test size(W, 1) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4))
 end
@@ -280,12 +298,12 @@ end
     Acopy = copy(A)
     W, P = left_orth(A, labels_A, labels_W, labels_P)
     @test A == Acopy # should not have altered initial array
-    A′ = contract(labels_A, W, (labels_W..., :w), P, (:w, labels_P...))
+    A′ = contractalign(labels_A, W, (labels_W..., :w), P, (:w, labels_P...))
     @test A ≈ A′
     @test size(W, 3) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4))
 
     W, P = left_orth(A, (2, 1), (4, 3))
-    @test A ≈ contract(labels_A, W, (labels_W..., :w), P, (:w, labels_P...))
+    @test A ≈ contractalign(labels_A, W, (labels_W..., :w), P, (:w, labels_P...))
 end
 
 @testset "Right orth ($T)" for T in elts
@@ -297,12 +315,12 @@ end
     Acopy = copy(A)
     P, W = right_orth(A, labels_A, labels_P, labels_W)
     @test A == Acopy # should not have altered initial array
-    A′ = contract(labels_A, P, (labels_P..., :w), W, (:w, labels_W...))
+    A′ = contractalign(labels_A, P, (labels_P..., :w), W, (:w, labels_W...))
     @test A ≈ A′
     @test size(W, 1) == min(size(A, 1) * size(A, 2), size(A, 3) * size(A, 4))
 
     P, W = right_orth(A, (2, 1), (4, 3))
-    @test A ≈ contract(labels_A, P, (labels_P..., :w), W, (:w, labels_W...))
+    @test A ≈ contractalign(labels_A, P, (labels_P..., :w), W, (:w, labels_W...))
 end
 
 # one (identity tensor)
diff --git a/test/test_matricize.jl b/test/test_matricize.jl
index 0582b4b..1b2dd58 100644
--- a/test/test_matricize.jl
+++ b/test/test_matricize.jl
@@ -108,7 +108,7 @@ end
     # for the identity destination bipermutation and the gather/scatter route otherwise.
     a1 = randn(rng, 2, 3, 5)
     a2 = randn(rng, 5, 3, 2)
-    ref = TensorAlgebra.contract((:i, :j, :k, :l), a1, (:i, :j, :m), a2, (:m, :k, :l))
+    ref = TensorAlgebra.contractalign((:i, :j, :k, :l), a1, (:i, :j, :m), a2, (:m, :k, :l))
     for labels in ((:i, :j, :k, :l), (:k, :l, :i, :j), (:k, :i, :l, :j))
         perm = map(l -> findfirst(==(l), (:i, :j, :k, :l)), labels)
         dest = randn(rng, map(d -> size(ref, d), perm)...)
diff --git a/test/test_mooncakeext.jl b/test/test_mooncakeext.jl
index 840b5f3..999f028 100644
--- a/test/test_mooncakeext.jl
+++ b/test/test_mooncakeext.jl
@@ -2,7 +2,7 @@ using Mooncake: Mooncake
 using Random: Random
 using TensorAlgebra: BiTuple, ContractAlgorithm, DefaultContractAlgorithm, Matricize,
     allocate_output, biperm, biperms, check_input, contract, contract!, contract_labels,
-    contractadd!, default_contract_algorithm, select_contract_algorithm
+    contractadd!, contractpermadd!, default_contract_algorithm, select_contract_algorithm
 using Test: @test, @testset
 
 @testset "MooncakeExt" begin
@@ -65,7 +65,7 @@ using Test: @test, @testset
     @testset "contract" begin
         α = true
         β = false
-        @testset "contractadd! (BiTuple)" begin
+        @testset "contractpermadd! (BiTuple)" begin
             dest = randn(elt, (2, 2))
             a1 = randn(elt, (2, 2))
             a2 = randn(elt, (2, 2))
@@ -73,7 +73,7 @@ using Test: @test, @testset
             biperm1 = BiTuple((1,), (2,))
             biperm2 = BiTuple((1,), (2,))
             Mooncake.TestUtils.test_rule(
-                rng, contractadd!, dest, biperm_dest.t1, biperm_dest.t2,
+                rng, contractpermadd!, dest, biperm_dest.t1, biperm_dest.t2,
                 a1, biperm1.t1, biperm1.t2, a2, biperm2.t1, biperm2.t2, α, β;
                 atol, rtol, mode, is_primitive
             )

From 678b3ba2e37e3de0b171185c1a27c2bab6e4c116 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 16 Sep 2026 01:35:06 -0400
Subject: [PATCH 16/17] Declare contractadd! public

It was the one rung of the labels ladder that was neither exported nor public, though it is as much a user-facing entry point as the rest.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 src/TensorAlgebra.jl | 2 +-
 test/test_exports.jl | 3 ++-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/TensorAlgebra.jl b/src/TensorAlgebra.jl
index b3b04fd..d1248a0 100644
--- a/src/TensorAlgebra.jl
+++ b/src/TensorAlgebra.jl
@@ -10,7 +10,7 @@ export contract, contract!, contractalign, dual, eig_full, eig_trunc, eig_vals,
 if VERSION >= v"1.11.0-DEV.469"
     eval(
         Meta.parse(
-            "public allocate_output, biperm, bipartition, cat_similar, check_input, concatenate, concatenate!, ContractAlgorithm, contractopadd!, contractperm, contractperm!, contractpermadd!, contractpermopadd!, data, datatype, directsum, flattenlinear, is_output_view, label_type, matricize, matricizeop, matricizeop!, matricizeopcopy, matricizeopview, output_axes, select_algorithm, default_algorithm, permutedims, permutedims!, scalar, similar_map, TensorOperationsAlgorithm, to_range, tr, tryflattenlinear, ungrade, zero!, scale!, permuteddims, PermutedDims"
+            "public allocate_output, biperm, bipartition, cat_similar, check_input, concatenate, concatenate!, ContractAlgorithm, contractadd!, contractopadd!, contractperm, contractperm!, contractpermadd!, contractpermopadd!, data, datatype, directsum, flattenlinear, is_output_view, label_type, matricize, matricizeop, matricizeop!, matricizeopcopy, matricizeopview, output_axes, select_algorithm, default_algorithm, permutedims, permutedims!, scalar, similar_map, TensorOperationsAlgorithm, to_range, tr, tryflattenlinear, ungrade, zero!, scale!, permuteddims, PermutedDims"
         )
     )
 end
diff --git a/test/test_exports.jl b/test/test_exports.jl
index 89f67df..40e6848 100644
--- a/test/test_exports.jl
+++ b/test/test_exports.jl
@@ -39,7 +39,8 @@ using Test: @test, @testset
             exports,
             [
                 :allocate_output, :biperm, :bipartition, :cat_similar, :check_input,
-                :concatenate, :concatenate!, :ContractAlgorithm, :contractopadd!,
+                :concatenate, :concatenate!, :ContractAlgorithm, :contractadd!,
+                :contractopadd!,
                 :contractperm, :contractperm!, :contractpermadd!, :contractpermopadd!,
                 :data,
                 :datatype, :directsum,

From 6100505d196d547d72def8fc57082319ab8cf8c5 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Wed, 16 Sep 2026 09:17:57 -0400
Subject: [PATCH 17/17] Drop the algorithm-positional non-mutating contractions

A contraction algorithm does not choose how the output is allocated, so these two entry points were a second way to contract that bypassed `allocate_output`. The algorithm stays a keyword above the in-place primitive.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 .../TensorAlgebraTensorOperationsExt.jl       | 35 -------------------
 1 file changed, 35 deletions(-)

diff --git a/ext/TensorAlgebraTensorOperationsExt/TensorAlgebraTensorOperationsExt.jl b/ext/TensorAlgebraTensorOperationsExt/TensorAlgebraTensorOperationsExt.jl
index 79c84f3..6017a1c 100644
--- a/ext/TensorAlgebraTensorOperationsExt/TensorAlgebraTensorOperationsExt.jl
+++ b/ext/TensorAlgebraTensorOperationsExt/TensorAlgebraTensorOperationsExt.jl
@@ -21,41 +21,6 @@ end
 # Using TensorOperations backends as TensorAlgebra implementations
 # ----------------------------------------------------------------
 
-# not in-place
-function TA.contractperm(
-        algorithm::TensorOperationsAlgorithm,
-        perm_dest_codomain, perm_dest_domain,
-        a1::AbstractArray, perm1_codomain, perm1_domain,
-        a2::AbstractArray, perm2_codomain, perm2_domain
-    )
-    permblocks1 = Tuple.((perm1_codomain, perm1_domain))
-    permblocks2 = Tuple.((perm2_codomain, perm2_domain))
-    permblocks_dest = Tuple.((perm_dest_codomain, perm_dest_domain))
-    conj1, conj2 = false, false
-    α = true
-    return TO.tensorcontract(
-        a1, permblocks1, conj1, a2, permblocks2, conj2,
-        permblocks_dest, α, backend(algorithm), allocator(algorithm)
-    )
-end
-
-function TA.contractalign(
-        algorithm::TensorOperationsAlgorithm,
-        labels_dest,
-        a1::AbstractArray, labels1,
-        a2::AbstractArray, labels2
-    )
-    permblocks1, permblocks2, permblocks_dest =
-        TO.contract_indices(labels1, labels2, labels_dest)
-    conj1, conj2 = false, false
-    α = true
-    return TO.tensorcontract(
-        a1, permblocks1, conj1, a2, permblocks2, conj2,
-        permblocks_dest, α, backend(algorithm), allocator(algorithm)
-    )
-end
-
-# in-place
 function TA.contractpermopadd!(
         algorithm::TensorOperationsAlgorithm,
         a_dest, perm_dest_codomain, perm_dest_domain,