diff --git a/Project.toml b/Project.toml
index 50f87a01..689ff5e5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "TensorAlgebra"
 uuid = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
-version = "0.20.0"
+version = "0.20.0-DEV"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
diff --git a/ext/TensorAlgebraTensorKitExt.jl b/ext/TensorAlgebraTensorKitExt.jl
index f386638a..40075a08 100644
--- a/ext/TensorAlgebraTensorKitExt.jl
+++ b/ext/TensorAlgebraTensorKitExt.jl
@@ -66,10 +66,10 @@ TensorAlgebra.sum(t::AbstractTensorMap; kwargs...) = Base.sum(convert(Array, t);
 # exactly what TensorKit's `similar(t, T, codomain, domain)` wants, so build the two
 # `ProductSpace`s directly.
 function TensorAlgebra.similar_map(
-        a::AbstractTensorMap, ::Type{T}, codomain_axes, domain_axes
+        a::AbstractTensorMap, ::Type{T}, axes_codomain, axes_domain
     ) where {T}
     S = spacetype(a)
-    return similar(a, T, ProductSpace{S}(codomain_axes...), ProductSpace{S}(domain_axes...))
+    return similar(a, T, ProductSpace{S}(axes_codomain...), ProductSpace{S}(axes_domain...))
 end
 
 # A plain-array prototype with native (space) axes is the operator/state construction case: `raw`
@@ -81,22 +81,22 @@ end
 # entries below can read it from whichever of codomain/domain is non-empty and share one builder,
 # mirroring `_map_homspace` and the map constructors.
 function similar_tensormap(
-        raw::AbstractArray, ::Type{T}, ::Type{S}, codomain_axes, domain_axes
+        raw::AbstractArray, ::Type{T}, ::Type{S}, axes_codomain, axes_domain
     ) where {T, S <: ElementarySpace}
     A = Base.promote_op(similar, typeof(raw), Type{T}, Int)
-    return TensorMapWithStorage{T, A}(undef, _map_homspace(S, codomain_axes, domain_axes))
+    return TensorMapWithStorage{T, A}(undef, _map_homspace(S, axes_codomain, axes_domain))
 end
 function TensorAlgebra.similar_map(
         raw::AbstractArray, ::Type{T},
-        codomain_axes::Tuple{S, Vararg{S}}, domain_axes::Tuple{Vararg{S}}
+        axes_codomain::Tuple{S, Vararg{S}}, axes_domain::Tuple{Vararg{S}}
     ) where {T, S <: ElementarySpace}
-    return similar_tensormap(raw, T, S, codomain_axes, domain_axes)
+    return similar_tensormap(raw, T, S, axes_codomain, axes_domain)
 end
 function TensorAlgebra.similar_map(
         raw::AbstractArray, ::Type{T},
-        codomain_axes::Tuple{}, domain_axes::Tuple{S, Vararg{S}}
+        axes_codomain::Tuple{}, axes_domain::Tuple{S, Vararg{S}}
     ) where {T, S <: ElementarySpace}
-    return similar_tensormap(raw, T, S, codomain_axes, domain_axes)
+    return similar_tensormap(raw, T, S, axes_codomain, axes_domain)
 end
 
 # ===============================  zeros_map / randn_map / rand_map  ========================
@@ -107,32 +107,32 @@ end
 # `S` is passed to `_map_homspace` explicitly so the two dispatch entries per constructor can
 # read it from whichever of the codomain/domain is non-empty and share one builder; an empty
 # axis tuple gives the unit space `ProductSpace{S}()`.
-function _map_homspace(::Type{S}, codomain_axes, domain_axes) where {S <: ElementarySpace}
-    return ProductSpace{S}(codomain_axes...) ← ProductSpace{S}(domain_axes...)
+function _map_homspace(::Type{S}, axes_codomain, axes_domain) where {S <: ElementarySpace}
+    return ProductSpace{S}(axes_codomain...) ← ProductSpace{S}(axes_domain...)
 end
 function TensorAlgebra.zeros_map(
-        ::Type{T}, codomain_axes::Tuple{S, Vararg{S}}, domain_axes::Tuple{Vararg{S}}
+        ::Type{T}, axes_codomain::Tuple{S, Vararg{S}}, axes_domain::Tuple{Vararg{S}}
     ) where {T, S <: ElementarySpace}
-    return TensorKit.zeros(T, _map_homspace(S, codomain_axes, domain_axes))
+    return TensorKit.zeros(T, _map_homspace(S, axes_codomain, axes_domain))
 end
 function TensorAlgebra.zeros_map(
-        ::Type{T}, codomain_axes::Tuple{}, domain_axes::Tuple{S, Vararg{S}}
+        ::Type{T}, axes_codomain::Tuple{}, axes_domain::Tuple{S, Vararg{S}}
     ) where {T, S <: ElementarySpace}
-    return TensorKit.zeros(T, _map_homspace(S, codomain_axes, domain_axes))
+    return TensorKit.zeros(T, _map_homspace(S, axes_codomain, axes_domain))
 end
 for (f, g) in ((:randn_map, :randn), (:rand_map, :rand))
     @eval begin
         function TensorAlgebra.$f(
                 rng::AbstractRNG, ::Type{T},
-                codomain_axes::Tuple{S, Vararg{S}}, domain_axes::Tuple{Vararg{S}}
+                axes_codomain::Tuple{S, Vararg{S}}, axes_domain::Tuple{Vararg{S}}
             ) where {T, S <: ElementarySpace}
-            return TensorKit.$g(rng, T, _map_homspace(S, codomain_axes, domain_axes))
+            return TensorKit.$g(rng, T, _map_homspace(S, axes_codomain, axes_domain))
         end
         function TensorAlgebra.$f(
                 rng::AbstractRNG, ::Type{T},
-                codomain_axes::Tuple{}, domain_axes::Tuple{S, Vararg{S}}
+                axes_codomain::Tuple{}, axes_domain::Tuple{S, Vararg{S}}
             ) where {T, S <: ElementarySpace}
-            return TensorKit.$g(rng, T, _map_homspace(S, codomain_axes, domain_axes))
+            return TensorKit.$g(rng, T, _map_homspace(S, axes_codomain, axes_domain))
         end
     end
 end
@@ -142,7 +142,7 @@ end
 # `TensorMap` is not an `AbstractArray`, so the generic `copyto!` default does not apply; delegate
 # to TensorKit's `project_symmetric!`, which fills the symmetry-allowed blocks from the dense data
 # and discards any component outside the block structure. Composed with the map constructors above,
-# this makes `project(dense, codomain_axes, domain_axes)` build a `TensorMap` from a dense matrix.
+# this makes `project(dense, axes_codomain, axes_domain)` build a `TensorMap` from a dense matrix.
 # `project_symmetric!` requires a matching dense size, so reshape `src` to `size(dest)` first (a
 # no-op when the ranks already match); this lets a lower-rank `src` omit trailing length-1 axes,
 # matching the generic `projectto!`, and rejects a genuine shape mismatch.
@@ -168,14 +168,14 @@ end
 # whichever side is non-empty, the same two-entry split `similar_map` uses. `project`'s allocation
 # stays generic (strict `similar_map`); only the `project_aux` derivation is `TensorMap`-specific.
 function TensorAlgebra.infer_aux_space(
-        raw::AbstractArray, codomain_axes::Tuple{S, Vararg{S}}, domain_axes::Tuple{Vararg{S}}
+        raw::AbstractArray, axes_codomain::Tuple{S, Vararg{S}}, axes_domain::Tuple{Vararg{S}}
     ) where {S <: ElementarySpace}
-    return infer_aux_space_tensormap(raw, S, codomain_axes, domain_axes)
+    return infer_aux_space_tensormap(raw, S, axes_codomain, axes_domain)
 end
 function TensorAlgebra.infer_aux_space(
-        raw::AbstractArray, codomain_axes::Tuple{}, domain_axes::Tuple{S, Vararg{S}}
+        raw::AbstractArray, axes_codomain::Tuple{}, axes_domain::Tuple{S, Vararg{S}}
     ) where {S <: ElementarySpace}
-    return infer_aux_space_tensormap(raw, S, codomain_axes, domain_axes)
+    return infer_aux_space_tensormap(raw, S, axes_codomain, axes_domain)
 end
 
 # The space of `raw`'s trailing auxiliary axis, derived so the projected result is
@@ -184,11 +184,11 @@ end
 # follows, so the aux slices must appear in that order. The result may span several sectors (a
 # direct-sum, MPO-style virtual leg).
 function infer_aux_space_tensormap(
-        raw, ::Type{S}, codomain_axes, domain_axes
+        raw, ::Type{S}, axes_codomain, axes_domain
     ) where {S <: ElementarySpace}
-    aux_dim = length(codomain_axes) + length(domain_axes) + 1
+    aux_dim = length(axes_codomain) + length(axes_domain) + 1
     aux_length = size(raw, aux_dim)
-    content = fuse(codomain_axes..., dual.(domain_axes)...)
+    content = fuse(axes_codomain..., dual.(axes_domain)...)
     # Probe the surplus axis slice by slice: a slice keeps the aux axis (width `dim(s)`), so its
     # rank matches the candidate axes exactly and `tryproject` allocates, fills, and round-trip-
     # verifies without re-entering the derivation branch. This builds one `TensorMap` per candidate
@@ -197,7 +197,7 @@ function infer_aux_space_tensormap(
     function slice_is_covariant(r, s)
         slice = selectdim(raw, aux_dim, r)
         return !isnothing(
-            TensorAlgebra.tryproject(slice, codomain_axes, (domain_axes..., S(s => 1)))
+            TensorAlgebra.tryproject(slice, axes_codomain, (axes_domain..., S(s => 1)))
         )
     end
     seccounts = Pair{TensorKit.sectortype(S), Int}[]
@@ -264,10 +264,10 @@ end
 # codomain-facing (un-dualized), which is exactly TensorKit's domain convention, so they build the
 # domain `ProductSpace` directly.
 function TensorAlgebra.unmatricize(
-        ::TensorKitMatricize, m::AbstractTensorMap, codomain_axes, domain_axes
+        ::TensorKitMatricize, m::AbstractTensorMap, axes_codomain, axes_domain
     )
     S = spacetype(m)
-    dest = ProductSpace{S}(codomain_axes...) ← ProductSpace{S}(domain_axes...)
+    dest = ProductSpace{S}(axes_codomain...) ← ProductSpace{S}(axes_domain...)
     space(m) == dest ||
         throw(ArgumentError("`unmatricize` space `$dest` does not match `$(space(m))`"))
     return m
diff --git a/src/TensorAlgebra.jl b/src/TensorAlgebra.jl
index df3f74a4..2853d1f3 100644
--- a/src/TensorAlgebra.jl
+++ b/src/TensorAlgebra.jl
@@ -23,7 +23,6 @@ include("permutedimsadd.jl")
 include("matricize.jl")
 include("concatenate.jl")
 include("directsum.jl")
-include("diagonal.jl")
 include("dual.jl")
 include("to_range.jl")
 include("contract/contractalgorithm.jl")
@@ -31,6 +30,7 @@ include("contract/contract.jl")
 include("contract/contract_labels.jl")
 include("contract/biperms.jl")
 include("contract/allocate_output.jl")
+include("diagonal.jl")
 include("contract/contract_matricize.jl")
 include("factorizations.jl")
 include("matrixfunctions.jl")
diff --git a/src/contract/allocate_output.jl b/src/contract/allocate_output.jl
index 994e21a9..c09f6716 100644
--- a/src/contract/allocate_output.jl
+++ b/src/contract/allocate_output.jl
@@ -64,7 +64,13 @@ function output_axes(
     axes_codomain, _ = bipartition(axes(a1), perm1_codomain, perm1_domain)
     _, axes_domain = bipartition(axes(a2), perm2_codomain, perm2_domain)
     axes_uncontracted = (axes_codomain..., axes_domain...)
-    return bipartition(axes_uncontracted, perm_dest_codomain, perm_dest_domain)
+    axes_codomain_dest, axes_domain_dest = bipartition(
+        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).
+    return axes_codomain_dest, conj.(axes_domain_dest)
 end
 
 # TODO: Use `ArrayLayouts`-like `MulAdd` object,
@@ -84,14 +90,21 @@ function allocate_output(
         perm2_codomain,
         perm2_domain
     )
-    codomain_axes_dest, domain_axes_dest = output_axes(
+    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
     )
-    T = promote_type(eltype(a1), eltype(a2))
-    # `domain_axes_dest` come straight from `axes(a2)` (stored/dualized convention), so
-    # un-dualize them into `similar_map`'s codomain-facing convention.
-    return zero!(similar_map(a1, T, codomain_axes_dest, conj.(domain_axes_dest)))
+    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))
 end
diff --git a/src/diagonal.jl b/src/diagonal.jl
index 23296a49..bbb62cd3 100644
--- a/src/diagonal.jl
+++ b/src/diagonal.jl
@@ -5,10 +5,9 @@ using LinearAlgebra: Diagonal
 # an operation is still diagonal. These methods hook the lowest-level primitives, so the
 # convenience wrappers built on them (`bipermutedims`, `permutedimsadd!`, `add!`, and the
 # matrix functions, which all route through `bipermutedimsopadd!` and `allocate_output`)
-# preserve `Diagonal`. Structure is given up (via the generic reshape path) only where the
-# result genuinely is not diagonal: vectorizing matricizations and `Diagonal`/dense mixing
-# (the latter falls back to Base's dense `similar`, since `contract` allocates from a flat
-# axis tuple, not a `BiTuple`).
+# preserve `Diagonal`. Structure is given up only where the result genuinely is not diagonal:
+# vectorizing matricizations, a non-`{1,1}` bond-split `unmatricize`, `Diagonal`/dense mixing,
+# and `contract` patterns other than the single-contracted-leg matmul (`Diagonal * Diagonal`).
 
 # Permuting the two axes of a square `Diagonal` (identity or transpose) leaves it
 # unchanged, so the lazy permutation is the matrix itself.
@@ -45,9 +44,33 @@ end
 # A `Diagonal` is already a matrix; the `(1 codomain, 1 domain)` matricization is the identity
 # reshape, so return it directly (maybe-alias, matching `matricize`'s general contract).
 matricize(::ReshapeMatricize, a::Diagonal, ::Val{1}) = a
+# 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.
 function unmatricize(
         ::ReshapeMatricize, m::Diagonal,
-        ::Tuple{<:AbstractUnitRange}, ::Tuple{<:AbstractUnitRange}
+        axes_codomain::Tuple{<:AbstractUnitRange}, axes_domain::Tuple{<:AbstractUnitRange}
     )
+    check_input(unmatricize, m, axes_codomain, axes_domain)
     return m
 end
+# Any other split is a genuine bond-split (for example `unmatricize(D[4×4], (2, 2), (4,))`) whose
+# result is not representable as a `Diagonal`, so densify and reshape like a dense matrix.
+# `copyto!(similar(m, axes(m)), m)` densifies while preserving `m`'s array backend (a plain
+# `Array` would force the result onto the CPU).
+function unmatricize(
+        style::ReshapeMatricize, m::Diagonal, axes_codomain::Tuple, axes_domain::Tuple
+    )
+    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}
+    )
+    return Diagonal(zero!(similar(a1.diag, T, (only(axes_codomain),))))
+end
diff --git a/src/factorizations.jl b/src/factorizations.jl
index a4be3ed1..11a38a4d 100644
--- a/src/factorizations.jl
+++ b/src/factorizations.jl
@@ -1,8 +1,9 @@
 using LinearAlgebra: LinearAlgebra
 using MatrixAlgebraKit: MatrixAlgebraKit
 
