Updates for Diagonal and rectangular array tests

Author: kshyatt

src/tensors/factorizations/diagonal.jl

@@ -8,17 +8,34 @@ for f in (:svd_compact, :svd_full, :svd_trunc, :svd_vals, :qr_compact, :qr_full,
     @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!)
+for f! in (: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 (:qr_full!, :qr_compact!)
+    @eval function initialize_output(::typeof($f!), d::AbstractTensorMap,
+                                     ::DiagonalAlgorithm)
+        return d, similar(d)
+    end
+    # to avoid ambiguities
+    @eval function initialize_output(::typeof($f!), d::AdjointTensorMap,
+                                     ::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
+    # to avoid ambiguities
+    @eval function initialize_output(::typeof($f!), d::AdjointTensorMap,
+                                     ::DiagonalAlgorithm)
+        return similar(d), d
+    end
 end
 
 for f! in

src/tensors/factorizations/factorizations.jl

@@ -61,10 +61,10 @@ const RealOrComplexFloat = Union{AbstractFloat,Complex{<:AbstractFloat}}
 #------------------------------#
 # LinearAlgebra overloads
 #------------------------------#
-LinearAlgebra.svdvals(t::AbstractTensorMap) = diagview(svd_vals(t))
+#LinearAlgebra.svdvals(t::AbstractTensorMap)  = diagview(svd_vals(t))
 LinearAlgebra.svdvals!(t::AbstractTensorMap) = diagview(svd_vals!(t))
-LinearAlgebra.eigvals(t::AbstractTensorMap) = diagview(eigvals(t))
-LinearAlgebra.eigvals!(t::AbstractTensorMap) = diagview(eigvals!(t))
+#LinearAlgebra.eigvals(t::AbstractTensorMap)  = diagview(eig_vals(t))
+LinearAlgebra.eigvals!(t::AbstractTensorMap; kwargs...) = diagview(eig_vals!(t))
 
 #--------------------------------------------------#
 # Checks for hermiticity and positive definiteness #

src/tensors/factorizations/matrixalgebrakit.jl

@@ -45,7 +45,7 @@ for f! in (:qr_compact!, :qr_full!,
 end
 
 # Handle these separately because single output instead of tuple
-for f! in (:qr_null!, :lq_null!, :svd_vals!, :eig_vals!, :eigh_vals!)
+for f! in (:qr_null!, :lq_null!)
     @eval function $f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm)
         check_input($f!, t, N, alg)
 
@@ -60,6 +60,22 @@ for f! in (:qr_null!, :lq_null!, :svd_vals!, :eig_vals!, :eigh_vals!)
     end
 end
 
+# Handle these separately because single output instead of tuple
+for f! in (:svd_vals!, :eig_vals!, :eigh_vals!)
+    @eval function $f!(t::AbstractTensorMap, N, alg::AbstractAlgorithm)
+        check_input($f!, t, N, alg)
+
+        foreachblock(t, N) do _, (b, n)
+            n′ = $f!(b, n.diag, alg)
+            # deal with the case where the output is not the same as the input
+            n.diag === n′ || copyto!(n, diagview(n′))
+            return nothing
+        end
+
+        return N
+    end
+end
+
 # Singular value decomposition
 # ----------------------------
 const _T_USVᴴ = Tuple{<:AbstractTensorMap,<:AbstractTensorMap,<:AbstractTensorMap}
@@ -101,11 +117,19 @@ end
 function check_input(::typeof(svd_vals!), t::AbstractTensorMap, S::SectorDict,
                      ::AbstractAlgorithm)
     @check_scalar S t real
-    V_cod = infimum(fuse(codomain(t)), fuse(domain(t)))
+    V_cod = V_dom = infimum(fuse(codomain(t)), fuse(domain(t)))
     @check_space(S, V_cod ← V_dom)
     return nothing
 end
 
+function check_input(::typeof(svd_vals!), t::AbstractTensorMap, D::DiagonalTensorMap,
+                     ::AbstractAlgorithm)
+    @check_scalar D t real
+    V_cod = V_dom = infimum(fuse(codomain(t)), fuse(domain(t)))
+    @check_space(D, V_cod ← V_dom)
+    return nothing
+end
+
 function initialize_output(::typeof(svd_full!), t::AbstractTensorMap, ::AbstractAlgorithm)
     V_cod = fuse(codomain(t))
     V_dom = fuse(domain(t))

test/factorizations.jl

@@ -27,9 +27,8 @@ for V in spacelist
             W = V1 ⊗ V2
             @testset for T in (Float32, ComplexF64)
                 # Test both a normal tensor and an adjoint one.
-                ts = (rand(T, W, W'), rand(T, W, W')')
+                ts = (rand(T, W, W'), rand(T, W, W')', rand(T, V1, W'), rand(T, V1, W')')
                 @testset for t in ts
-                    # test squares and rectangles here
                     @testset "leftorth with $alg" for alg in
                                                       (TensorKit.LAPACK_HouseholderQR(),
                                                        TensorKit.LAPACK_HouseholderQR(;
@@ -40,10 +39,10 @@ for V in spacelist
                                                        TensorKit.PolarViaSVD(TensorKit.LAPACK_DivideAndConquer()),
                                                        TensorKit.LAPACK_QRIteration(),
                                                        TensorKit.LAPACK_DivideAndConquer())
+                        (codomain(t) ≾ domain(t)) && alg isa TensorKit.PolarViaSVD && continue
                         Q, R = @constinferred leftorth(t; alg=alg)
                         @test isisometry(Q)
-                        tQR = Q * R
-                        @test tQR ≈ t
+                        @test Q * R ≈ t
                     end
                     @testset "leftnull with $alg" for alg in
                                                       (TensorKit.LAPACK_HouseholderQR(),
@@ -61,6 +60,7 @@ for V in spacelist
                                                         TensorKit.PolarViaSVD(TensorKit.LAPACK_DivideAndConquer()),
                                                         TensorKit.LAPACK_QRIteration(),
                                                         TensorKit.LAPACK_DivideAndConquer())
+                        (domain(t) ≾ codomain(t)) && alg isa TensorKit.PolarViaSVD && continue
                         L, Q = @constinferred rightorth(t; alg=alg)
                         @test isisometry(Q; side=:right)
                         @test L * Q ≈ t
@@ -80,28 +80,28 @@ for V in spacelist
                         @test isisometry(V; side=:right)
                         @test U * S * V ≈ t
 
-                        s = LinearAlgebra.svdvals(t)
+                        s  = LinearAlgebra.svdvals(t)
                         s′ = LinearAlgebra.diag(S)
                         for (c, b) in s
                             @test b ≈ s′[c]
                         end
-                        s = LinearAlgebra.svdvals(t')
+                        s  = LinearAlgebra.svdvals(t')
                         s′ = LinearAlgebra.diag(S')
                         for (c, b) in s
                             @test b ≈ s′[c]
                         end
                     end
                     @testset "cond and rank" begin
-                        d1 = dim(codomain(t))
-                        d2 = dim(domain(t))
+                        d1   = dim(codomain(t))
+                        d2   = dim(domain(t))
                         @test rank(t) == min(d1, d2)
-                        M = leftnull(t)
-                        @test rank(M) == max(d1, d2) - min(d1, d2)
-                        t3 = unitary(T, V1 ⊗ V2, V1 ⊗ V2)
+                        M    = leftnull(t)
+                        @test rank(M) + rank(t) == d1
+                        t3   = unitary(T, V1 ⊗ V2, V1 ⊗ V2)
                         @test cond(t3) ≈ one(real(T))
                         @test rank(t3) == dim(V1 ⊗ V2)
-                        t4 = randn(T, V1 ⊗ V2, V1 ⊗ V2)
-                        t4 = (t4 + t4') / 2
+                        t4   = randn(T, V1 ⊗ V2, V1 ⊗ V2)
+                        t4   = (t4 + t4') / 2
                         vals = LinearAlgebra.eigvals(t4)
                         λmax = maximum(s -> maximum(abs, s), values(vals))
                         λmin = minimum(s -> minimum(abs, s), values(vals))