fix exponentiation and add test, rename tensor concatenation

Author: Jutho

src/tensors/linalg.jl

@@ -12,15 +12,19 @@ function Base.:-(t1::AbstractTensorMap, t2::AbstractTensorMap)
     return axpy!(-one(T), t2, copyto!(similar(t1, T), t1))
 end
 
-Base.:*(t::AbstractTensorMap, α::Number) = mul!(similar(t, promote_type(eltype(t), typeof(α))), t, α)
-Base.:*(α::Number, t::AbstractTensorMap) = mul!(similar(t, promote_type(eltype(t), typeof(α))), α, t)
+Base.:*(t::AbstractTensorMap, α::Number) =
+    mul!(similar(t, promote_type(eltype(t), typeof(α))), t, α)
+Base.:*(α::Number, t::AbstractTensorMap) =
+    mul!(similar(t, promote_type(eltype(t), typeof(α))), α, t)
 Base.:/(t::AbstractTensorMap, α::Number) = *(t, one(α)/α)
 Base.:\(α::Number, t::AbstractTensorMap) = *(t, one(α)/α)
 
 LinearAlgebra.normalize!(t::AbstractTensorMap, p::Real = 2) = rmul!(t, inv(norm(t, p)))
-LinearAlgebra.normalize(t::AbstractTensorMap, p::Real = 2) = mul!(similar(t), t, inv(norm(t, p)))
+LinearAlgebra.normalize(t::AbstractTensorMap, p::Real = 2) =
+    mul!(similar(t), t, inv(norm(t, p)))
 
-Base.:*(t1::AbstractTensorMap, t2::AbstractTensorMap) = mul!(similar(t1, promote_type(eltype(t1),eltype(t2)), codomain(t1)←domain(t2)), t1, t2)
+Base.:*(t1::AbstractTensorMap, t2::AbstractTensorMap) =
+    mul!(similar(t1, promote_type(eltype(t1),eltype(t2)), codomain(t1)←domain(t2)), t1, t2)
 Base.exp(t::AbstractTensorMap) = exp!(copy(t))
 
 # Special purpose constructors
@@ -159,15 +163,26 @@ end
 
 # TensorMap exponentation:
 function exp!(t::TensorMap)
-    domain(t) == codomain(t) || error("Exponentional of a tensor only exist when domain == codomain.")
+    domain(t) == codomain(t) ||
+        error("Exponentional of a tensor only exist when domain == codomain.")
     for (c,b) in blocks(t)
-        copyto!(b, exp!(b))
+        copyto!(b, LinearAlgebra.exp!(b))
     end
     return t
 end
 
-# hcat and vcat of tensors
-function Base.hcat(t1::AbstractTensorMap{S,N₁,1}, t2::AbstractTensorMap{S,N₁,1}) where {S,N₁}
+# # concatenate tensors
+# function concatenate(t1::AbstractTensorMap{S}, ts::AbstractTensorMap{S}...; direction::Symbol) where {S}
+#     if direction = :domain
+#         numin(t1) == 1 || throw(SpaceMismatch("concatenation along domain"))
+#         for t2 in ts
+#             domain(t1) == domain(t2) || throw(SpaceMismatch())
+#         end
+#         t = t1
+#         for t2 in ts
+#             cat(t1)
+
+function catdomain(t1::AbstractTensorMap{S,N₁,1}, t2::AbstractTensorMap{S,N₁,1}) where {S,N₁}
     codomain(t1) == codomain(t2) || throw(SpaceMismatch())
 
     V1, = domain(t1)