-# Each factorization reconstructs its factors with `unmatricize`, reading the freshly created
-# bond axis off the factor itself: it is the factor's last axis on a codomain factor
+# Each factorization reconstructs its outer factors with `unmatricize`, while a spectrum factor
+# (`S`, `D`) is returned bare at the matrix level. The freshly created bond axis is read off the
+# factor itself: it is the factor's last axis on a codomain factor
 # (`[group…, bond]`) and its first axis on a domain factor (`[bond, group…]`), on every backend
 # (a fusing backend returns a rank-2 factor, a `TensorMap` keeps the group's original legs). The
 # bond is dualized to codomain-facing form (`conj`, a no-op on a dense axis) when it lands on the
@@ -261,7 +262,7 @@ for f in (:svd_compact, :svd_full)
             U, S, Vᴴ = MatrixAlgebraKit.$f(A_mat; kwargs...)
             axes_codomain, axes_domain = bipartition_axes(axes(A), ndims_codomain)
             return unmatricize(style, U, axes_codomain, (conj(axes(U, ndims(U))),)),
-                unmatricize(style, S, (axes(S, 1),), (conj(axes(S, 2)),)),
+                S,
                 unmatricize(style, Vᴴ, (axes(Vᴴ, 1),), axes_domain)
         end
         function $f(A, ndims_codomain::Val; kwargs...)
