further progress on DiagonalTensorMap [skip ci]

Author: Jutho

src/spaces/homspace.jl

@@ -129,6 +129,9 @@ Return the `HomSpace` obtained by permuting the indices of the domain and codoma
 according to the permutation `p₁` and `p₂` respectively.
 """
 function permute(W::HomSpace{S}, (p₁, p₂)::Index2Tuple{N₁,N₂}) where {S,N₁,N₂}
+    p = (p₁..., p₂...)
+    TupleTools.isperm(p) && length(p) == numind(W) ||
+        throw(ArgumentError("$((p₁, p₂)) is not a valid permutation for $(W)"))
     cod = ProductSpace{S,N₁}(map(n -> W[n], p₁))
     dom = ProductSpace{S,N₂}(map(n -> dual(W[n]), p₂))
     return cod ← dom

src/tensors/diagtensor.jl

@@ -14,14 +14,14 @@ struct DiagonalTensorMap{T,S<:IndexSpace,A<:DenseVector{T}} <: AbstractTensorMap
     function DiagonalTensorMap{T,S,A}(data::A,
                                       dom::S) where {T,S<:IndexSpace,A<:DenseVector{T}}
         T ⊆ field(S) || @warn("scalartype(data) = $T ⊈ $(field(S)))", maxlog = 1)
-        return DiagonalTensorMap{T,S,A}(data, dom)
+        return new{T,S,A}(data, dom)
     end
 end
 reduceddim(V::IndexSpace) = sum(c -> dim(V, c), sectors(V); init=0)
 
 # Basic methods for characterising a tensor:
 #--------------------------------------------
-space(t::DiagonalTensorMap) = t.domain ← t.domain
+space(d::DiagonalTensorMap) = d.domain ← d.domain
 
 """
     storagetype(::Union{T,Type{T}}) where {T<:TensorMap} -> Type{A<:DenseVector}
@@ -58,36 +58,49 @@ end
 
 # Special case adjoint:
 #-----------------------
-Base.adjoint(t::DiagonalTensorMap{<:Real}) = t
-Base.adjoint(t::DiagonalTensorMap{<:Complex}) = DiagonalTensorMap(conj(t.data), t.domain)
+Base.adjoint(d::DiagonalTensorMap{<:Real}) = d
+Base.adjoint(d::DiagonalTensorMap{<:Complex}) = DiagonalTensorMap(conj(d.data), d.domain)
 
-# Efficient copy constructors
-#-----------------------------
-Base.copy(t::DiagonalTensorMap) = typeof(t)(copy(t.data), t.domain)
+# Efficient copy constructors and TensorMap converters
+#-----------------------------------------------------
+Base.copy(d::DiagonalTensorMap) = typeof(d)(copy(d.data), d.domain)
 
-function Base.complex(t::DiagonalTensorMap)
-    if scalartype(t) <: Complex
-        return t
-    else
-        return DiagonalTensorMap(complex(t.data), t.domain)
+function Base.copy!(t::AbstractTensorMap, d::DiagonalTensorMap)
+    space(t) == space(d) || throw(SpaceMismatch())
+    for (c, b) in blocks(d)
+        copy!(block(t, c), b)
+    end
+    return t
+end
+TensorMap(d::DiagonalTensorMap) = copy!(similar(d), d)
+Base.convert(::Type{TensorMap}, d::DiagonalTensorMap) = TensorMap(d)
+
+# Complex, real and imaginary parts
+#-----------------------------------
+for f in (:real, :imag, :complex)
+    @eval begin
+        function Base.$f(d::DiagonalTensorMap)
+            return DiagonalTensorMap($f(d.data), d.domain)
+        end
     end
 end
 
 # Getting and setting the data at the block level
 #-------------------------------------------------
-blocksectors(t::DiagonalTensorMap) = blocksectors(t.domain)
+blocksectors(d::DiagonalTensorMap) = blocksectors(d.domain)
 
-function block(t::DiagonalTensorMap, s::Sector)
-    sectortype(t) == typeof(s) || throw(SectorMismatch())
+function block(d::DiagonalTensorMap, s::Sector)
+    sectortype(d) == typeof(s) || throw(SectorMismatch())
     offset = 0
-    for c in sectors(t)
+    dom = domain(d)[1]
+    for c in sectors(dom)
         if c < s
-            offset += dim(t, c)
+            offset += dim(dom, c)
         elseif c == s
