Generalize 2-arg eigen towards AbstractMatrix (#46797)

Author: dkarrasch

stdlib/LinearAlgebra/src/bunchkaufman.jl

@@ -197,7 +197,7 @@ julia> S.L*S.D*S.L' - A[S.p, S.p]
 ```
 """
 bunchkaufman(A::AbstractMatrix{T}, rook::Bool=false; check::Bool = true) where {T} =
-    bunchkaufman!(copymutable_oftype(A, typeof(sqrt(oneunit(T)))), rook; check = check)
+    bunchkaufman!(eigencopy_oftype(A, typeof(sqrt(oneunit(T)))), rook; check = check)
 
 BunchKaufman{T}(B::BunchKaufman) where {T} =
     BunchKaufman(convert(Matrix{T}, B.LD), B.ipiv, B.uplo, B.symmetric, B.rook, B.info)

stdlib/LinearAlgebra/src/diagonal.jl

@@ -752,7 +752,7 @@ function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true, sortby::Union{
     λ = eigvals(D)
     if !isnothing(sortby)
         p = sortperm(λ; alg=QuickSort, by=sortby)
-        λ = λ[p] # make a copy, otherwise this permutes D.diag
+        λ = λ[p]
         evecs = zeros(Td, size(D))
         @inbounds for i in eachindex(p)
             evecs[p[i],i] = one(Td)
@@ -762,6 +762,29 @@ function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true, sortby::Union{
     end
     Eigen(λ, evecs)
 end
+function eigen(Da::Diagonal, Db::Diagonal; sortby::Union{Function,Nothing}=nothing)
+    if any(!isfinite, Da.diag) || any(!isfinite, Db.diag)
+        throw(ArgumentError("matrices contain Infs or NaNs"))
+    end
+    if any(iszero, Db.diag)
+        throw(ArgumentError("right-hand side diagonal matrix is singular"))
+    end
+    return GeneralizedEigen(eigen(Db \ Da; sortby)...)
+end
+function eigen(A::AbstractMatrix, D::Diagonal; sortby::Union{Function,Nothing}=nothing)
+    if any(iszero, D.diag)
+        throw(ArgumentError("right-hand side diagonal matrix is singular"))
+    end
+    if size(A, 1) == size(A, 2) && isdiag(A)
+        return eigen(Diagonal(A), D; sortby)
+    elseif ishermitian(A)
+        S = promote_type(eigtype(eltype(A)), eltype(D))
+        return eigen!(eigencopy_oftype(Hermitian(A), S), Diagonal{S}(D); sortby)
+    else
+        S = promote_type(eigtype(eltype(A)), eltype(D))
+        return eigen!(eigencopy_oftype(A, S), Diagonal{S}(D); sortby)
+    end
+end
 
 #Singular system
 svdvals(D::Diagonal{<:Number}) = sort!(abs.(D.diag), rev = true)

stdlib/LinearAlgebra/src/eigen.jl

@@ -233,16 +233,20 @@ true
 ```
 """
 function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T
-    AA = copymutable_oftype(A, eigtype(T))
-    isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
-    return eigen!(AA; permute=permute, scale=scale, sortby=sortby)
+    isdiag(A) && return eigen(Diagonal{eigtype(T)}(diag(A)); sortby)
+    ishermitian(A) && return eigen!(eigencopy_oftype(Hermitian(A), eigtype(T)); sortby)
+    AA = eigencopy_oftype(A, eigtype(T))
+    return eigen!(AA; permute, scale, sortby)
 end
 function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where {T <: Union{Float16,Complex{Float16}}}
-    AA = copymutable_oftype(A, eigtype(T))
-    isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
-    A = eigen!(AA; permute, scale, sortby)
-    values = convert(AbstractVector{isreal(A.values) ? Float16 : Complex{Float16}}, A.values)
-    vectors = convert(AbstractMatrix{isreal(A.vectors) ? Float16 : Complex{Float16}}, A.vectors)
+    isdiag(A) && return eigen(Diagonal{eigtype(T)}(diag(A)); sortby)
+    E = if ishermitian(A)
+        eigen!(eigencopy_oftype(Hermitian(A), eigtype(T)); sortby)
+    else
+        eigen!(eigencopy_oftype(A, eigtype(T)); permute, scale, sortby)
+    end
+    values = convert(AbstractVector{isreal(E.values) ? Float16 : Complex{Float16}}, E.values)
+    vectors = convert(AbstractMatrix{isreal(E.vectors) ? Float16 : Complex{Float16}}, E.vectors)
     return Eigen(values, vectors)
 end
 eigen(x::Number) = Eigen([x], fill(one(x), 1, 1))
@@ -333,7 +337,7 @@ julia> eigvals(diag_matrix)
 ```
 """
 eigvals(A::AbstractMatrix{T}; kws...) where T =
-    eigvals!(copymutable_oftype(A, eigtype(T)); kws...)
+    eigvals!(eigencopy_oftype(A, eigtype(T)); kws...)
 
 """
 For a scalar input, `eigvals` will return a scalar.
@@ -507,12 +511,20 @@ true
 ```
 """
 function eigen(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
-    S = promote_type(eigtype(TA),TB)
-    eigen!(copymutable_oftype(A, S), copymutable_oftype(B, S); kws...)
+    S = promote_type(eigtype(TA), TB)
+    eigen!(eigencopy_oftype(A, S), eigencopy_oftype(B, S); kws...)
 end
-
 eigen(A::Number, B::Number) = eigen(fill(A,1,1), fill(B,1,1))
 
+"""
+    LinearAlgebra.eigencopy_oftype(A::AbstractMatrix, ::Type{S})
+
+Creates a dense copy of `A` with eltype `S` by calling `copy_similar(A, S)`.
+In the case of `Hermitian` or `Symmetric` matrices additionally retains the wrapper,
+together with the `uplo` field.
+"""
+eigencopy_oftype(A, S) = copy_similar(A, S)
+
 """
     eigvals!(A, B; sortby) -> values
 
@@ -586,8 +598,8 @@ julia> eigvals(A,B)
 ```
 """
 function eigvals(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
-    S = promote_type(eigtype(TA),TB)
-    return eigvals!(copymutable_oftype(A, S), copymutable_oftype(B, S); kws...)
+    S = promote_type(eigtype(TA), TB)
+    return eigvals!(eigencopy_oftype(A, S), eigencopy_oftype(B, S); kws...)
 end
 
 """

stdlib/LinearAlgebra/src/symmetriceigen.jl

@@ -1,13 +1,16 @@
 # This file is a part of Julia. License is MIT: https://julialang.org/license
 
+# preserve HermOrSym wrapper
+eigencopy_oftype(A::Hermitian, S) = Hermitian(copy_similar(A, S), sym_uplo(A.uplo))
+eigencopy_oftype(A::Symmetric, S) = Symmetric(copy_similar(A, S), sym_uplo(A.uplo))
+
 # Eigensolvers for symmetric and Hermitian matrices
 eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing) =
     Eigen(sorteig!(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)..., sortby)...)
 
 function eigen(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing)
-    T = eltype(A)
-    S = eigtype(T)
-    eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), sortby=sortby)
+    S = eigtype(eltype(A))
+    eigen!(eigencopy_oftype(A, S), sortby=sortby)
 end
 
 eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, irange::UnitRange) =
@@ -31,9 +34,8 @@ The [`UnitRange`](@ref) `irange` specifies indices of the sorted eigenvalues to
     will be a *truncated* factorization.
 """
 function eigen(A::RealHermSymComplexHerm, irange::UnitRange)
-    T = eltype(A)
-    S = eigtype(T)
-    eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange)
+    S = eigtype(eltype(A))
+    eigen!(eigencopy_oftype(A, S), irange)
 end
 
 eigen!(A::RealHermSymComplexHerm{T,<:StridedMatrix}, vl::Real, vh::Real) where {T<:BlasReal} =
@@ -57,9 +59,8 @@ The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`
     will be a *truncated* factorization.
 """
 function eigen(A::RealHermSymComplexHerm, vl::Real, vh::Real)
-    T = eltype(A)
-    S = eigtype(T)
-    eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh)
+    S = eigtype(eltype(A))
+    eigen!(eigencopy_oftype(A, S), vl, vh)
 end
 
 function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby::Union{Function,Nothing}=nothing)
@@ -69,9 +70,8 @@ function eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby:
 end
 
 function eigvals(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing)
-    T = eltype(A)
-    S = eigtype(T)
-    eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), sortby=sortby)
+    S = eigtype(eltype(A))
+    eigvals!(eigencopy_oftype(A, S), sortby=sortby)
 end
 
 """
@@ -110,9 +110,8 @@ julia> eigvals(A)
 ```
 """
 function eigvals(A::RealHermSymComplexHerm, irange::UnitRange)
-    T = eltype(A)
-    S = eigtype(T)
-    eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange)
+    S = eigtype(eltype(A))
+    eigvals!(eigencopy_oftype(A, S), irange)
 end
 
 """
@@ -150,9 +149,8 @@ julia> eigvals(A)
 ```
 """
 function eigvals(A::RealHermSymComplexHerm, vl::Real, vh::Real)
-    T = eltype(A)
-    S = eigtype(T)
-    eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh)
+    S = eigtype(eltype(A))
+    eigvals!(eigencopy_oftype(A, S), vl, vh)
 end
 
 eigmax(A::RealHermSymComplexHerm{<:Real,<:StridedMatrix}) = eigvals(A, size(A, 1):size(A, 1))[1]
@@ -166,107 +164,25 @@ function eigen!(A::Hermitian{T,S}, B::Hermitian{T,S}; sortby::Union{Function,Not
     vals, vecs, _ = LAPACK.sygvd!(1, 'V', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))
     GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
 end
-
 function eigen!(A::RealHermSymComplexHerm{T,S}, B::AbstractMatrix{T}; sortby::Union{Function,Nothing}=nothing) where {T<:Number,S<:StridedMatrix}
+    return _choleigen!(A, B, sortby)
+end
+function eigen!(A::StridedMatrix{T}, B::Union{RealHermSymComplexHerm{T},Diagonal{T}}; sortby::Union{Function,Nothing}=nothing) where {T<:Number}
+    return _choleigen!(A, B, sortby)
+end
+function _choleigen!(A, B, sortby)
     U = cholesky(B).U
     vals, w = eigen!(UtiAUi!(A, U))
     vecs = U \ w
     GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
 end
 
-# Perform U' \ A / U in-place.
-UtiAUi!(As::Symmetric, Utr::UpperTriangular) = Symmetric(_UtiAsymUi!(As.uplo, parent(As), parent(Utr)), sym_uplo(As.uplo))
-UtiAUi!(As::Hermitian, Utr::UpperTriangular) = Hermitian(_UtiAsymUi!(As.uplo, parent(As), parent(Utr)), sym_uplo(As.uplo))
-UtiAUi!(As::Symmetric, Udi::Diagonal) = Symmetric(_UtiAsymUi_diag!(As.uplo, parent(As), Udi), sym_uplo(As.uplo))
-UtiAUi!(As::Hermitian, Udi::Diagonal) = Hermitian(_UtiAsymUi_diag!(As.uplo, parent(As), Udi), sym_uplo(As.uplo))
-
-# U is upper triangular
-function _UtiAsymUi!(uplo, A, U)
-    n = size(A, 1)
-    μ⁻¹ = 1 / U[1, 1]
-    αμ⁻² = A[1, 1] * μ⁻¹' * μ⁻¹
-
-    # Update (1, 1) element
-    A[1, 1] = αμ⁻²
-    if n > 1
-        Unext = view(U, 2:n, 2:n)
-
-        if uplo === 'U'
-            # Update submatrix
-            for j in 2:n, i in 2:j
-                A[i, j] = (
-                    A[i, j]
-                    - μ⁻¹' * U[1, j] * A[1, i]'
-                    - μ⁻¹ * A[1, j] * U[1, i]'
-                    + αμ⁻² * U[1, j] * U[1, i]'
-                )
-            end
-
-            # Update vector
-            for j in 2:n
-                A[1, j] = A[1, j] * μ⁻¹' - U[1, j] * αμ⁻²
-            end
-            ldiv!(view(A', 2:n, 1), UpperTriangular(Unext)', view(A', 2:n, 1))
-        else
-            # Update submatrix
-            for j in 2:n, i in 2:j
-                A[j, i] = (
-                    A[j, i]
-                    - μ⁻¹ * A[i, 1]' * U[1, j]'
-                    - μ⁻¹' * U[1, i] * A[j, 1]
-                    + αμ⁻² * U[1, i] * U[1, j]'
-                )
-            end
-
-            # Update vector
-            for j in 2:n
-                A[j, 1] = A[j, 1] * μ⁻¹ - U[1, j]' * αμ⁻²
-            end
-            ldiv!(view(A, 2:n, 1), UpperTriangular(Unext)', view(A, 2:n, 1))
-        end
-
-        # Recurse
-        _UtiAsymUi!(uplo, view(A, 2:n, 2:n), Unext)
-    end
-
-    return A
-end
+# Perform U' \ A / U in-place, where U::Union{UpperTriangular,Diagonal}
+UtiAUi!(A::StridedMatrix, U) = _UtiAUi!(A, U)
+UtiAUi!(A::Symmetric, U) = Symmetric(_UtiAUi!(copytri!(parent(A), A.uplo), U), sym_uplo(A.uplo))
+UtiAUi!(A::Hermitian, U) = Hermitian(_UtiAUi!(copytri!(parent(A), A.uplo, true), U), sym_uplo(A.uplo))
 
-# U is diagonal
-function _UtiAsymUi_diag!(uplo, A, U)
-    n = size(A, 1)
-    μ⁻¹ = 1 / U[1, 1]
-    αμ⁻² = A[1, 1] * μ⁻¹' * μ⁻¹
-
-    # Update (1, 1) element
-    A[1, 1] = αμ⁻²
-    if n > 1
-        Unext = view(U, 2:n, 2:n)
-
-        if uplo === 'U'
-            # No need to update any submatrix when U is diagonal
-
-            # Update vector
-            for j in 2:n
-                A[1, j] = A[1, j] * μ⁻¹'
-            end
-            ldiv!(view(A', 2:n, 1), Diagonal(Unext)', view(A', 2:n, 1))
-        else
-            # No need to update any submatrix when U is diagonal
-
-            # Update vector
-            for j in 2:n
-                A[j, 1] = A[j, 1] * μ⁻¹
-            end
-            ldiv!(view(A, 2:n, 1), Diagonal(Unext)', view(A, 2:n, 1))
-        end
-
-        # Recurse
-        _UtiAsymUi!(uplo, view(A, 2:n, 2:n), Unext)
-    end
-
-    return A
-end
+_UtiAUi!(A, U) = rdiv!(ldiv!(U', A), U)
 
 function eigvals!(A::HermOrSym{T,S}, B::HermOrSym{T,S}; sortby::Union{Function,Nothing}=nothing) where {T<:BlasReal,S<:StridedMatrix}
     vals = LAPACK.sygvd!(1, 'N', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))[1]