@@ -278,7 +279,7 @@ function svd_trunc(style::MatricizeStyle, A, ndims_codomain::Val; kwargs...)
     U, S, Vᴴ, ϵ = MatrixAlgebraKit.svd_trunc(A_mat; kwargs...)
     axes_codomain, axes_domain = bipartition_axes(axes(A), ndims_codomain)
     return unmatricize(style, U, axes_codomain, (conj(axes(U, ndims(U))),)),
-        unmatricize(style, S, (axes(S, 1),), (conj(axes(S, 2)),)),
+        S,
         unmatricize(style, Vᴴ, (axes(Vᴴ, 1),), axes_domain),
         ϵ
 end
@@ -287,14 +288,15 @@ function svd_trunc(A, ndims_codomain::Val; kwargs...)
 end
 
 # Eigendecomposition: `D` is the rank × rank spectrum and `V` carries the codomain axes plus a
-# trailing rank axis. Both are unmatricized back to the array type, as in `svd_*`.
+# trailing rank axis. `D` is returned bare (its axis is the internal bond, so there is nothing to
+# unfold); `V` is unmatricized back to the array type, as in `svd_*`.
 for f in (:eigh_full, :eig_full, :eigh_trunc, :eig_trunc)
     @eval begin
         function $f(style::MatricizeStyle, A, ndims_codomain::Val; kwargs...)
             A_mat = matricize(style, A, ndims_codomain)
             D, V = MatrixAlgebraKit.$f(A_mat; kwargs...)
             axes_codomain = first(bipartition(axes(A), ndims_codomain))
-            return unmatricize(style, D, (axes(D, 1),), (conj(axes(D, 2)),)),
+            return D,
                 unmatricize(style, V, axes_codomain, (conj(axes(V, ndims(V))),))
         end
         function $f(A, ndims_codomain::Val; kwargs...)
@@ -812,8 +814,8 @@ one
 function one!!(style::MatricizeStyle, A, ndims_codomain::Val; kwargs...)
     A_mat = matricize(style, A, ndims_codomain)
     MatrixAlgebraKit.one!(A_mat)
-    codomain_axes, domain_axes = bipartition_axes(axes(A), ndims_codomain)
-    return unmatricize(style, A_mat, codomain_axes, domain_axes)
+    axes_codomain, axes_domain = bipartition_axes(axes(A), ndims_codomain)
+    return unmatricize(style, A_mat, axes_codomain, axes_domain)
 end
 function one!!(A, ndims_codomain::Val; kwargs...)
     return one!!(MatricizeStyle(A), A, ndims_codomain; kwargs...)
diff --git a/src/linearbroadcasted.jl b/src/linearbroadcasted.jl
index 9d4bdeac..1c8e9174 100644
--- a/src/linearbroadcasted.jl
+++ b/src/linearbroadcasted.jl
@@ -145,7 +145,9 @@ arguments(a::AddBroadcasted) = addends(a)
 # Mul — lazy matrix multiplication (standalone, not LinearBroadcasted)
 # ---------------------------------------------------------------------------- #
 
