Implement sorting for Eigen

base/combinatorics.jl

@@ -71,6 +71,35 @@ isperm(p::Tuple{}) = true
 isperm(p::Tuple{Int}) = p[1] == 1
 isperm(p::Tuple{Int,Int}) = ((p[1] == 1) & (p[2] == 2)) | ((p[1] == 2) & (p[2] == 1))
 
+#swap columns i and j of a, in-place
+function swapcols!(a::AbstractMatrix, i, j)
+    i == j && return
+    cols = indices(a,2)
+    (i in cols && j in cols) || throw(BoundsError())
+    for k in indices(a,1)
+        @inbounds a[k,i],a[k,j] = a[k,j],a[k,i]
+    end
+end
+# like permute!! applied to each row of a, in-place in a (overwriting p).
+function permutecols!!(a::AbstractMatrix, p::AbstractVector{<:Integer})
+    count = 0
+    start = 0
+    while count < length(p)
+        ptr = start = findnext(!iszero, p, start+1)
+        next = p[start]
+        count += 1
+        while next != start
+            swapcols!(a, ptr, next)
+            p[ptr] = 0
+            ptr = next
+            next = p[next]
+            count += 1
+        end
+        p[ptr] = 0
+    end
+    a
+end
+
 function permute!!(a, p::AbstractVector{<:Integer})
     count = 0
     start = 0

base/linalg/eigen.jl

@@ -27,6 +27,53 @@ function getindex(A::Union{Eigen,GeneralizedEigen}, d::Symbol)
     throw(KeyError(d))
 end
 
+eigsortby(λ::Real) = λ
+eigsortby(λ::Complex) = reim(λ)
+
+"""
+    sort!(F::Base.LinAlg.Eigen; kw...)
+
+Sort the eigenvectors and eigenvalues of an eigenfactorization `F` in place using the eigenvalues
+for comparison. See [`sort`](@ref) for a variant that returns a sorted copy leaving `F` itself
+unmodified.
+"""
+function sort!(F::Eigen; by=eigsortby, kw...)
+    if !issorted(F[:values], by=by, kw...)
+        perm = sortperm(F[:values]; by=by, kw...)
+        permute!(F[:values], perm)
+        Base.permutecols!!(F[:vectors], perm)
+    end
+    return F
+end
+
+"""
+    sort(F::Base.LinAlg.Eigen; kw...)
+
+Sort the eigenvectors and eigenvalues of an eigenfactorization `F` using the eigenvalues for
+comparison. See [`sort!`](@ref) for a variant that sorts `F` in place. See
+[`sort(::AbstractVector)`](@ref) for a description of possible keyword arguments.
+
+# Examples
+```jldoctest
+julia> F1 = eigfact([4.0 1.0; 0 1.0])
+Base.LinAlg.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}([4.0, 1.0], [1.0 -0.316228; 0.0 0.948683])
+
+julia> F1[:values]
+2-element Array{Float64,1}:
+ 4.0
+ 1.0
+
+julia> F2 = sort(F1)
+Base.LinAlg.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}([1.0, 4.0], [-0.316228 1.0; 0.948683 0.0])
+
+julia> F2[:values]
+2-element Array{Float64,1}:
+ 1.0
+ 4.0
+```
+"""
+sort(F::Eigen; kw...) = sort!(copy(F), kw...)
+
 isposdef(A::Union{Eigen,GeneralizedEigen}) = isreal(A.values) && all(x -> x > 0, A.values)
 
 """
@@ -435,3 +482,5 @@ convert(::Type{AbstractMatrix}, F::Eigen) = F.vectors * Diagonal(F.values) / F.v
 convert(::Type{AbstractArray}, F::Eigen) = convert(AbstractMatrix, F)
 convert(::Type{Matrix}, F::Eigen) = convert(Array, convert(AbstractArray, F))
 convert(::Type{Array}, F::Eigen) = convert(Matrix, F)
+
+copy(F::Eigen) = Eigen(copy(F.values), copy(F.vectors))

base/linalg/linalg.jl

@@ -15,7 +15,7 @@ import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, as
     csch, eltype, exp, eye, findmax, findmin, fill!, floor, getindex, hcat, imag, indices,
     inv, isapprox, isone, IndexStyle, kron, length, log, map, ndims, oneunit, parent,
     power_by_squaring, print_matrix, promote_rule, real, round, sec, sech, setindex!, show, similar,
-    sin, sincos, sinh, size, sqrt, tan, tanh, transpose, trunc, typed_hcat, vec
+    sin, sincos, sinh, size, sort, sort!, sqrt, tan, tanh, transpose, trunc, typed_hcat, vec
 using Base: hvcat_fill, iszero, IndexLinear, _length, promote_op, promote_typeof,
     @propagate_inbounds, @pure, reduce, typed_vcat
 # We use `_length` because of non-1 indices; releases after julia 0.5

doc/src/stdlib/linalg.md

@@ -141,6 +141,8 @@ Base.LinAlg.eigs(::Any, ::Any)
 Base.LinAlg.svds
 Base.LinAlg.peakflops
 Base.LinAlg.stride1
+Base.sort(::Base.LinAlg.Eigen)
+Base.sort!(::Base.LinAlg.Eigen)
 ```
 
 ## Low-level matrix operations

doc/src/stdlib/sort.md

@@ -107,8 +107,8 @@ can be specified via the `lt` keyword.
 ## Sorting Functions
 
 ```@docs
-Base.sort!
-Base.sort
+Base.sort!(::AbstractVector)
+Base.sort(::AbstractVector)
 Base.sortperm
 Base.Sort.sortperm!
 Base.Sort.sortrows

test/linalg/eigen.jl

@@ -105,6 +105,16 @@ end
     end
 end
 
+let aa = rand(10, 10)
+    bb = aa'aa
+    for a in (aa, bb)
+        f = sort(eigfact(a))
+        @test issorted(f[:values], by=Base.LinAlg.eigsortby)
+        @test a ≈ f[:vectors] * Diagonal(f[:values]) / f[:vectors]
+    end
+end
+
+
 # test a matrix larger than 140-by-140 for #14174
 let aa = rand(200, 200)
     for a in (aa, view(aa, 1:n, 1:n))