Get eigen working with Array storage

NDTensors/src/NDTensors.jl

@@ -145,6 +145,7 @@ include("arraystorage/arraystorage/tensor/permutedims.jl")
 include("arraystorage/arraystorage/tensor/mul.jl")
 include("arraystorage/arraystorage/tensor/contract.jl")
 include("arraystorage/arraystorage/tensor/qr.jl")
+include("arraystorage/arraystorage/tensor/eigen.jl")
 include("arraystorage/arraystorage/tensor/svd.jl")
 
 # DiagonalArray storage

NDTensors/src/arraystorage/arraystorage/tensor/eigen.jl

@@ -0,0 +1,71 @@
+# TODO: Rewrite this function to be more modern:
+# 1. List keyword arguments in function signature.
+# 2. Output `Spectrum` as a keyword argument that gets overwritten.
+# 3. Dispatch on `alg`.
+# 4. Remove keyword argument deprecations.
+# 5. Make this into two layers, one that handles indices and one that works with `Matrix`.
+# 6. Use `eltype` instead of `where`.
+function eigen(
+  T::Hermitian{ElT,<:ArrayStorageTensor{ElT,2,IndsT}}; kwargs...
+) where {ElT<:Union{Real,Complex},IndsT}
+  # Keyword argument deprecations
+  use_absolute_cutoff = false
+  if haskey(kwargs, :absoluteCutoff)
+    @warn "In svd, keyword argument absoluteCutoff is deprecated in favor of use_absolute_cutoff"
+    use_absolute_cutoff = get(kwargs, :absoluteCutoff, use_absolute_cutoff)
+  end
+  use_relative_cutoff = true
+  if haskey(kwargs, :doRelCutoff)
+    @warn "In svd, keyword argument doRelCutoff is deprecated in favor of use_relative_cutoff"
+    use_relative_cutoff = get(kwargs, :doRelCutoff, use_relative_cutoff)
+  end
+
+  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)
+  maxdim::Int = get(kwargs, :maxdim, minimum(dims(T)))
+  mindim::Int = get(kwargs, :mindim, 1)
+  cutoff::Union{Nothing,Float64} = get(kwargs, :cutoff, 0.0)
+  use_absolute_cutoff::Bool = get(kwargs, :use_absolute_cutoff, use_absolute_cutoff)
+  use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff)
+
+  matrixT = matrix(T)
+  ## TODO Here I am calling parent to ensure that the correct `any` function
+  ## is envoked for non-cpu matrices
+  if any(!isfinite, parent(matrixT))
+    throw(
+      ArgumentError(
+        "Trying to perform the eigendecomposition of a matrix containing NaNs or Infs"
+      ),
+    )
+  end
+
+  DM, VM = eigen(matrixT)
+
+  # Sort by largest to smallest eigenvalues
+  # TODO: Replace `cpu` with `leaf_parenttype` dispatch.
+  p = sortperm(cpu(DM); rev=true, by=abs)
+  DM = DM[p]
+  VM = VM[:, p]
+
+  if truncate
+    DM, truncerr, _ = truncate!!(
+      DM; mindim, maxdim, cutoff, use_absolute_cutoff, use_relative_cutoff, kwargs...
+    )
+    dD = length(DM)
+    if dD < size(VM, 2)
+      VM = VM[:, 1:dD]
+    end
+  else
+    dD = length(DM)
+    truncerr = 0.0
+  end
+  spec = Spectrum(DM, truncerr)
+
+  # Make the new indices to go onto V
+  l = eltype(IndsT)(dD)
+  r = eltype(IndsT)(dD)
+  Vinds = IndsT((dag(ind(T, 2)), dag(r)))
+  Dinds = IndsT((l, dag(r)))
+  V = tensor(VM, Vinds)
+  D = tensor(Diag(DM), Dinds)
+  return D, V, spec
+end

NDTensors/src/arraystorage/arraystorage/tensor/svd.jl

@@ -1,12 +1,12 @@
 # TODO: Rewrite this function to be more modern:
 # 1. List keyword arguments in function signature.
-# 2. Remove `Dense` and `Diag`.
-# 3. Output `Spectrum` as a keyword argument that gets overwritten.
-# 4. Dispatch on `alg`.
-# 5. Remove keyword argument deprecations.
-# 6. Make this into two layers, one that handles indices and one that works with `Matrix`.
+# 2. Output `Spectrum` as a keyword argument that gets overwritten.
+# 3. Dispatch on `alg`.
+# 4. Remove keyword argument deprecations.
+# 5. Make this into two layers, one that handles indices and one that works with `Matrix`.
+# 6. Use `eltype` instead of `where`.
 """
-    svd(T::DenseTensor{<:Number,2}; kwargs...)
+    svd(T::ArrayStorageTensor{<:Number,2}; kwargs...)
 
 svd of an order-2 DenseTensor
 """

NDTensors/src/blocksparse/linearalgebra.jl

@@ -207,7 +207,7 @@ _eigen_eltypes(T::Hermitian{ElT,<:BlockSparseMatrix{ElT}}) where {ElT} = real(El
 _eigen_eltypes(T::BlockSparseMatrix{ElT}) where {ElT} = complex(ElT), complex(ElT)
 
 function eigen(
-  T::Union{Hermitian{ElT,<:BlockSparseMatrix{ElT}},BlockSparseMatrix{ElT}}; kwargs...
+  T::Union{Hermitian{ElT,<:Tensor{ElT,2,<:BlockSparse}},Tensor{ElT,2,<:BlockSparse}}; kwargs...
 ) where {ElT<:Union{Real,Complex}}
   truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)