-            r = offset .+ (1:dim(t, c))
-            return Diagonal(view(t.data, r))
+            r = offset .+ (1:dim(dom, c))
+            return Diagonal(view(d.data, r))
         else # s not in sectors(t)
-            return Diagonal(view(t.data, 1:0))
+            return Diagonal(view(d.data, 1:0))
         end
     end
 end
@@ -96,18 +109,130 @@ end
 
 # Indexing and getting and setting the data at the subblock level
 #-----------------------------------------------------------------
-@inline function Base.getindex(t::DiagonalTensorMap,
+@inline function Base.getindex(d::DiagonalTensorMap,
                                f₁::FusionTree{I,1},
                                f₂::FusionTree{I,1}) where {I<:Sector}
     s = f₁.uncoupled[1]
-    s == f₁.uncoulped == f₂.uncoupled[1] == f₂.uncoupled || throw(SectorMismatch())
-    return block(t, s)
+    s == f₁.coupled == f₂.uncoupled[1] == f₂.coupled || throw(SectorMismatch())
+    return block(d, s)
     # TODO: do we want a StridedView here? Then we need to allocate a new matrix.
 end
 
-function Base.setindex!(t::TensorMap,
+function Base.setindex!(d::DiagonalTensorMap,
                         v,
                         f₁::FusionTree{I,1},
                         f₂::FusionTree{I,1}) where {I<:Sector}
-    return copy!(getindex(t, f₁, f₂), v)
+    return copy!(getindex(d, f₁, f₂), v)
+end
+
+function Base.getindex(d::DiagonalTensorMap)
+    sectortype(d) === Trivial || throw(SectorMismatch())
+    return Diagonal(d.data)
+end
+
+# Index manipulations
+# -------------------
+function has_shared_permute(d::DiagonalTensorMap, (p₁, p₂)::Index2Tuple)
+    if p₁ === codomainind(d) && p₂ === domainind(d)
+        return true
+    elseif BraidingStyle(sectortype(d)) isa Bosonic
+        # TODO: is this always correct? transpose has no effect for bosonic sectors?
+        return p₁ === domainind(d) && p₂ === codomainind(d)
+    else
+        return false
+    end
+end
+
+function permute(d::DiagonalTensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂};
+                 copy::Bool=false) where {N₁,N₂}
+    if !copy
+        if p₁ === codomainind(d) && p₂ === domainind(d)
+            return d
+        elseif has_shared_permute(d, (p₁, p₂)) # tranpose for bosonic sectors
+            return DiagonalTensorMap(d.data, dual(d.domain))
+        end
+    end
+    # general case
+    space′ = permute(space(d), (p₁, p₂))
+    @inbounds begin
+        return permute!(similar(d, space′), d, (p₁, p₂))
+    end
+end
+
+# VectorInterface
+# ---------------
+function VectorInterface.zerovector(d::DiagonalTensorMap, ::Type{S}) where {S<:Number}
+    return DiagonalTensorMap(zerovector(d.data, S), d.domain)
 end
+function VectorInterface.add(ty::DiagonalTensorMap, tx::DiagonalTensorMap,
+                             α::Number, β::Number)
+    domain(ty) == domain(tx) || throw(SpaceMismatch("$(space(ty)) ≠ $(space(tx))"))
+    T = VectorInterface.promote_add(ty, tx, α, β)
+    return add!(scale!(zerovector(ty, T), ty, β), tx, α) # zerovector instead of similar preserves diagonal structure
+end
+
+# Linear Algebra and factorizations
+# ---------------------------------
+function one!(d::DiagonalTensorMap)
+    fill!(d.data, one(eltype(d.data)))
+    return d
+end
+function Base.one(d::DiagonalTensorMap)
+    return DiagonalTensorMap(one.(d.data), d.domain)
+end
+function Base.zero(d::DiagonalTensorMap)
+    return DiagonalTensorMap(zero.(d.data), d.domain)
+end
+
+function eig!(d::DiagonalTensorMap)
+    return d, one(d)
+end
+function eigh!(d::DiagonalTensorMap{<:Real})
+    return d, one(d)
+end
+function eigh!(d::DiagonalTensorMap{<:Complex})
+    # TODO: should this test for hermiticity? `eigh!(::TensorMap)` also does not do this.
+    return DiagonalTensorMap(real(d.data), d.domain), one(d)
+end
+
+function leftorth!(d::DiagonalTensorMap; kwargs...)
+    return one(d), d # TODO: this is only correct for `alg = QR()` or `alg = QL()`
+end
+function rightorth!(d::DiagonalTensorMap; kwargs...)
+    return d, one(d) # TODO: this is only correct for `alg = LQ()` or `alg = RQ()`
+end
+# not much to do here:
+leftnull!(d::DiagonalTensorMap; kwargs...) = leftnull!(TensorMap(d); kwargs...)
+rightnull!(d::DiagonalTensorMap; kwargs...) = rightnull!(TensorMap(d); kwargs...)
+
+function tsvd!(d::DiagonalTensorMap; trunc=NoTruncation(), p::Real=2, alg=SDD())
+    return _tsvd!(d, alg, trunc, p)
+end
+# helper function
+function _compute_svddata!(d::DiagonalTensorMap, alg::Union{SVD,SDD})
+    InnerProductStyle(d) === EuclideanInnerProduct() || throw_invalid_innerproduct(:tsvd!)
+    I = sectortype(d)
+    dims = SectorDict{I,Int}()
+    generator = Base.Iterators.map(blocks(d)) do (c, b)
+        lb = length(b.diag)
+        U = zerovector!(similar(b.diag, lb, lb))
+        V = zerovector!(similar(b.diag, lb, lb))
+        p = sortperm(b.diag; by=abs, rev=true)
+        for (i, pi) in enumerate(p)
+            U[pi, i] = MatrixAlgebra.safesign(b.diag[pi])
+            V[i, pi] = 1
+        end
+        Σ = abs.(view(b.diag, p))
+        dims[c] = lb
+        return c => (U, Σ, V)
+    end
+    SVDdata = SectorDict(generator)
+    return SVDdata, dims
+end
+
+# matrix functions
+for f in
+    (:exp, :cos, :sin, :tan, :cot, :cosh, :sinh, :tanh, :coth, :atan, :acot, :asinh, :sqrt,
+     :log, :asin, :acos, :acosh, :atanh, :acoth)
+    @eval Base.$f(d::DiagonalTensorMap) = DiagonalTensorMap($f.(d.data), d.domain)
+end

src/tensors/indexmanipulations.jl

@@ -31,12 +31,23 @@ function permute(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂};
                  copy::Bool=false) where {N₁,N₂}
     space′ = permute(space(t), (p₁, p₂))
     # share data if possible