@@ -182,7 +197,7 @@ function Base.hcat(t1::AbstractTensorMap{S,N₁,1}, t2::AbstractTensorMap{S,N₁
     end
     return t
 end
-function Base.vcat(t1::AbstractTensorMap{S,1,N₂}, t2::AbstractTensorMap{S,1,N₂}) where {S,N₂}
+function catcodomain(t1::AbstractTensorMap{S,1,N₂}, t2::AbstractTensorMap{S,1,N₂}) where {S,N₂}
     domain(t1) == domain(t2) || throw(SpaceMismatch())
 
     V1, = codomain(t1)

src/tensors/tensoroperations.jl

@@ -131,7 +131,7 @@ function add!(α, tsrc::AbstractTensorMap{S}, β, tdst::AbstractTensorMap{S,N₁
                 end
             end
         else # debugging is easier this way
-            @inbounds for (f1,f2) in fusiontrees(tsrc)
+            for (f1,f2) in fusiontrees(tsrc)
                 _addabelianblock!(α, tsrc, β, tdst, p1, p2, f1, f2)
             end
         end
@@ -145,46 +145,49 @@ function add!(α, tsrc::AbstractTensorMap{S}, β, tdst::AbstractTensorMap{S,N₁
         elseif β != 1
             mul!(tdst, β, tdst)
         end
-        @inbounds for (f1,f2) in fusiontrees(tsrc)
+        for (f1,f2) in fusiontrees(tsrc)
             for ((f1′,f2′), coeff) in permute(f1, f2, p1, p2)
-                for i in p2
+                @inbounds for i in p2
                     if i <= n && !isdual(cod[i])
                         b = f1.uncoupled[i]
                         coeff *= frobeniusschur(b) #*fermionparity(b)
                     end
                 end
-                for i in p1
+                @inbounds for i in p1
                     if i > n && isdual(dom[i-n])
                         b = f2.uncoupled[i-n]
                         coeff /= frobeniusschur(b) #*fermionparity(b)
                     end
                 end
-                axpy!(α*coeff, permutedims(tsrc[f1,f2], pdata), tdst[f1′,f2′])
+                @inbounds axpy!(α*coeff, permutedims(tsrc[f1,f2], pdata), tdst[f1′,f2′])
             end
         end
     end
     return tdst
 end
 
-@inbounds function _addabelianblock!(α, tsrc::AbstractTensorMap{S}, β, tdst::AbstractTensorMap{S,N₁,N₂}, p1::IndexTuple{N₁}, p2::IndexTuple{N₂}, f1::FusionTree, f2::FusionTree)  where {S,N₁,N₂}
+function _addabelianblock!(α, tsrc::AbstractTensorMap{S},
+                            β, tdst::AbstractTensorMap{S,N₁,N₂},
+                            p1::IndexTuple{N₁}, p2::IndexTuple{N₂},
+                            f1::FusionTree, f2::FusionTree) where {S,N₁,N₂}
     cod = codomain(tsrc)
     dom = domain(tsrc)
     n = length(cod)
     (f1′,f2′), coeff = first(permute(f1, f2, p1, p2))
-    for i in p2
+    @inbounds for i in p2
         if i <= n && !isdual(cod[i])
             b = f1.uncoupled[i]
             coeff *= frobeniusschur(b) #*fermionparity(b)
         end
     end
-    for i in p1
+    @inbounds for i in p1
         if i > n && isdual(dom[i-n])
             b = f2.uncoupled[i-n]
             coeff /= frobeniusschur(b) #*fermionparity(b)
         end
     end
     pdata = (p1...,p2...)
-    axpby!(α*coeff, permutedims(tsrc[f1,f2], pdata), β, tdst[f1′,f2′])
+    @inbounds axpby!(α*coeff, permutedims(tsrc[f1,f2], pdata), β, tdst[f1′,f2′])
 end
 
 

test/tensors.jl

@@ -159,4 +159,14 @@ VSU₂ = (ℂ[SU₂](0=>1, 1//2=>1, 1=>2),
             end
         end
     end
+    @testset "Exponentiation" begin
+        W = V1 ⊗ V2 ⊗ V3
+        for T in (Float32, Float64, ComplexF32, ComplexF64)
+            t = TensorMap(rand, T, W, W)
+            s = dim(W)
+            expt = @inferred exp(t)
+            @test reshape(convert(Array, expt), (s,s)) ≈
+                    exp(reshape(convert(Array, t), (s,s)))
+        end
+    end
 end