Add named methods for the Hermitian-safe matrix functions

Project.toml

@@ -1,6 +1,6 @@
 name = "ITensorBase"
 uuid = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
-version = "0.10.8"
+version = "0.10.9"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
@@ -54,7 +54,7 @@ OMEinsumContractionOrders = "1"
 OrderedCollections = "1.6"
 Random = "1.10"
 SimpleTraits = "0.9.4"
-TensorAlgebra = "0.16.5"
+TensorAlgebra = "0.17"
 TensorKit = "0.17"
 TensorOperations = "5.3.1"
 TermInterface = "2"

docs/Project.toml

@@ -16,5 +16,5 @@ ITensorBase = "0.10"
 ITensorFormatter = "0.2.27"
 Literate = "2"
 MatrixAlgebraKit = "0.2, 0.3, 0.4, 0.5, 0.6"
-TensorAlgebra = "0.16"
+TensorAlgebra = "0.17"
 Test = "1.10"

src/abstractnamedtensor.jl

@@ -218,11 +218,14 @@ Base.zero(a::AbstractNamedTensor) = nameddimsof(a, zero(unnamed(a)))
 # `AbstractArray` fallback did through `axes`).
 Base.CartesianIndices(a::AbstractNamedTensor) = CartesianIndices(axes(a))
 
-# Forward `conj` to the underlying so that graded axes flip their sector arrows.
-# The default `AbstractArray` fallback would broadcast `conj` over elements without
-# touching the axes, which silently changes the contraction convention for tensors
-# with graded (dual-tagged) axes.
-Base.conj(a::AbstractNamedTensor) = nameddimsof(a, conj(unnamed(a)))
+# Eager conjugation reuses the linear-combination broadcast machinery: `conj.(a)` lowers to a
+# `ConjBroadcasted` leaf that the op primitives absorb (folding `conj` into their `op` flag) and
+# materializes into an all-codomain destination with dualized axes. This keeps one conjugation
+# mechanism shared by the lazy `conj.` path and eager `conj`, so a non-`AbstractArray` backend
+# (a `TensorMap`) needs no separate eager hook. The default `AbstractArray` `conj` would instead
+# map over elements without dualizing the axes, silently changing the contraction convention for
+# graded tensors.
+Base.conj(a::AbstractNamedTensor) = conj.(a)
 
 # `LinearAlgebra.normalize` infers result eltype via `typeof(first(a)/nrm)`, which
 # scalar-indexes block-structured storage. `a / norm(a, p)` already preserves names.
@@ -399,14 +402,17 @@ end
 Base.axes(a::AbstractNamedTensor, dimname::Name) = axes(a, dim(a, dimname))
 Base.size(a::AbstractNamedTensor, dimname::Name) = size(a, dim(a, dimname))
 
+# Lowered through `TensorAlgebra.similar_map` (all-codomain, so identical to
+# `similar(parent, elt, axes)` for dense) so non-`AbstractArray` backends whose `similar`
+# wants a map-shaped space (e.g. a `TensorMap`) allocate through their own overload.
 function similar_nameddims(a::AbstractNamedTensor, elt::Type, ax)
     return nameddims(
-        similar(unnamed(a), elt, unnamed.(Tuple(ax))),
+        TensorAlgebra.similar_map(unnamed(a), elt, unnamed.(Tuple(ax)), ()),
         name.(ax)
     )
 end
 function similar_nameddims(a::AbstractArray, elt::Type, ax)
-    return nameddims(similar(a, elt, unnamed.(Tuple(ax))), name.(ax))
+    return nameddims(TensorAlgebra.similar_map(a, elt, unnamed.(Tuple(ax)), ()), name.(ax))
 end
 
 # Base.similar gets the eltype at compile time.
@@ -715,8 +721,10 @@ function Base.getindex(
     )
     return getindex(a, to_indices(a, (I1, Irest...))...)
 end
+# Routed through `TensorAlgebra.scalar` (default `a[]`) so non-`AbstractArray` backends
+# without trivial-sector scalar indexing (e.g. a `TensorMap`) can overload it.
 function Base.getindex(a::AbstractNamedTensor)
-    return getindex(unnamed(a))
+    return TensorAlgebra.scalar(unnamed(a))
 end
 # Linear indexing.
 function Base.getindex(a::AbstractNamedTensor, I::Int)
@@ -1191,9 +1199,10 @@ end
 # `sum` is routed to the underlying data rather than left to fall back on the
 # `mapreduce` method above because some array types (such as graded arrays) define
 # `Base.sum` directly but not the general `mapreduce`, so the unnamed `sum` is the