-# Same as `LinearAlgebra.matprod`, but duplicated here since it is private.
+# The element type of a matrix-product accumulation, modeled on `LinearAlgebra.matprod` (private
+# there, so inlined): a sum of products, so `promote_op(matprod, ...)` widens exactly as an
+# accumulating inner loop would.
 matprod(x, y) = x * y + x * y
 
 struct Mul{A, B}
diff --git a/src/matricize.jl b/src/matricize.jl
index 5b5a9bde..ed6cd1c0 100644
--- a/src/matricize.jl
+++ b/src/matricize.jl
@@ -209,18 +209,18 @@ function matricizeopperm(
 end
 
 # ====================================  unmatricize  =======================================
-# Split form: `codomain_axes` and `domain_axes` are the destination axes for the codomain and
+# 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`
 # (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
 # disambiguate axis tuples from permutation tuples regardless of how unconstrained `m` and the
 # axes are.
-function unmatricize(style::MatricizeStyle, m, codomain_axes, domain_axes)
-    return throw(MethodError(unmatricize, (style, m, codomain_axes, domain_axes)))
+function unmatricize(style::MatricizeStyle, m, axes_codomain, axes_domain)
+    return throw(MethodError(unmatricize, (style, m, axes_codomain, axes_domain)))
 end
-function unmatricize(m, codomain_axes, domain_axes)
-    return unmatricize(MatricizeStyle(m), m, codomain_axes, domain_axes)
+function unmatricize(m, axes_codomain, axes_domain)
+    return unmatricize(MatricizeStyle(m), m, axes_codomain, axes_domain)
 end
 
 # Split `axes` into its codomain and domain groups like `bipartition`, but present the domain
@@ -228,8 +228,8 @@ end
 # take. The domain axes `bipartition` reads off `axes(a)` are in the stored (dualized) form, so
 # this bridges from `axes(a)` to the `unmatricize` axis convention (a no-op on dense axes).
 function bipartition_axes(t::Tuple, split...)
-    codomain_axes, domain_axes = bipartition(t, split...)
-    return codomain_axes, conj.(domain_axes)
+    axes_codomain, axes_domain = bipartition(t, split...)
+    return axes_codomain, conj.(axes_domain)
 end
 
 # Inverse-bipermutation form: split `axes_dest` into codomain/domain groups reordered by the
@@ -253,8 +253,8 @@ function unmatricizeperm(
     invbiperm = BiTuple(invperm_codomain, invperm_domain)
     length(axes_dest) == length(invbiperm) ||
         throw(ArgumentError("axes do not match permutation"))
-    codomain_axes, domain_axes = bipartition_axes(axes_dest, invbiperm)
-    a12 = unmatricize(style, m, codomain_axes, domain_axes)
+    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
@@ -272,8 +272,8 @@ function unmatricizeperm!(
     invbiperm = BiTuple(invperm_codomain, invperm_domain)
     ndims(a_dest) == length(invbiperm) ||
         throw(ArgumentError("destination does not match permutation"))
-    codomain_axes, domain_axes = bipartition_axes(axes(a_dest), invbiperm)
-    a_perm = unmatricize(style, m, codomain_axes, domain_axes)
+    axes_codomain, axes_domain = bipartition_axes(axes(a_dest), invbiperm)
+    a_perm = unmatricize(style, m, axes_codomain, axes_domain)
     biperm_dest = BiTuple(Tuple(invperm(invbiperm)), Val(length_codomain(invbiperm)))
     return bipermutedims!(a_dest, a_perm, biperm_dest)
 end
@@ -317,8 +317,20 @@ end
 function matricizepermaliases(style::ReshapeMatricize, perm_codomain, perm_domain)
     return matricizekind(style, perm_codomain, perm_domain) != PermuteMatricizeKind
 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
+# would reshape silently.
+function check_input(::typeof(unmatricize), m, axes_codomain, axes_domain)
+    (
+        ndims(m) == 2 &&
+            size(m, 1) == prod(length, axes_codomain; init = 1) &&
+            size(m, 2) == prod(length, axes_domain; init = 1)
+    ) || throw(DimensionMismatch("`unmatricize` axes do not match the matrix size"))
+    return nothing
+end
 # A dense reshape ignores the codomain/domain split: it just reshapes to the concatenated axes.
-# `conj` re-dualizes the codomain-facing `domain_axes` into stored form, a no-op on a dense axis.
-function unmatricize(style::ReshapeMatricize, m, codomain_axes, domain_axes)
-    return reshape(m, (codomain_axes..., conj.(domain_axes)...))
+# `conj` re-dualizes the codomain-facing `axes_domain` into stored form, a no-op on a dense axis.
+function unmatricize(style::ReshapeMatricize, m, axes_codomain, axes_domain)
+    check_input(unmatricize, m, axes_codomain, axes_domain)
+    return reshape(m, (axes_codomain..., conj.(axes_domain)...))
 end
diff --git a/src/matrixfunctions.jl b/src/matrixfunctions.jl
index 85240c08..5ccdc460 100644
--- a/src/matrixfunctions.jl
+++ b/src/matrixfunctions.jl
@@ -36,8 +36,8 @@ for f in MATRIX_FUNCTIONS
         function $f(style::MatricizeStyle, a, ndims_codomain::Val; kwargs...)
             a_mat = matricize(style, a, ndims_codomain)
             fa_mat = Base.$f(a_mat; kwargs...)
-            codomain_axes, domain_axes = bipartition_axes(axes(a), ndims_codomain)
-            return unmatricize(style, fa_mat, codomain_axes, domain_axes)
+            axes_codomain, axes_domain = bipartition_axes(axes(a), ndims_codomain)
+            return unmatricize(style, fa_mat, axes_codomain, axes_domain)
         end
         function $f(a, ndims_codomain::Val; kwargs...)
             return $f(MatricizeStyle(a), a, ndims_codomain; kwargs...)
diff --git a/src/permutedimsadd.jl b/src/permutedimsadd.jl
index dd85d666..c94510d6 100644
--- a/src/permutedimsadd.jl
+++ b/src/permutedimsadd.jl
@@ -227,10 +227,10 @@ function permutedims(a, perm)
     return permutedims!(dest, a, perm)
 end
 function permutedims(a, perm_codomain, perm_domain)
-    codomain_axes = map(p -> axes(a, p), perm_codomain)
+    axes_codomain = map(p -> axes(a, p), perm_codomain)
     # `similar_map` dualizes the domain axes it is given (as it does for `similar_map`
     # itself), so pass them pre-conjugated to land back on the source's own axes.
-    domain_axes = map(p -> conj(axes(a, p)), perm_domain)
-    dest = similar_map(a, eltype(a), codomain_axes, domain_axes)
+    axes_domain = map(p -> conj(axes(a, p)), perm_domain)
+    dest = similar_map(a, eltype(a), axes_codomain, axes_domain)
     return permutedims!(dest, a, perm_codomain, perm_domain)
 end
diff --git a/src/projectto.jl b/src/projectto.jl
index 25c445ae..60a3418a 100644
--- a/src/projectto.jl
+++ b/src/projectto.jl
@@ -28,20 +28,20 @@ function projectto!(dest, src)
 end
 
 """
-    allocate_project(raw, codomain_axes, domain_axes) -> dest
+    allocate_project(raw, axes_codomain, axes_domain) -> dest
 
 Allocate the destination that projecting `raw` onto
-`codomain_axes`/`domain_axes` fills. This is a backend customization point
+`axes_codomain`/`axes_domain` fills. This is a backend customization point
 (with [`projectto!`](@ref) and [`is_projected`](@ref)); the default is plain
-`similar_map(raw, codomain_axes, domain_axes)`.
+`similar_map(raw, axes_codomain, axes_domain)`.
 
 `project` projects into exactly the given axes, so `raw` must not have more
 axes than they account for. To append a derived flux-carrying auxiliary axis
 for a charge-shifting operator or a non-invariant state, use
 [`project_aux`](@ref) instead.
 """
-function allocate_project(raw, codomain_axes, domain_axes)
-    nphys = length(codomain_axes) + length(domain_axes)
+function allocate_project(raw, axes_codomain, axes_domain)
+    nphys = length(axes_codomain) + length(axes_domain)
     ndims(raw) <= nphys || throw(
         ArgumentError(
             "`project` projects into exactly the given axes and does not derive an auxiliary \
@@ -49,25 +49,25 @@ function allocate_project(raw, codomain_axes, domain_axes)
             append a derived flux-carrying leg, or pass the axis explicitly."
         )
     )
-    return similar_map(raw, codomain_axes, domain_axes)
+    return similar_map(raw, axes_codomain, axes_domain)
 end
 
 """
-    unchecked_project(raw, codomain_axes, domain_axes) -> dest
+    unchecked_project(raw, axes_codomain, axes_domain) -> dest
     unchecked_project(raw, axes) -> dest
 
 Project `raw` into a symmetry-restricted array shaped as a map from
-`domain_axes` to `codomain_axes`, without checking which components are
+`axes_domain` to `axes_codomain`, without checking which components are
 discarded: entries of `raw` outside the symmetry-allowed structure are
 dropped without inspection. Most callers want [`project`](@ref), which
 verifies that nothing was discarded, or [`tryproject`](@ref), its nullable
 sibling. All three derive from the backend customization points: this one is
-`projectto!(allocate_project(raw, codomain_axes, domain_axes), raw)`. The
+`projectto!(allocate_project(raw, axes_codomain, axes_domain), raw)`. The
 two-argument form takes a flat list of `axes` and is equivalent to an empty
 domain.
 """
-function unchecked_project(raw, codomain_axes, domain_axes)
-    return projectto!(allocate_project(raw, codomain_axes, domain_axes), raw)
+function unchecked_project(raw, axes_codomain, axes_domain)
+    return projectto!(allocate_project(raw, axes_codomain, axes_domain), raw)
 end
 # The flat all-codomain (state) form: a list of `axes` with an empty domain.
 unchecked_project(raw, axes) = unchecked_project(raw, axes, ())
@@ -103,7 +103,7 @@ Inverse of [`project`](@ref): recover the dense array that `project` maps to `a`
 codomain/domain split `ndims_codomain` as a `Val`. The default is `convert(Array, a)`; a backend
 that changes basis in `project` overloads this to undo that change, so that
 
-    unproject(project(raw, codomain_axes, domain_axes), Val(length(codomain_axes))) ≈ raw
+    unproject(project(raw, axes_codomain, axes_domain), Val(length(axes_codomain))) ≈ raw
 """
 unproject(a, ::Val) = convert(Array, a)
 
@@ -126,33 +126,33 @@ function project!(dest, src; kwargs...)
 end
 
 """
-    project(raw, codomain_axes, domain_axes; kwargs...) -> dest
+    project(raw, axes_codomain, axes_domain; kwargs...) -> dest
     project(raw, axes; kwargs...) -> dest
 
 Project `raw` into a symmetry-restricted array shaped as a map from
-`domain_axes` to `codomain_axes`, verifying that only a negligible component
+`axes_domain` to `axes_codomain`, verifying that only a negligible component
 of `raw` is discarded and throwing an `InexactError` otherwise (keyword
 arguments are forwarded to the `isapprox` tolerance check; the default
 tolerances are subject to change in future versions). See
 [`tryproject`](@ref) for a nullable version and [`unchecked_project`](@ref)
 for the unchecked projection this derives from.
 
-`raw` must not have more axes than `codomain_axes`/`domain_axes` account for:
+`raw` must not have more axes than `axes_codomain`/`axes_domain` account for:
 `project` projects into exactly the given axes. To append a derived
 flux-carrying auxiliary axis (for a charge-shifting operator or a
 non-invariant state), use [`project_aux`](@ref). The two-argument form takes a
 flat list of `axes` and is equivalent to an empty domain.
 """
-function project(raw, codomain_axes, domain_axes; kwargs...)
-    dest = unchecked_project(raw, codomain_axes, domain_axes)
-    is_projected(dest, raw, Val(length(codomain_axes)); kwargs...) ||
+function project(raw, axes_codomain, axes_domain; kwargs...)
+    dest = unchecked_project(raw, axes_codomain, axes_domain)
+    is_projected(dest, raw, Val(length(axes_codomain)); kwargs...) ||
         throw(InexactError(:project, typeof(dest), raw))
     return dest
 end
 project(raw, axes; kwargs...) = project(raw, axes, (); kwargs...)
 
 """
-    tryproject(raw, codomain_axes, domain_axes; kwargs...) -> Union{dest, Nothing}
+    tryproject(raw, axes_codomain, axes_domain; kwargs...) -> Union{dest, Nothing}
     tryproject(raw, axes; kwargs...) -> Union{dest, Nothing}
 
 Like [`project`](@ref), but return `nothing` instead of throwing when more
@@ -165,14 +165,14 @@ flux-carrying leg:
 
 Keyword arguments are forwarded to the `isapprox` tolerance check.
 """
-function tryproject(raw, codomain_axes, domain_axes; kwargs...)
-    dest = unchecked_project(raw, codomain_axes, domain_axes)
-    return is_projected(dest, raw, Val(length(codomain_axes)); kwargs...) ? dest : nothing
+function tryproject(raw, axes_codomain, axes_domain; kwargs...)
+    dest = unchecked_project(raw, axes_codomain, axes_domain)
+    return is_projected(dest, raw, Val(length(axes_codomain)); kwargs...) ? dest : nothing
 end
 tryproject(raw, axes; kwargs...) = tryproject(raw, axes, (); kwargs...)
 
 """
-    infer_aux_space(raw, codomain_axes, domain_axes) -> aux
+    infer_aux_space(raw, axes_codomain, axes_domain) -> aux
 
 Derive the auxiliary axis the `*_aux` projection verbs append as the last
 domain axis, so the projected result is symmetry-allowed. `raw` carries the
@@ -182,16 +182,16 @@ from `raw`, while a symmetric backend reads it from the sector structure (a
 graded backend derives per-slice sectors, the `TensorMap` backend scans the
 `codomain ⊗ conj(domain)` content).
 """
-function infer_aux_space(raw, codomain_axes, domain_axes)
-    return axes(raw, length(codomain_axes) + length(domain_axes) + 1)
+function infer_aux_space(raw, axes_codomain, axes_domain)
+    return axes(raw, length(axes_codomain) + length(axes_domain) + 1)
 end
 
 # Reshape a physical-rank `raw` up to one trailing slice axis, derive the auxiliary space, and
-# return the `(raw, codomain_axes, domain_axes)` triple to forward to a projection verb, with the
+# return the `(raw, axes_codomain, axes_domain)` triple to forward to a projection verb, with the
 # aux appended to the domain. A rank beyond one surplus axis is an error. Shared by the three
 # `*_aux` verbs below.
-function project_aux_args(raw, codomain_axes, domain_axes)
-    nphys = length(codomain_axes) + length(domain_axes)
+function project_aux_args(raw, axes_codomain, axes_domain)
+    nphys = length(axes_codomain) + length(axes_domain)
     nphys <= ndims(raw) <= nphys + 1 || throw(
         ArgumentError(
             "`project_aux` expected a rank-$nphys or rank-$(nphys + 1) input for $nphys given \
@@ -199,12 +199,12 @@ function project_aux_args(raw, codomain_axes, domain_axes)
         )
     )
     slices = ndims(raw) == nphys ? reshape(raw, (size(raw)..., 1)) : raw
-    aux = infer_aux_space(slices, codomain_axes, domain_axes)
-    return slices, codomain_axes, (domain_axes..., aux)
+    aux = infer_aux_space(slices, axes_codomain, axes_domain)
+    return slices, axes_codomain, (axes_domain..., aux)
 end
 
 """
-    project_aux(raw, codomain_axes, domain_axes; kwargs...) -> dest
+    project_aux(raw, axes_codomain, axes_domain; kwargs...) -> dest
     project_aux(raw, axes; kwargs...) -> dest
 
 Project `raw` and append a derived auxiliary domain axis carrying its flux,
@@ -217,32 +217,32 @@ multiplet as laid out by `stack`). Like `project`, it verifies that only a
 negligible component is discarded; see [`unchecked_project_aux`](@ref) and
 [`tryproject_aux`](@ref) for the unchecked and nullable siblings.
 """
-function project_aux(raw, codomain_axes, domain_axes; kwargs...)
-    return project(project_aux_args(raw, codomain_axes, domain_axes)...; kwargs...)
+function project_aux(raw, axes_codomain, axes_domain; kwargs...)
+    return project(project_aux_args(raw, axes_codomain, axes_domain)...; kwargs...)
 end
 project_aux(raw, axes; kwargs...) = project_aux(raw, axes, (); kwargs...)
 
 """
-    unchecked_project_aux(raw, codomain_axes, domain_axes) -> dest
+    unchecked_project_aux(raw, axes_codomain, axes_domain) -> dest
     unchecked_project_aux(raw, axes) -> dest
 
 The unchecked sibling of [`project_aux`](@ref): derive and append the auxiliary
 axis, then project without verifying which components are discarded.
 """
-function unchecked_project_aux(raw, codomain_axes, domain_axes)
-    return unchecked_project(project_aux_args(raw, codomain_axes, domain_axes)...)
+function unchecked_project_aux(raw, axes_codomain, axes_domain)
+    return unchecked_project(project_aux_args(raw, axes_codomain, axes_domain)...)
 end
 unchecked_project_aux(raw, axes) = unchecked_project_aux(raw, axes, ())
 
 """
-    tryproject_aux(raw, codomain_axes, domain_axes; kwargs...) -> Union{dest, Nothing}
+    tryproject_aux(raw, axes_codomain, axes_domain; kwargs...) -> Union{dest, Nothing}
     tryproject_aux(raw, axes; kwargs...) -> Union{dest, Nothing}
 
 The nullable sibling of [`project_aux`](@ref): derive and append the auxiliary
 axis, returning `nothing` instead of throwing when more than a negligible
 component of `raw` would be discarded.
 """
-function tryproject_aux(raw, codomain_axes, domain_axes; kwargs...)
-    return tryproject(project_aux_args(raw, codomain_axes, domain_axes)...; kwargs...)
+function tryproject_aux(raw, axes_codomain, axes_domain; kwargs...)
+    return tryproject(project_aux_args(raw, axes_codomain, axes_domain)...; kwargs...)
 end
 tryproject_aux(raw, axes; kwargs...) = tryproject_aux(raw, axes, (); kwargs...)
diff --git a/src/similar_map.jl b/src/similar_map.jl
index 4cd22f54..bdb24d76 100644
--- a/src/similar_map.jl
+++ b/src/similar_map.jl
@@ -1,13 +1,13 @@
 using Random: Random, AbstractRNG
 
 """
-    similar_map(prototype, [T,] codomain_axes, domain_axes) -> M
+    similar_map(prototype, [T,] axes_codomain, axes_domain) -> M
 
-Allocate an array shaped as a linear map from `domain_axes` to `codomain_axes`
+Allocate an array shaped as a linear map from `axes_domain` to `axes_codomain`
 with element type `T` (defaulting to `eltype(prototype)`), using `prototype` to
 determine the array backend. The domain axes are given un-dualized (codomain
 facing) and stored dual, so the default is
-`similar(prototype, T, (codomain_axes..., conj.(domain_axes)...))`. `conj`
+`similar(prototype, T, (axes_codomain..., conj.(axes_domain)...))`. `conj`
 dualizes a graded axis and is a no-op on a dense axis. Backends with map-shaped
 storage (e.g. a `TensorMap`) overload this to build the codomain/domain directly.
 
@@ -24,11 +24,11 @@ julia> eltype(M), size(M)
 (Float32, (2, 3, 4, 5))
 ```
 """
-function similar_map(prototype, ::Type{T}, codomain_axes, domain_axes) where {T}
-    return similar(prototype, T, (codomain_axes..., conj.(domain_axes)...))
+function similar_map(prototype, ::Type{T}, axes_codomain, axes_domain) where {T}
+    return similar(prototype, T, (axes_codomain..., conj.(axes_domain)...))
 end
-function similar_map(prototype, codomain_axes, domain_axes)
-    return similar_map(prototype, eltype(prototype), codomain_axes, domain_axes)
+function similar_map(prototype, axes_codomain, axes_domain)
+    return similar_map(prototype, eltype(prototype), axes_codomain, axes_domain)
 end
 
 """
@@ -70,20 +70,20 @@ for (f, g) in ((:randn, :randn), (:rand, :rand))
 end
 
 """
-    zeros_map([T,] codomain_axes, domain_axes) -> M
-    ones_map([T,] codomain_axes, domain_axes) -> M
-    randn_map([rng,] [T,] codomain_axes, domain_axes) -> M
-    rand_map([rng,] [T,] codomain_axes, domain_axes) -> M
-    fill_map(v, codomain_axes, domain_axes) -> M
+    zeros_map([T,] axes_codomain, axes_domain) -> M
+    ones_map([T,] axes_codomain, axes_domain) -> M
+    randn_map([rng,] [T,] axes_codomain, axes_domain) -> M
+    rand_map([rng,] [T,] axes_codomain, axes_domain) -> M
+    fill_map(v, axes_codomain, axes_domain) -> M
 
-Construct an array shaped as a linear map from `domain_axes` to `codomain_axes`,
+Construct an array shaped as a linear map from `axes_domain` to `axes_codomain`,
 filled with zeros (`zeros_map`), ones (`ones_map`), normally-distributed values
 (`randn_map`), uniformly-distributed values (`rand_map`), or the value `v` (`fill_map`),
 with element type `T` (defaulting to `Float64`; `fill_map` takes it from `v`). These are the
 value-filling companions of [`similar_map`](@ref): the domain axes are given un-dualized
 (codomain facing) and stored dual, so the default flattens to the axis-friendly
 [`zeros`](@ref)/[`ones`](@ref)/[`randn`](@ref)/[`rand`](@ref)/[`fill`](@ref) over
-`(codomain_axes..., conj.(domain_axes)...)` (`conj` dualizes a graded axis and is a
+`(axes_codomain..., conj.(axes_domain)...)` (`conj` dualizes a graded axis and is a
 no-op on a dense one). Backends with map-shaped storage (e.g. a `TensorMap`) overload
 these to build the codomain/domain directly.
 """
@@ -97,31 +97,31 @@ function fill_map end
 @doc (@doc zeros_map) rand_map
 @doc (@doc zeros_map) fill_map
 
-zeros_map(codomain_axes, domain_axes) = zeros_map(Float64, codomain_axes, domain_axes)
-function zeros_map(::Type{T}, codomain_axes, domain_axes) where {T}
-    return zeros(T, (codomain_axes..., conj.(domain_axes)...))
+zeros_map(axes_codomain, axes_domain) = zeros_map(Float64, axes_codomain, axes_domain)
+function zeros_map(::Type{T}, axes_codomain, axes_domain) where {T}
+    return zeros(T, (axes_codomain..., conj.(axes_domain)...))
 end
-ones_map(codomain_axes, domain_axes) = ones_map(Float64, codomain_axes, domain_axes)
-function ones_map(::Type{T}, codomain_axes, domain_axes) where {T}
-    return ones(T, (codomain_axes..., conj.(domain_axes)...))
+ones_map(axes_codomain, axes_domain) = ones_map(Float64, axes_codomain, axes_domain)
+function ones_map(::Type{T}, axes_codomain, axes_domain) where {T}
+    return ones(T, (axes_codomain..., conj.(axes_domain)...))
 end
-function fill_map(value, codomain_axes, domain_axes)
-    return fill(value, (codomain_axes..., conj.(domain_axes)...))
+function fill_map(value, axes_codomain, axes_domain)
+    return fill(value, (axes_codomain..., conj.(axes_domain)...))
 end
 
 for f in (:randn_map, :rand_map)
     g = Symbol(chopsuffix(String(f), "_map"))
     @eval begin
-        $f(codomain_axes, domain_axes) =
-            $f(Random.default_rng(), codomain_axes, domain_axes)
-        function $f(rng::AbstractRNG, codomain_axes, domain_axes)
-            return $f(rng, Float64, codomain_axes, domain_axes)
+        $f(axes_codomain, axes_domain) =
+            $f(Random.default_rng(), axes_codomain, axes_domain)
+        function $f(rng::AbstractRNG, axes_codomain, axes_domain)
+            return $f(rng, Float64, axes_codomain, axes_domain)
         end
-        function $f(::Type{T}, codomain_axes, domain_axes) where {T}
-            return $f(Random.default_rng(), T, codomain_axes, domain_axes)
+        function $f(::Type{T}, axes_codomain, axes_domain) where {T}
+            return $f(Random.default_rng(), T, axes_codomain, axes_domain)
         end
-        function $f(rng::AbstractRNG, ::Type{T}, codomain_axes, domain_axes) where {T}
-            return $g(rng, T, (codomain_axes..., conj.(domain_axes)...))
+        function $f(rng::AbstractRNG, ::Type{T}, axes_codomain, axes_domain) where {T}
+            return $g(rng, T, (axes_codomain..., conj.(axes_domain)...))
         end
     end
 end
diff --git a/test/test_basics.jl b/test/test_basics.jl
index 100e65e6..9e285f86 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -156,14 +156,19 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         unmatricizeperm!(a, m, invperm_codomain, invperm_domain)
         @test a ≈ a1
 
-        a = unmatricize(m, (), axes0)
+        a = unmatricize(reshape(a0, 1, 120), (), axes0)
         @test eltype(a) === elt
         @test a ≈ a0
 
-        a = unmatricize(m, axes0, ())
+        a = unmatricize(reshape(a0, 120, 1), axes0, ())
         @test eltype(a) === elt
         @test a ≈ a0
 
+        # The rows must be the fused codomain and the columns the fused domain, so a
+        # matrix with the right element count but the wrong split is rejected.
+        @test_throws DimensionMismatch unmatricize(m, (), axes0)
+        @test_throws DimensionMismatch unmatricize(m, axes0, ())
+
         m = randn(elt, 1, 1)
         a = unmatricize(m, (), ())
         @test a isa Array{elt, 0}
@@ -173,6 +178,14 @@ TensorAlgebra.label_type(::Type{OptInLabel}) = Int
         @test_throws ArgumentError unmatricizeperm!(m, m, (1, 2), (3,))
     end
 
+    @testset "contract eltype widens like a matrix product" begin
+        a1 = ones(Bool, (2, 2))
+        a2 = ones(Bool, (2, 2))
+        a_dest, = contract(a1, ("i", "k"), a2, ("k", "j"))
+        @test eltype(a_dest) === Int
+        @test a_dest == fill(2, (2, 2))
+    end
+
     alg_tensoroperations = ContractAlgorithm(TensorOperations.StridedBLAS())
     @testset "contract (eltype1=$elt1, eltype2=$elt2)" for elt1 in elts, elt2 in elts
         elt_dest = promote_type(elt1, elt2)
diff --git a/test/test_diagonal.jl b/test/test_diagonal.jl
index 29795c6d..c492918a 100644
--- a/test/test_diagonal.jl
+++ b/test/test_diagonal.jl
@@ -1,6 +1,6 @@
 using LinearAlgebra: Diagonal, diag
 using TensorAlgebra: TensorAlgebra
-using Test: @test, @testset
+using Test: @test, @test_throws, @testset
 
 @testset "Diagonal TensorAlgebra interface (eltype=$elt)" for elt in (Float64, ComplexF64)
     d = Diagonal(elt[2, 3, 4])
@@ -45,12 +45,30 @@ using Test: @test, @testset
         @test m === d
     end
 
-    @testset "unmatricize round-trips a Diagonal" begin
+    @testset "unmatricize round-trips a Diagonal on its own {1,1} axes" begin
         ax = axes(d, 1)
         back = TensorAlgebra.unmatricize(TensorAlgebra.ReshapeMatricize(), d, (ax,), (ax,))
         @test back === d
     end
 
+    @testset "unmatricize densifies a genuine bond-split" begin
+        d4 = Diagonal(elt[1, 2, 3, 4])
+        axes_codomain = (Base.OneTo(2), Base.OneTo(2))
+        axes_domain = (Base.OneTo(4),)
+        t = TensorAlgebra.unmatricize(
+            TensorAlgebra.ReshapeMatricize(), d4, axes_codomain, axes_domain
+        )
+        @test !(t isa Diagonal)
+        @test t == reshape(Array(d4), 2, 2, 4)
+    end
+
+    @testset "unmatricize errors on a mismatched {1,1} split" begin
+        wrong = Base.OneTo(length(diag(d)) + 1)
+        @test_throws DimensionMismatch TensorAlgebra.unmatricize(
+            TensorAlgebra.ReshapeMatricize(), d, (wrong,), (wrong,)
+        )
+    end
+
     @testset "matrix functions preserve Diagonal" begin
         dp = Diagonal(elt[4, 9, 16])
         s = TensorAlgebra.sqrt(dp, ("i", "j"), ("i",), ("j",))
@@ -61,18 +79,40 @@ using Test: @test, @testset
         @test e ≈ exp(dp)
     end
 
-    @testset "contract densifies (Diagonal is an input structure, not an output one)" begin
+    @testset "contract stays Diagonal on the matmul pattern, densifies otherwise" begin
         d2 = Diagonal(elt[10, 20, 30])
-        # One contracted leg: a matrix product, materialized dense.
+        # One contracted leg: the matmul/endomorphism pattern stays Diagonal.
         c2, = TensorAlgebra.contract(d, ("i", "k"), d2, ("k", "j"))
-        @test !(c2 isa Diagonal)
+        @test c2 isa Diagonal
         @test c2 ≈ d * d2
+        # All transpose variants of the single-contracted-leg pattern stay Diagonal.
+        for (l1, l2) in (
+                (("i", "k"), ("j", "k")),
+                (("k", "i"), ("k", "j")),
+                (("k", "i"), ("j", "k")),
+            )
+            ct, = TensorAlgebra.contract(d, l1, d2, l2)
+            @test ct isa Diagonal
+            @test ct ≈ d * d2
+        end
         # Both legs contracted: a scalar.
         c0, = TensorAlgebra.contract(d, ("i", "j"), d2, ("i", "j"))
         @test ndims(c0) == 0
         @test c0[] ≈ sum(diag(d) .* diag(d2))
-        # No contracted legs: a rank-4 outer product.
+        # No contracted legs: a rank-4 outer product, densified.
         c4, = TensorAlgebra.contract(d, ("i", "j"), d2, ("k", "l"))
+        @test !(c4 isa Diagonal)
         @test ndims(c4) == 4
+        # Diagonal times dense: densifies.
+        a = reshape(elt[1:9;], 3, 3)
+        cda, = TensorAlgebra.contract(d, ("i", "k"), a, ("k", "j"))
+        @test !(cda isa Diagonal)
+        @test cda ≈ d * a
+        # Mixed element types promote through the product.
+        dc = Diagonal(ComplexF64[1 + im, 2, 3im])
+        cmix, = TensorAlgebra.contract(d, ("i", "k"), dc, ("k", "j"))
+        @test cmix isa Diagonal
+        @test eltype(cmix) === ComplexF64
+        @test cmix ≈ d * dc
     end
 end
diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl
index 7e4673f7..a933ee0a 100644
--- a/test/test_factorizations.jl
+++ b/test/test_factorizations.jl
@@ -90,7 +90,7 @@ end
     D, V = eig_full(A, labels_A, labels_V, labels_V′)
     @test A == Acopy # should not have altered initial array
     @test eltype(D) == eltype(V) && eltype(D) <: Complex
-    # `D` is unmatricized like the other factors; the `(1, 1)` unmatricize preserves `Diagonal`.
+    # `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))
diff --git a/test/test_similar_map.jl b/test/test_similar_map.jl
index 7446b6e2..009dbd10 100644
--- a/test/test_similar_map.jl
+++ b/test/test_similar_map.jl
@@ -46,7 +46,7 @@ end
     end
 end
 
