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))