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,45 @@ 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 also `sort!(::AbstractVector)`.
+
+# Examples
+```jldoctest
+julia> F = 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> F[:values]
+2-element Array{Float64,1}:
+ 4.0
+ 1.0
+
+julia> sort!(F)
+Base.LinAlg.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}([1.0, 4.0], [-0.316228 1.0; 0.948683 0.0])
+
+julia> F[:values]
+2-element Array{Float64,1}:
+ 1.0
+ 4.0
+```
+"""
+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::Eigen; kw...) = sort!(deepcopy(F), kw...)
+
 isposdef(A::Union{Eigen,GeneralizedEigen}) = isreal(A.values) && all(x -> x > 0, A.values)
 
 """

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

test/linalg/eigen.jl

@@ -105,6 +105,12 @@ end
  end
 end
 
+let a = rand(10, 10)
+ f = sort(eigfact(a))
+ @test a ≈ f[:vectors] * Diagonal(f[:values]) / f[:vectors]
+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))