+    if !copy && p₁ === codomainind(t) && p₂ === domainind(t)
+        return t
+    end
+    # general case
+    @inbounds begin
+        return permute!(similar(t, space′), t, (p₁, p₂))
+    end
+end
+function permute(t::TensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂}; copy::Bool=false) where {N₁,N₂}
+    @show (p₁, p₂)
+    space′ = permute(space(t), (p₁, p₂))
+    # share data if possible
     if !copy
         if p₁ === codomainind(t) && p₂ === domainind(t)
             return t
         elseif has_shared_permute(t, (p₁, p₂))
-            return TensorMap(reshape(t.data, dim(codomain(space′)), dim(domain(space′))),
-                             codomain(space′), domain(space′))
+            return TensorMap(t.data, space′)
         end
     end
     # general case
@@ -53,6 +64,9 @@ function permute(t::AbstractTensorMap, p::IndexTuple; copy::Bool=false)
     return permute(t, (p, ()); copy)
 end
 
+function has_shared_permute(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple)
+    return (p₁ === codomainind(t) && p₂ === domainind(t))
+end
 function has_shared_permute(t::TensorMap, (p₁, p₂)::Index2Tuple)
     if p₁ === codomainind(t) && p₂ === domainind(t)
         return true

src/tensors/vectorinterface.jl

@@ -23,7 +23,7 @@ VectorInterface.zerovector!!(t::AbstractTensorMap) = zerovector!(t)
 #-------------------------
 function VectorInterface.scale(t::AbstractTensorMap, α::Number)
     T = VectorInterface.promote_scale(t, α)
-    return scale!(similar(t, T), t, α)
+    return scale!(zerovector(t, T), t, α)
 end
 function VectorInterface.scale!(t::AbstractTensorMap, α::Number)
     for (c, b) in blocks(t)
@@ -74,7 +74,7 @@ function VectorInterface.add!(ty::AbstractTensorMap, tx::AbstractTensorMap,
                               α::Number, β::Number)
     space(ty) == space(tx) || throw(SpaceMismatch("$(space(ty)) ≠ $(space(tx))"))
     for ((cy, by), (cx, bx)) in zip(blocks(ty), blocks(tx))
-        add!(StridedView(by), StridedView(bx), α, β)
+        add!(by, bx, α, β)
     end
     return ty
 end