-# path that works for them.
+# path that works for them. Routed through `TensorAlgebra.sum` so non-iterable
+# backends (e.g. a `TensorMap`) can overload it.
 function Base.sum(a::AbstractNamedTensor; kwargs...)
-    return sum(unnamed(a); kwargs...)
+    return TensorAlgebra.sum(unnamed(a); kwargs...)
 end
 
 function LinearAlgebra.promote_leaf_eltypes(a::AbstractNamedTensor)

src/broadcast.jl

@@ -27,18 +27,19 @@ function broadcasted_unnamed(a::AbstractNamedTensor, names)
     # common case for the rest) needs no permutation, avoiding a `getperm` allocation and the
     # identity `permuteddims` wrapper. Skipping it makes a small add several times slower.
     dimnames(a) == names && return unnamed(a)
-    return _broadcast_permuteddims_to(unnamed(a), getperm(dimnames(a), names))
+    return _broadcast_permuteddims(unnamed(a), getperm(dimnames(a), names))
 end
 # Broadcasting-only alignment: unlike the public `unnamed(a, names)` (which returns a
 # `Base.PermutedDimsArray`, a full array), this wraps in `TensorAlgebra.PermutedDims`, which stores
 # the permutation in a field rather than a type parameter, so it builds cheaply and type-stably
 # from the runtime permutation and is a broadcast leaf the linear-combination fold absorbs via
 # `bipermutedimsopadd!`. `PermutedDims` has almost no array interface, so it stays confined to this
 # hot path and is never handed back to users. Function barrier: `unnamed(a)` is abstractly typed,
-# so dispatching on the concrete array makes `ndims` a compile-time constant for the inferrable
-# `ntuple(…, Val(ndims))` permutation.
-@noinline function _broadcast_permuteddims_to(array::AbstractArray, perm)
-    return TA.PermutedDims(array, ntuple(i -> perm[i], Val(ndims(array))))
+# so dispatching on the concrete array makes the rank a compile-time constant for the inferrable
+# `ntuple(…, Val(ndims))` permutation. The rank comes from `TensorAlgebra.ndims`, which also
+# covers non-`AbstractArray` backends like a `TensorMap`.
+@noinline function _broadcast_permuteddims(array, perm)
+    return TA.PermutedDims(array, ntuple(i -> perm[i], Val(TA.ndims(array))))
 end
 function broadcasted_unnamed(bc::Broadcasted, names)
     return broadcasted(bc.f, Base.Fix2(broadcasted_unnamed, names).(bc.args)...)
@@ -85,10 +86,17 @@ end
 # call re-specializes on the concrete runtime types and everything below is type-stable
 # (`eltype(lb)` is now inferrable, no runtime `promote_op`). Inlining the body into `copy`
 # instead costs one extra allocation per call.
+#
+# Allocate from `axes(lb)`, the flattened expression's axes, rather than the prototype's own:
+# an axis-changing operand (a `conj` leaf dualizes its axes) makes them differ, and the
+# destination must match the expression. All axes go in the codomain (empty domain), the
+# all-codomain output convention `@tensor` uses for an unbipartitioned left-hand side; on a
+# non-bipartitioned backend (a dense array) `similar_map` with an empty domain is a plain
+# `similar` over `axes(lb)`.
 function _copy_unnamed(bc_unnamed, prototype)
     lb = TA.tryflattenlinear(bc_unnamed)
     isnothing(lb) && return copy(bc_unnamed)
-    return copyto!(similar(prototype, eltype(lb)), lb)
+    return copyto!(TA.similar_map(prototype, eltype(lb), axes(lb), ()), lb)
 end
 
 # `Base.Broadcast.materialize!` otherwise reconstructs the broadcast over `axes(dest)` and

src/index.jl

@@ -33,8 +33,12 @@ function IndexName(
     )
     return IndexName(id, to_tags(tags), plev)
 end
-function uniquename(rng::AbstractRNG, n::IndexName)
-    return setid(n, uuid4(rng))
+# `uniquename` on an existing `IndexName` keeps its tags and prime level, minting only a
+# fresh id (the legacy `sim`). The type form drops them: a factorization bond or a fresh
+# operator leg has no relationship to any seed's decoration, so its callers pass the name
+# type to opt out of inheriting it.
+function uniquename(rng::AbstractRNG, name::IndexName)
+    return IndexName(rng; tags = tags_stored(name), plev = plev(name))
 end
 function uniquename(rng::AbstractRNG, ::Type{<:IndexName})
     return IndexName(rng)