Fix named tensor inner products

Project.toml

@@ -1,6 +1,6 @@
 name = "ITensorBase"
 uuid = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
-version = "0.16.1"
+version = "0.16.2"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
@@ -34,13 +34,13 @@ ITensorBaseTensorKitExt = "TensorKit"
 Accessors = "0.1.39"
 Adapt = "4.1.1"
 ArrayLayouts = "1.11"
-GradedArrays = "0.17"
+GradedArrays = "0.17.3"
 LinearAlgebra = "1.10"
 MatrixAlgebraKit = "0.2, 0.3, 0.4, 0.5, 0.6"
 Mooncake = "0.4.202, 0.5"
 Random = "1.10"
 SimpleTraits = "0.9.4"
-TensorAlgebra = "0.23.1"
+TensorAlgebra = "0.23.2"
 TensorKit = "0.17"
 TensorKitSectors = "0.3.9"
 UUIDs = "1.10"

src/linearalgebra.jl

@@ -1,4 +1,5 @@
 using LinearAlgebra: LinearAlgebra as LA
+using TensorAlgebra: TensorAlgebra as TA
 
 # We overload `LinearAlgebra.norm` because the LinearAlgebra.jl AbstractArray definition
 # uses scalar indexing:
@@ -29,10 +30,8 @@ for f! in [:mul!, :div!]
     end
 end
 
-# We overload `LienarAlgebra.dot` because the LinearAlgebra.jl AbstractArray definition
-# uses scalar indexing:
-# https://github.com/JuliaLang/LinearAlgebra.jl/blob/3a4fdad7f608928ecb4b41e76b1e9ecacd058444/src/generic.jl#L919-L1009
-# which isn't friendly for named arrays wrapping GPU arrays.
+# The Hilbert–Schmidt pairing is tr(a1' * a2) after aligning the matrix representations.
+# Contracting conj(a1) with a2 inserts fermionic parity signs on dual legs.
 function LA.dot(a1::AbstractNamedTensor, a2::AbstractNamedTensor)
-    return (conj(a1) * a2)[]
+    return TA.dot(unnamed(a1), names(a1), unnamed(a2), names(a2))
 end

test/test_linearalgebra.jl

@@ -1,6 +1,11 @@
 import LinearAlgebra as LA
-using ITensorBase: ITensorBase, Named, unname, unnamed
+using GradedArrays: SU2, Z2, fSU2, fZ2, gradedrange
+using ITensorBase: ITensorBase, Named, NamedTensor, unname, unnamed
+using StableRNGs: StableRNG
+using TensorAlgebra: bipermutedims
+using TensorKit: TensorKit
 using Test: @test, @testset
+using VectorInterface: VectorInterface as VI
 
 @testset "LinearAlgebra (eltype=$(elt))" for elt in
     (Float32, Float64, Complex{Float32})
@@ -16,3 +21,37 @@ using Test: @test, @testset
     @test unnamed(LA.ldiv!(2, copy(a))) ≈ 2 \ unnamed(a)
     @test LA.dot(a, b) ≈ LA.dot(unnamed(a), unname(b, ITensorBase.names(a)))
 end
+
+@testset "Graded inner product ($G, $T, split=$n)" for (G, g) in (
+            ("Z2", gradedrange([Z2(0) => 2, Z2(1) => 1])),
+            ("fZ2", gradedrange([fZ2(false) => 2, fZ2(true) => 1])),
+            ("SU2", gradedrange([SU2(0) => 2, SU2(1 // 2) => 1, SU2(1) => 1])),
+            ("fSU2", gradedrange([fSU2(0) => 2, fSU2(1 // 2) => 1, fSU2(1) => 1])),
+        ),
+        T in (Float64, ComplexF64),
+        n in 0:3
+
+    rng = StableRNG(123)
+    cod = ntuple(_ -> g, n)
+    dom = ntuple(_ -> g, 3 - n)
+    x = randn(rng, T, cod, dom)
+    y = randn(rng, T, cod, dom)
+    a = NamedTensor(x, (:i, :j, :k))
+    b = NamedTensor(y, (:i, :j, :k))
+    expected = LA.dot(TensorKit.TensorMap(x), TensorKit.TensorMap(y))
+    @test LA.dot(a, a) ≈ LA.norm(a)^2
+    @test LA.dot(a, b) ≈ expected
+    @test VI.inner(a, b) ≈ expected
+    @test LA.dot(a, b) ≈ conj(LA.dot(b, a))
+    for m in 0:3
+        repartitioned =
+            NamedTensor(bipermutedims(y, Tuple(1:m), Tuple((m + 1):3)), (:i, :j, :k))
+        @test repartitioned ≈ b
+        @test LA.dot(a, repartitioned) ≈ expected
+        perm = (3, 1, 2)
+        reordered = NamedTensor(bipermutedims(y, perm[1:m], perm[(m + 1):3]), (:k, :i, :j))
+        @test reordered ≈ b
+        @test LA.dot(a, reordered) ≈ expected
+        @test LA.dot(reordered, a) ≈ conj(expected)
+    end
+end