-# The dense map constructors flatten `(codomain_axes..., conj.(domain_axes)...)`; `conj` is a
+# The dense map constructors flatten `(axes_codomain..., conj.(axes_domain)...)`; `conj` is a
 # no-op on a dense axis, so the shape is the concatenation of codomain and domain lengths.
 @testset "map construction ($T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
     cod = (Base.OneTo(2), Base.OneTo(3))
diff --git a/test/test_tensorkitext.jl b/test/test_tensorkitext.jl
index 74d223a7..557b4983 100644
--- a/test/test_tensorkitext.jl
+++ b/test/test_tensorkitext.jl
@@ -56,13 +56,13 @@ using Test: @test, @test_throws, @testset
         B = Rep[U₁](0 => 1, -1 => 2)
         C1 = Rep[U₁](0 => 2)
         t = randn(rng, elt, A1 ⊗ A2, B ⊗ C1)
-        codomain_axes = (space(t, 1), space(t, 2))
+        axes_codomain = (space(t, 1), space(t, 2))
         # `unmatricize` takes the domain axes codomain-facing (un-dualized), so pass `B`, `C1`
         # directly rather than the dualized `space(t, 3)`, `space(t, 4)`.
-        domain_axes = (B, C1)
+        axes_domain = (B, C1)
         m = matricize(t, Val(2))
         @test space(m) == space(t)
-        back = unmatricize(m, codomain_axes, domain_axes)
+        back = unmatricize(m, axes_codomain, axes_domain)
         @test back ≈ t
     end