Comments

Author: kshyatt

ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl

@@ -49,27 +49,11 @@ MatrixAlgebraKit.initialize_output(::typeof(eigh_vals!), A::AbstractMatrix, ::GL
 
 function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix, DV, ::GLA_QRIteration)
     eigval, eigvec = eigen!(Hermitian(A); sortby = real)
-    D, V = DV
-    if isnothing(D)
-        D = Diagonal(eigval::AbstractVector{real(eltype(A))})
-    else
-        copyto!(D, Diagonal(eigval::AbstractVector{real(eltype(A))}))
-    end
-    if isnothing(V)
-        V = eigvec::AbstractMatrix{eltype(A)}
-    else
-        copyto!(V, eigvec::AbstractMatrix{eltype(A)})
-    end
-    return D, V
+    return Diagonal(eigval::AbstractVector{real(eltype(A))}), eigvec::AbstractMatrix{eltype(A)}
 end
 
 function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix, D, ::GLA_QRIteration)
-    if isnothing(D)
-        D = eigvals!(Hermitian(A); sortby = real)
-    else
-        copyto!(D, eigvals!(Hermitian(A); sortby = real))
-    end
-    return D
+    return eigvals!(Hermitian(A); sortby = real)
 end
 
 function MatrixAlgebraKit.default_qr_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{Float16, ComplexF16, BigFloat, Complex{BigFloat}}}}

test/testsuite/eigh.jl

@@ -26,18 +26,16 @@ function test_eigh_full(
     )
     summary_str = testargs_summary(T, sz)
     return @testset "eigh_full! $summary_str" begin
-        A = instantiate_matrix(T, sz)
-        A = (A + A') / 2
+        A = project_hermitian!(instantiate_matrix(T, sz))
         Ac = deepcopy(A)
 
         D, V = @testinferred eigh_full(A)
         @test A * V ≈ V * D
         @test isunitary(V)
         @test all(isreal, D)
 
-        D2, V2 = eigh_full!(copy(A), (D, V))
-        @test D2 === D
-        @test V2 === V
+        D2, V2 = eigh_full!(Ac, (D, V))
+        @test A * V2 ≈ V2 * D2
 
         D3 = @testinferred eigh_vals(A)
         @test D ≈ Diagonal(D3)
@@ -51,18 +49,16 @@ function test_eigh_full_algs(
     )
     summary_str = testargs_summary(T, sz)
     return @testset "eigh_full! algorithm $alg $summary_str" for alg in algs
-        A = instantiate_matrix(T, sz)
-        A = (A + A') / 2
+        A = project_hermitian!(instantiate_matrix(T, sz))
         Ac = deepcopy(A)
 
         D, V = @testinferred eigh_full(A; alg)
         @test A * V ≈ V * D
         @test isunitary(V)
         @test all(isreal, D)
 
-        D2, V2 = eigh_full!(copy(A), (D, V); alg)
-        @test D2 === D
-        @test V2 === V
+        D2, V2 = eigh_full!(Ac, (D, V); alg)
+        @test A * V2 ≈ V2 * D2
 
         D3 = @testinferred eigh_vals(A; alg)
         @test D ≈ Diagonal(D3)
@@ -76,9 +72,7 @@ function test_eigh_trunc(
     )
     summary_str = testargs_summary(T, sz)
     return @testset "eigh_trunc! $summary_str" begin
-        A = instantiate_matrix(T, sz)
-        A = A * A'
-        A = (A + A') / 2
+        A = project_hermitian!(instantiate_matrix(T, sz))
         Ac = deepcopy(A)
         if !(T <: Diagonal)
 
@@ -132,8 +126,7 @@ function test_eigh_trunc(
             Ddiag = similar(A, real(eltype(T)), m4)
             copyto!(Ddiag, real(eltype(T))[0.9, 0.3, 0.1, 0.01])
             D = Diagonal(Ddiag)
-            A = V * D * V'
-            A = (A + A') / 2
+            A = project_hermitian!(V * D * V')
             alg = TruncatedAlgorithm(MatrixAlgebraKit.default_eigh_algorithm(A), truncrank(2))
             D2, V2, ϵ2 = @testinferred eigh_trunc(A; alg)
             @test diagview(D2) ≈ diagview(D)[1:2]
@@ -155,9 +148,7 @@ function test_eigh_trunc_algs(
     )
     summary_str = testargs_summary(T, sz)
     return @testset "eigh_trunc! algorithm $alg $summary_str" for alg in algs
-        A = instantiate_matrix(T, sz)
-        A = A * A'
-        A = (A + A') / 2
+        A = project_hermitian!(instantiate_matrix(T, sz))
         Ac = deepcopy(A)
 
         m = size(A, 1)
@@ -172,8 +163,7 @@ function test_eigh_trunc_algs(
         Ddiag = similar(A, real(eltype(T)), m4)
         copyto!(Ddiag, real(eltype(T))[0.9, 0.3, 0.1, 0.01])
         D = Diagonal(Ddiag)
-        A = V * D * V'
-        A = (A + A') / 2
+        A = project_hermitian!(V * D * V')
         truncalg = TruncatedAlgorithm(alg, truncrank(2))
         D2, V2, ϵ2 = @testinferred eigh_trunc(A; alg = truncalg)
         @test diagview(D2) ≈ diagview(D)[1:2]