some simplifications and extensions

Author: Jutho

src/factorizations/matrixalgebrakit.jl

@@ -81,7 +81,8 @@ end
 
 function MAK.initialize_output(::typeof(svd_vals!), t::AbstractTensorMap, alg::AbstractAlgorithm)
     V_cod = infimum(fuse(codomain(t)), fuse(domain(t)))
-    return diagview(DiagonalTensorMap{real(scalartype(t))}(undef, V_cod))
+    T = real(scalartype(t))
+    return SectorVector{T}(undef, V_cod)
 end
 
 # Eigenvalue decomposition
@@ -105,13 +106,13 @@ end
 function MAK.initialize_output(::typeof(eigh_vals!), t::AbstractTensorMap, alg::AbstractAlgorithm)
     V_D = fuse(domain(t))
     T = real(scalartype(t))
-    return diagview(DiagonalTensorMap{Tc}(undef, V_D))
+    return SectorVector{T}(undef, V_D)
 end
 
 function MAK.initialize_output(::typeof(eig_vals!), t::AbstractTensorMap, alg::AbstractAlgorithm)
     V_D = fuse(domain(t))
     Tc = complex(scalartype(t))
-    return diagview(DiagonalTensorMap{Tc}(undef, V_D))
+    return SectorVector{Tc}(undef, V_cod)
 end
 
 # QR decomposition

src/factorizations/truncation.jl

@@ -78,8 +78,7 @@ end
 function MAK.truncate(
         ::typeof(left_null!), (U, S)::NTuple{2, AbstractTensorMap}, strategy::TruncationStrategy
     )
-    extended_S = SectorVector{eltype(S)}(undef, fuse(codomain(U)))
-    fill!(extended_S.data, zero(eltype(extended_S.data)))
+    extended_S = zerovector!(SectorVector{eltype(S)}(undef, fuse(codomain(U))))
     for (c, b) in blocks(S)
         copyto!(extended_S[c], diagview(b)) # copyto! since `b` might be shorter
     end
@@ -92,8 +91,7 @@ end
 function MAK.truncate(
         ::typeof(right_null!), (S, Vᴴ)::NTuple{2, AbstractTensorMap}, strategy::TruncationStrategy
     )
-    extended_S = SectorVector{eltype(S)}(undef, fuse(domain(Vᴴ)))
-    fill!(extended_S.data, zero(eltype(extended_S.data)))
+    extended_S = zerovector!(SectorVector{eltype(S)}(undef, fuse(domain(Vᴴ))))
     for (c, b) in blocks(S)
         copyto!(extended_S[c], diagview(b)) # copyto! since `b` might be shorter
     end
@@ -184,12 +182,12 @@ function _sort_and_perm(values::SectorVector; by = identity, rev::Bool = false)
     perms = SectorDict(
         (
                 begin
-                    p = sortperm(d; by, rev)
+                    p = sortperm(v; by, rev)
                     vs = values_sorted[c]
-                    vs .= d[p]
+                    vs .= view(v, p)
                     c => p
                 end
-            ) for (c, d) in pairs(values)
+            ) for (c, v) in pairs(values)
     )
     return values_sorted, perms
 end
@@ -300,5 +298,5 @@ function MAK.truncation_error!(values::SectorVector, ind)
         v = values[c]
         v[ind_c] .= zero(eltype(v))
     end
-    return TensorKit._norm(pairs(values), 2, zero(real(eltype(valtype(values)))))
+    return norm(values)
 end

src/tensors/sectorvector.jl

@@ -57,3 +57,49 @@ Base.similar(v::SectorVector, V::ElementarySpace) = SectorVector(undef, V)
 blocksectors(v::SectorVector) = keys(v)
 blocks(v::SectorVector) = pairs(v)
 block(v::SectorVector{T, I, A}, c::I) where {T, I, A} = Base.getindex(v, c)
+
+# VectorInterface and LinearAlgebra interface
+# ----------------------------------------------
+VectorInterface.zerovector(v::SectorVector, ::Type{T}) where {T} = 
+    SectorVector(zerovector(parent(v), T), v.structure)
+VectorInterface.zerovector!(v::SectorVector) = (zerovector!(parent(v)); return v)
+VectorInterface.zerovector!!(v::SectorVector) = (zerovector!!(parent(v)); return v)
+
+VectorInterface.scale(v::SectorVector, α) = SectorVector(scale(parent(v), α), v.structure)
+VectorInterface.scale!(v::SectorVector, α) = (scale!(parent(v), α); return v)
+VectorInterface.scale!!(v::SectorVector, α) = (scale!!(parent(v), α); return v)
+
+function VectorInterface.add(v1::SectorVector, v2::SectorVector, α = One(), β = One())
+    SectorVector(add(parent(v1), parent(v2), α, β), v1.structure)
+end
+function VectorInterface.add!(v1::SectorVector, v2::SectorVector, α = One(), β = One())
+    add!(parent(v1), parent(v2), α, β)
+    return v1
+end
+function VectorInterface.add!!(v1::SectorVector, v2::SectorVector, α = One(), β = One())
+    add!!(parent(v1), parent(v2), α, β)
+    return v1
+end
+
+function VectorInterface.inner(v1::SectorVector, v2::SectorVector)
+    v1.structure == v2.structure || throw(SpaceMismatch("Sector structures do not match"))
+    I = sectortype(v1)
+    if FusionStyle(I) isa UniqueFusion # all quantum dimensions are one
+        return inner(parent(v1), parent(v2))
+    else
+        T = VectorInterface.promote_inner(v1, v2)
+        s = zero(T)
+        for c in blocksectors(v1)
+            b1 = block(v1, c)
+            b2 = block(v2, c)
+            s += convert(T, dim(c)) * inner(b1, b2)
+        end
+    end
+    return s
+end
+
+LinearAlgebra.dot(v1::SectorVector, v2::SectorVector) = inner(v1, v2)
+
+function LinearAlgebra.norm(v::SectorVector, p::Real = 2)
+    return _norm(blocks(v), p, float(zero(real(scalartype(v)))))
+end