Implement eig_full!

Author: lkdvos

Project.toml

@@ -6,6 +6,7 @@ version = "0.14.5"
 [deps]
 LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -31,6 +32,7 @@ Combinatorics = "1"
 FiniteDifferences = "0.12"
 LRUCache = "1.0.2"
 LinearAlgebra = "1"
+MatrixAlgebraKit = "0.1.1"
 PackageExtensionCompat = "1"
 Random = "1"
 SparseArrays = "1"

src/TensorKit.jl

@@ -100,6 +100,8 @@ using TensorOperations: TensorOperations, @tensor, @tensoropt, @ncon, ncon
 using TensorOperations: IndexTuple, Index2Tuple, linearize, AbstractBackend
 const TO = TensorOperations
 
+using MatrixAlgebraKit: MatrixAlgebraKit as MAK
+
 using LRUCache
 
 using TensorKitSectors
@@ -194,6 +196,7 @@ include("tensors/treetransformers.jl")
 include("tensors/indexmanipulations.jl")
 include("tensors/diagonal.jl")
 include("tensors/truncation.jl")
+include("tensors/matrixalgebrakit.jl")
 include("tensors/factorizations.jl")
 include("tensors/braidingtensor.jl")
 

src/tensors/matrixalgebrakit.jl

@@ -0,0 +1,49 @@
+function MAK.copy_input(::typeof(MAK.eig_full), t::AbstractTensorMap)
+    return copy_oftype(t, factorisation_scalartype(MAK.eig_full!, t))
+end
+
+function factorisation_scalartype(::typeof(MAK.eig_full!), t::AbstractTensorMap)
+    T = complex(scalartype(t))
+    return promote_type(ComplexF32, typeof(zero(T) / sqrt(abs2(one(T)))))
+end
+
+function MAK.check_input(::typeof(MAK.eig_full!), t::AbstractTensorMap, (D, V))
+    domain(t) == codomain(t) ||
+        throw(ArgumentError("Eigenvalue decomposition requires square input tensor"))
+    Tc = complex(scalartype(t))
+
+    (D isa DiagonalTensorMap &&
+     scalartype(D) == Tc &&
+     fuse(domain(t)) == space(D, 1)) ||
+        throw(ArgumentError("`eig_full!` requires diagonal tensor D with isomorphic domain and complex `scalartype`"))
+
+    V isa AbstractTensorMap &&
+        scalartype(V) == Tc &&
+        space(V) == (codomain(t) ← codomain(D)) ||
+        throw(ArgumentError("`eig_full!` requires square tensor V with isomorphic domain and complex `scalartype`"))
+
+    return nothing
+end
+
+function MAK.initialize_output(::typeof(MAK.eig_full!), t::AbstractTensorMap,
+                               ::MAK.LAPACK_EigAlgorithm)
+    Tc = complex(scalartype(t))
+    V_diag = fuse(domain(t))
+    return DiagonalTensorMap{Tc}(undef, V_diag), similar(t, Tc, domain(t) ← V_diag)
+end
+
+function MAK.eig_full!(t::AbstractTensorMap, (D, V), alg::MAK.LAPACK_EigAlgorithm)
+    MAK.check_input(MAK.eig_full!, t, (D, V))
+    foreachblock(t, D, V) do (_, (b, d, v))
+        d′, v′ = MAK.eig_full!(b, (d, v), alg)
+        # deal with the case where the output is not the same as the input
+        d === d′ || copyto!(d, d′)
+        v === v′ || copyto!(v, v′)
+        return nothing
+    end
+    return D, V
+end
+
+function MAK.default_eig_algorithm(::TensorMap{<:LinearAlgebra.BlasFloat}; kwargs...)
+    return MAK.LAPACK_Expert(; kwargs...)
+end