clean up handling DiagonalTensorMap

Author: lkdvos

Project.toml

@@ -33,7 +33,7 @@ Combinatorics = "1"
 FiniteDifferences = "0.12"
 LRUCache = "1.0.2"
 LinearAlgebra = "1"
-MatrixAlgebraKit = "0.3"
+MatrixAlgebraKit = "0.3.1"
 OhMyThreads = "0.8.0"
 PackageExtensionCompat = "1"
 Random = "1"

src/tensors/factorizations/diagonal.jl

@@ -0,0 +1,78 @@
+# DiagonalTensorMap
+# -----------------
+_repack_diagonal(d::DiagonalTensorMap) = Diagonal(d.data)
+
+for f in [
+        :svd_compact, :svd_full, :svd_trunc, :svd_vals, :qr_compact, :qr_full, :qr_null,
+        :lq_compact, :lq_full, :lq_null, :eig_full, :eig_trunc, :eig_vals, :eigh_full,
+        :eigh_trunc, :eigh_vals, :left_polar, :right_polar,
+    ]
+    @eval copy_input(::typeof($f), d::DiagonalTensorMap) = copy(d)
+end
+
+for f! in (:qr_full!, :qr_compact!, :eig_full!, :eig_trunc!, :eigh_full!, :eigh_trunc!)
+    @eval function initialize_output(::typeof($f!), d::AbstractTensorMap, ::DiagonalAlgorithm)
+        return d, similar(d)
+    end
+end
+for f! in (:lq_full!, :lq_compact!)
+    @eval function initialize_output(::typeof($f!), d::AbstractTensorMap, ::DiagonalAlgorithm)
+        return similar(d), d
+    end
+end
+
+for f! in (:qr_full!, :qr_compact!, :lq_full!, :lq_compact!, :eig_full!, :eig_trunc!, :eigh_full!, :eigh_trunc!)
+    @eval function $f!(d::DiagonalTensorMap, F, alg::DiagonalAlgorithm)
+        check_input($f!, d, F, alg)
+        $f!(_repack_diagonal(d), _repack_diagonal.(F), alg)
+        return F
+    end
+end
+
+for f! in (:qr_full!, :qr_compact!)
+    @eval function check_input(::typeof($f!), d::AbstractTensorMap, (Q, R)::_T_QR, ::DiagonalAlgorithm)
+        @assert d isa DiagonalTensorMap
+        @assert Q isa DiagonalTensorMap && R isa DiagonalTensorMap
+        @check_scalar Q d
+        @check_scalar R d
+        @check_space(Q, space(d))
+        @check_space(R, space(d))
+
+        return nothing
+    end
+end
+
+for f! in (:lq_full!, :lq_compact!)
+    @eval function check_input(::typeof($f!), d::AbstractTensorMap, (L, Q)::_T_LQ, ::DiagonalAlgorithm)
+        @assert d isa DiagonalTensorMap
+        @assert Q isa DiagonalTensorMap && L isa DiagonalTensorMap
+        @check_scalar Q d
+        @check_scalar L d
+        @check_space(Q, space(d))
+        @check_space(L, space(d))
+
+        return nothing
+    end
+end
+
+# f_vals
+# ------
+
+for f! in (:eig_vals!, :eigh_vals!, :svd_vals!)
+    @eval function $f!(d::AbstractTensorMap, V, alg::DiagonalAlgorithm)
+        check_input($f!, d, V, alg)
+        $f!(_repack_diagonal(d), diagview(_repack_diagonal(V)), alg)
+        return V
+    end
+    @eval function initialize_output(::typeof($f!), d::DiagonalTensorMap, alg::DiagonalAlgorithm)
+        data = initialize_output($f!, _repack_diagonal(d), alg)
+        return DiagonalTensorMap(data, d.domain)
+    end
+end
+
+function check_input(::typeof(eig_vals!), t::AbstractTensorMap, D::DiagonalTensorMap,
+                     ::DiagonalAlgorithm)
+    @check_scalar D t
+    @check_space D space(t)
+    return nothing
+end

src/tensors/factorizations/factorizations.jl

@@ -21,7 +21,7 @@ using TensorOperations: Index2Tuple
 using MatrixAlgebraKit
 using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, TruncationStrategy,
                         NoTruncation, TruncationKeepAbove, TruncationKeepBelow,
-                        TruncationIntersection, TruncationKeepFiltered
+                        TruncationIntersection, TruncationKeepFiltered, DiagonalAlgorithm
 import MatrixAlgebraKit: default_algorithm,
                          copy_input, check_input, initialize_output,
                          qr_compact!, qr_full!, qr_null!, lq_compact!, lq_full!, lq_null!,
@@ -41,6 +41,7 @@ include("matrixalgebrakit.jl")
 include("truncation.jl")
 include("deprecations.jl")
 include("adjoint.jl")
+include("diagonal.jl")
 
 TensorKit.one!(A::AbstractMatrix) = MatrixAlgebraKit.one!(A)
 
@@ -55,59 +56,6 @@ end
 #------------------------------------------------------------------------------------------
 const RealOrComplexFloat = Union{AbstractFloat,Complex{<:AbstractFloat}}
 
-# DiagonalTensorMap
-# -----------------
-function leftorth!(d::DiagonalTensorMap; alg=QR(), kwargs...)
-    @assert alg isa Union{QR,QL}
-    return one(d), d # TODO: this is only correct for `alg = QR()` or `alg = QL()`
-end
-function rightorth!(d::DiagonalTensorMap; alg=LQ(), kwargs...)
-    @assert alg isa Union{LQ,RQ}
-    return d, one(d) # TODO: this is only correct for `alg = LQ()` or `alg = RQ()`
-end
-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] = 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
-
-eig!(d::DiagonalTensorMap) = d, one(d)
-eigh!(d::DiagonalTensorMap{<:Real}) = d, one(d)
-eigh!(d::DiagonalTensorMap{<:Complex}) = DiagonalTensorMap(real(d.data), d.domain), one(d)
-
-function LinearAlgebra.svdvals(d::DiagonalTensorMap)
-    return SectorDict(c => LinearAlgebra.svdvals(b) for (c, b) in blocks(d))
-end
-function LinearAlgebra.eigvals(d::DiagonalTensorMap)
-    return SectorDict(c => LinearAlgebra.eigvals(b) for (c, b) in blocks(d))
-end
-
-function LinearAlgebra.cond(d::DiagonalTensorMap, p::Real=2)
-    return LinearAlgebra.cond(Diagonal(d.data), p)
-end
 #------------------------------#
 # Singular value decomposition #
 #------------------------------#