diff --git a/NEWS.md b/NEWS.md
index 72681e1ba2b7b..976628e69a0be 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -704,9 +704,9 @@ Deprecated or removed
     `type` is fully deprecated to `mutable struct` ([#19157], [#20418]).
 
   * `lufact`, `schurfact`, `lqfact`, `qrfact`, `ldltfact`, `svdfact`,
-    `bkfact`, and `hessfact` have respectively been deprecated to
+    `bkfact`, `hessfact`, and `eigfact` have respectively been deprecated to
     `lu`, `schur`, `lq`, `qr`, `ldlt`, `svd`, `bunchkaufman`,
-    and `hessenberg` ([#27159]).
+    `hessenberg`, and `eigen` ([#27159]).
 
   * `lufact!`, `schurfact!`, `lqfact!`, `qrfact!`, `ldltfact!`, `svdfact!`,
     `bkfact!`, and `hessfact!` have respectively been deprecated to
diff --git a/stdlib/IterativeEigensolvers/src/IterativeEigensolvers.jl b/stdlib/IterativeEigensolvers/src/IterativeEigensolvers.jl
index 9b527c6112453..ca2f543a8dd8e 100644
--- a/stdlib/IterativeEigensolvers/src/IterativeEigensolvers.jl
+++ b/stdlib/IterativeEigensolvers/src/IterativeEigensolvers.jl
@@ -78,7 +78,7 @@ function _eigs(A, B;
     sym = !iscmplx && issymmetric(A) && issymmetric(B)
     nevmax = sym ? n-1 : n-2
     if nevmax <= 0
-        throw(ArgumentError("input matrix A is too small. Use eigfact instead."))
+        throw(ArgumentError("input matrix A is too small. Use eigen instead."))
     end
     if nev > nevmax
         @warn "Adjusting nev from $nev to $nevmax"
diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md
index 9d55519cdd55d..5b081e7aeefdf 100644
--- a/stdlib/LinearAlgebra/docs/src/index.md
+++ b/stdlib/LinearAlgebra/docs/src/index.md
@@ -337,7 +337,7 @@ LinearAlgebra.eigvals!
 LinearAlgebra.eigmax
 LinearAlgebra.eigmin
 LinearAlgebra.eigvecs
-LinearAlgebra.eigfact
+LinearAlgebra.eigen
 LinearAlgebra.eigfact!
 LinearAlgebra.hessenberg
 LinearAlgebra.hessenberg!
diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
index 87e420480cdaa..e68f3c2126676 100644
--- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl
+++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
@@ -81,7 +81,7 @@ export
     diagm,
     dot,
     eig,
-    eigfact,
+    eigen,
     eigfact!,
     eigmax,
     eigmin,
diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl
index 4a1f4dc5aa8d1..4ec329b8439f6 100644
--- a/stdlib/LinearAlgebra/src/bidiag.jl
+++ b/stdlib/LinearAlgebra/src/bidiag.jl
@@ -625,4 +625,4 @@ function eigvecs(M::Bidiagonal{T}) where T
     end
     Q #Actually Triangular
 end
-eigfact(M::Bidiagonal) = Eigen(eigvals(M), eigvecs(M))
+eigen(M::Bidiagonal) = Eigen(eigvals(M), eigvecs(M))
diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl
index 5477e3d7693df..312ef17fb3c0a 100644
--- a/stdlib/LinearAlgebra/src/dense.jl
+++ b/stdlib/LinearAlgebra/src/dense.jl
@@ -478,7 +478,7 @@ Compute the matrix exponential of `A`, defined by
 e^A = \\sum_{n=0}^{\\infty} \\frac{A^n}{n!}.
 ```
 
-For symmetric or Hermitian `A`, an eigendecomposition ([`eigfact`](@ref)) is
+For symmetric or Hermitian `A`, an eigendecomposition ([`eigen`](@ref)) is
 used, otherwise the scaling and squaring algorithm (see [^H05]) is chosen.
 
 [^H05]: Nicholas J. Higham, "The squaring and scaling method for the matrix exponential revisited", SIAM Journal on Matrix Analysis and Applications, 26(4), 2005, 1179-1193. [doi:10.1137/090768539](https://doi.org/10.1137/090768539)
@@ -602,7 +602,7 @@ the unique matrix ``X`` such that ``e^X = A`` and ``-\\pi < Im(\\lambda) < \\pi`
 the eigenvalues ``\\lambda`` of ``X``. If `A` has nonpositive eigenvalues, a nonprincipal
 matrix function is returned whenever possible.
 
-If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is
+If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is
 used, if `A` is triangular an improved version of the inverse scaling and squaring method is
 employed (see [^AH12] and [^AHR13]). For general matrices, the complex Schur form
 ([`schur`](@ref)) is computed and the triangular algorithm is used on the
@@ -660,7 +660,7 @@ If `A` has no negative real eigenvalues, compute the principal matrix square roo
 that is the unique matrix ``X`` with eigenvalues having positive real part such that
 ``X^2 = A``. Otherwise, a nonprincipal square root is returned.
 
-If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is
+If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is
 used to compute the square root. Otherwise, the square root is determined by means of the
 Björck-Hammarling method [^BH83], which computes the complex Schur form ([`schur`](@ref))
 and then the complex square root of the triangular factor.
@@ -732,7 +732,7 @@ end
 
 Compute the matrix cosine of a square matrix `A`.
 
-If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is used to
+If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to
 compute the cosine. Otherwise, the cosine is determined by calling [`exp`](@ref).
 
 # Examples
@@ -765,7 +765,7 @@ end
 
 Compute the matrix sine of a square matrix `A`.
 
-If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is used to
+If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to
 compute the sine. Otherwise, the sine is determined by calling [`exp`](@ref).
 
 # Examples
@@ -851,7 +851,7 @@ end
 
 Compute the matrix tangent of a square matrix `A`.
 
-If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is used to
+If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to
 compute the tangent. Otherwise, the tangent is determined by calling [`exp`](@ref).
 
 # Examples
@@ -924,7 +924,7 @@ end
 
 Compute the inverse matrix cosine of a square matrix `A`.
 
-If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is used to
+If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to
 compute the inverse cosine. Otherwise, the inverse cosine is determined by using
 [`log`](@ref) and [`sqrt`](@ref).  For the theory and logarithmic formulas used to compute
 this function, see [^AH16_1].
@@ -955,7 +955,7 @@ end
 
 Compute the inverse matrix sine of a square matrix `A`.
 
-If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is used to
+If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to
 compute the inverse sine. Otherwise, the inverse sine is determined by using [`log`](@ref)
 and [`sqrt`](@ref).  For the theory and logarithmic formulas used to compute this function,
 see [^AH16_2].
@@ -986,7 +986,7 @@ end
 
 Compute the inverse matrix tangent of a square matrix `A`.
 
-If `A` is symmetric or Hermitian, its eigendecomposition ([`eigfact`](@ref)) is used to
+If `A` is symmetric or Hermitian, its eigendecomposition ([`eigen`](@ref)) is used to
 compute the inverse tangent. Otherwise, the inverse tangent is determined by using
 [`log`](@ref).  For the theory and logarithmic formulas used to compute this function, see
 [^AH16_3].
diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl
index f513d8d3af35e..08208769bc9ea 100644
--- a/stdlib/LinearAlgebra/src/deprecated.jl
+++ b/stdlib/LinearAlgebra/src/deprecated.jl
@@ -1359,3 +1359,22 @@ export hessfact
 # deprecate hessfact! to hessenberg!
 export hessenberg!
 @deprecate(hessfact!(A::StridedMatrix{<:BlasFloat}), hessenberg!(A))
+
+# deprecate eigfact to eigen
+export eigfact
+@deprecate(eigfact(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T, eigen(A; permute=permute, scale=scale))
+@deprecate(eigfact(x::Number), eigen(x))
+@deprecate(eigfact(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) where {TA,TB}, eigen(A, B))
+@deprecate(eigfact(A::Number, B::Number), eigen(A, B))
+
+@deprecate(eigfact(A::SymTridiagonal{T}) where T, eigen(A))
+@deprecate(eigfact(A::SymTridiagonal{T}, irange::UnitRange) where T, eigen(A))
+@deprecate(eigfact(A::SymTridiagonal{T}, vl::Real, vu::Real) where T, eigen(A))
+
+@deprecate(eigfact(M::Bidiagonal), eigen(M))
+
+@deprecate(eigfact(A::RealHermSymComplexHerm), eigen(A))
+@deprecate(eigfact(A::RealHermSymComplexHerm, irange::UnitRange), eigen(A, irange))
+@deprecate(eigfact(A::RealHermSymComplexHerm, vl::Real, vh::Real), eigen(A, vl, vh))
+@deprecate(eigfact(A::AbstractTriangular), eigen(A))
+@deprecate(eigfact(D::Diagonal; permute::Bool=true, scale::Bool=true), eigen(D; permute=permute, scale=scale))
diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl
index 39e0c8451c119..5466c0c1266c8 100644
--- a/stdlib/LinearAlgebra/src/diagonal.jl
+++ b/stdlib/LinearAlgebra/src/diagonal.jl
@@ -445,7 +445,7 @@ end
 eigvals(D::Diagonal{<:Number}) = D.diag
 eigvals(D::Diagonal) = [eigvals(x) for x in D.diag] #For block matrices, etc.
 eigvecs(D::Diagonal) = Matrix{eltype(D)}(I, size(D))
-function eigfact(D::Diagonal; permute::Bool=true, scale::Bool=true)
+function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true)
     if any(!isfinite, D.diag)
         throw(ArgumentError("matrix contains Infs or NaNs"))
     end
diff --git a/stdlib/LinearAlgebra/src/eigen.jl b/stdlib/LinearAlgebra/src/eigen.jl
index 56890060f13f2..07bc5bafa691d 100644
--- a/stdlib/LinearAlgebra/src/eigen.jl
+++ b/stdlib/LinearAlgebra/src/eigen.jl
@@ -20,12 +20,24 @@ end
 GeneralizedEigen(values::AbstractVector{V}, vectors::AbstractMatrix{T}) where {T,V} =
     GeneralizedEigen{T,V,typeof(vectors),typeof(values)}(values, vectors)
 
+# iteration for destructuring into components
+Base.iterate(S::Union{Eigen,GeneralizedEigen}) = (S.values, Val(:vectors))
+Base.iterate(S::Union{Eigen,GeneralizedEigen}, ::Val{:vectors}) = (S.vectors, Val(:done))
+Base.iterate(S::Union{Eigen,GeneralizedEigen}, ::Val{:done}) = nothing
+
+# indexing for destructuring into components
+@inline function Base.getindex(S::Union{Eigen,GeneralizedEigen}, i::Integer)
+    i == 1 ? (return S.values) :
+    i == 2 ? (return S.vectors) :
+        throw(BoundsError(S, i))
+end
+
 isposdef(A::Union{Eigen,GeneralizedEigen}) = isreal(A.values) && all(x -> x > 0, A.values)
 
 """
     eigfact!(A, [B])
 
-Same as [`eigfact`](@ref), but saves space by overwriting the input `A` (and
+Same as [`eigen`](@ref), but saves space by overwriting the input `A` (and
 `B`), instead of creating a copy.
 """
 function eigfact!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T<:BlasReal
@@ -59,12 +71,14 @@ function eigfact!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) whe
 end
 
 """
-    eigfact(A; permute::Bool=true, scale::Bool=true) -> Eigen
+    eigen(A; permute::Bool=true, scale::Bool=true) -> Eigen
 
 Computes the eigenvalue decomposition of `A`, returning an `Eigen` factorization object `F`
 which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the
 matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.)
 
+Iterating the decomposition produces the components `F.values` and `F.vectors`.
+
 The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).
 
 For general nonsymmetric matrices it is possible to specify how the matrix is balanced
@@ -74,7 +88,7 @@ make rows and columns more equal in norm. The default is `true` for both options
 
 # Examples
 ```jldoctest
-julia> F = eigfact([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
+julia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
 Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
 eigenvalues:
 3-element Array{Float64,1}:
@@ -98,17 +112,22 @@ julia> F.vectors
  1.0  0.0  0.0
  0.0  1.0  0.0
  0.0  0.0  1.0
+
+julia> vals, vecs = F; # destructuring via iteration
+
+julia> vals == F.values && vecs == F.vectors
+true
 ```
 """
-function eigfact(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T
+function eigen(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T
     AA = copy_oftype(A, eigtype(T))
-    isdiag(AA) && return eigfact(Diagonal(AA), permute = permute, scale = scale)
+    isdiag(AA) && return eigen(Diagonal(AA), permute = permute, scale = scale)
     return eigfact!(AA, permute = permute, scale = scale)
 end
-eigfact(x::Number) = Eigen([x], fill(one(x), 1, 1))
+eigen(x::Number) = Eigen([x], fill(one(x), 1, 1))
 
 function eig(A::Union{Number, StridedMatrix}; permute::Bool=true, scale::Bool=true)
-    F = eigfact(A, permute=permute, scale=scale)
+    F = eigen(A, permute=permute, scale=scale)
     F.values, F.vectors
 end
 
@@ -118,7 +137,7 @@ end
     eig(A, permute::Bool=true, scale::Bool=true) -> D, V
 
 Computes eigenvalues (`D`) and eigenvectors (`V`) of `A`.
-See [`eigfact`](@ref) for details on the
+See [`eigen`](@ref) for details on the
 `irange`, `vl`, and `vu` arguments
 (for [`SymTridiagonal`](@ref), [`Hermitian`](@ref), and
 [`Symmetric`](@ref) matrices)
@@ -131,11 +150,11 @@ julia> eig([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
 ([1.0, 3.0, 18.0], [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0])
 ```
 
-`eig` is a wrapper around [`eigfact`](@ref), extracting all parts of the
-factorization to a tuple; where possible, using [`eigfact`](@ref) is recommended.
+`eig` is a wrapper around [`eigen`](@ref), extracting all parts of the
+factorization to a tuple; where possible, using [`eigen`](@ref) is recommended.
 """
 function eig(A::AbstractMatrix, args...)
-    F = eigfact(A, args...)
+    F = eigen(A, args...)
     F.values, F.vectors
 end
 
@@ -144,7 +163,7 @@ end
 
 Return a matrix `M` whose columns are the eigenvectors of `A`. (The `k`th eigenvector can
 be obtained from the slice `M[:, k]`.) The `permute` and `scale` keywords are the same as
-for [`eigfact`](@ref).
+for [`eigen`](@ref).
 
 # Examples
 ```jldoctest
@@ -156,7 +175,7 @@ julia> eigvecs([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
 ```
 """
 eigvecs(A::Union{Number, AbstractMatrix}; permute::Bool=true, scale::Bool=true) =
-    eigvecs(eigfact(A, permute=permute, scale=scale))
+    eigvecs(eigen(A, permute=permute, scale=scale))
 eigvecs(F::Union{Eigen, GeneralizedEigen}) = F.vectors
 
 eigvals(F::Union{Eigen, GeneralizedEigen}) = F.values
@@ -351,13 +370,15 @@ function eigfact!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasComplex
 end
 
 """
-    eigfact(A, B) -> GeneralizedEigen
+    eigen(A, B) -> GeneralizedEigen
 
 Computes the generalized eigenvalue decomposition of `A` and `B`, returning a
 `GeneralizedEigen` factorization object `F` which contains the generalized eigenvalues in
 `F.values` and the generalized eigenvectors in the columns of the matrix `F.vectors`.
 (The `k`th generalized eigenvector can be obtained from the slice `F.vectors[:, k]`.)
 
+Iterating the decomposition produces the components `F.values` and `F.vectors`.
+
 # Examples
 ```jldoctest
 julia> A = [1 0; 0 -1]
@@ -370,7 +391,7 @@ julia> B = [0 1; 1 0]
  0  1
  1  0
 
-julia> F = eigfact(A, B);
+julia> F = eigen(A, B);
 
 julia> F.values
 2-element Array{Complex{Float64},1}:
@@ -381,22 +402,27 @@ julia> F.vectors
 2×2 Array{Complex{Float64},2}:
   0.0-1.0im   0.0+1.0im
  -1.0-0.0im  -1.0+0.0im
+
+julia> vals, vecs = F; # destructuring via iteration
+
+julia> vals == F.values && vecs == F.vectors
+true
 ```
 """
-function eigfact(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) where {TA,TB}
+function eigen(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) where {TA,TB}
     S = promote_type(eigtype(TA),TB)
     return eigfact!(copy_oftype(A, S), copy_oftype(B, S))
 end
 
-eigfact(A::Number, B::Number) = eigfact(fill(A,1,1), fill(B,1,1))
+eigen(A::Number, B::Number) = eigen(fill(A,1,1), fill(B,1,1))
 
 """
     eig(A, B) -> D, V
 
 Computes generalized eigenvalues (`D`) and vectors (`V`) of `A` with respect to `B`.
 
-`eig` is a wrapper around [`eigfact`](@ref), extracting all parts of the
-factorization to a tuple; where possible, using [`eigfact`](@ref) is recommended.
+`eig` is a wrapper around [`eigen`](@ref), extracting all parts of the
+factorization to a tuple; where possible, using [`eigen`](@ref) is recommended.
 
 # Examples
 ```jldoctest
@@ -415,11 +441,11 @@ julia> eig(A, B)
 ```
 """
 function eig(A::AbstractMatrix, B::AbstractMatrix)
-    F = eigfact(A,B)
+    F = eigen(A,B)
     F.values, F.vectors
 end
 function eig(A::Number, B::Number)
-    F = eigfact(A,B)
+    F = eigen(A,B)
     F.values, F.vectors
 end
 
@@ -524,7 +550,7 @@ julia> eigvecs(A, B)
  -1.0-0.0im  -1.0+0.0im
 ```
 """
-eigvecs(A::AbstractMatrix, B::AbstractMatrix) = eigvecs(eigfact(A, B))
+eigvecs(A::AbstractMatrix, B::AbstractMatrix) = eigvecs(eigen(A, B))
 
 function show(io::IO, mime::MIME{Symbol("text/plain")}, F::Union{Eigen,GeneralizedEigen})
     println(io, summary(F))
diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl
index 8967f0919a994..3a0bcf2a2871a 100644
--- a/stdlib/LinearAlgebra/src/symmetric.jl
+++ b/stdlib/LinearAlgebra/src/symmetric.jl
@@ -476,7 +476,7 @@ inv(A::Symmetric{<:Any,<:StridedMatrix}) = Symmetric(_inv(A), Symbol(A.uplo))
 
 eigfact!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)...)
 
-function eigfact(A::RealHermSymComplexHerm)
+function eigen(A::RealHermSymComplexHerm)
     T = eltype(A)
     S = eigtype(T)
     eigfact!(S != T ? convert(AbstractMatrix{S}, A) : copy(A))
@@ -485,11 +485,13 @@ end
 eigfact!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, irange::UnitRange) = Eigen(LAPACK.syevr!('V', 'I', A.uplo, A.data, 0.0, 0.0, irange.start, irange.stop, -1.0)...)
 
 """
-    eigfact(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> Eigen
+    eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> Eigen
 
 Computes the eigenvalue decomposition of `A`, returning an `Eigen` factorization object `F`
-which contains the eigenvalues in `F[:values]` and the eigenvectors in the columns of the
-matrix `F[:vectors]`. (The `k`th eigenvector can be obtained from the slice `F[:vectors][:, k]`.)
+which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the
+matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.)
+
+Iterating the decomposition produces the components `F.values` and `F.vectors`.
 
 The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).
 
@@ -499,7 +501,7 @@ The `UnitRange` `irange` specifies indices of the sorted eigenvalues to search f
     If `irange` is not `1:n`, where `n` is the dimension of `A`, then the returned factorization
     will be a *truncated* factorization.
 """
-function eigfact(A::RealHermSymComplexHerm, irange::UnitRange)
+function eigen(A::RealHermSymComplexHerm, irange::UnitRange)
     T = eltype(A)
     S = eigtype(T)
     eigfact!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange)
@@ -509,11 +511,13 @@ eigfact!(A::RealHermSymComplexHerm{T,<:StridedMatrix}, vl::Real, vh::Real) where
     Eigen(LAPACK.syevr!('V', 'V', A.uplo, A.data, convert(T, vl), convert(T, vh), 0, 0, -1.0)...)
 
 """
-    eigfact(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> Eigen
+    eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> Eigen
 
 Computes the eigenvalue decomposition of `A`, returning an `Eigen` factorization object `F`
-which contains the eigenvalues in `F[:values]` and the eigenvectors in the columns of the
-matrix `F[:vectors]`. (The `k`th eigenvector can be obtained from the slice `F[:vectors][:, k]`.)
+which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the
+matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]`.)
+
+Iterating the decomposition produces the components `F.values` and `F.vectors`.
 
 The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`](@ref), and [`isposdef`](@ref).
 
@@ -523,7 +527,7 @@ The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`
     If [`vl`, `vu`] does not contain all eigenvalues of `A`, then the returned factorization
     will be a *truncated* factorization.
 """
-function eigfact(A::RealHermSymComplexHerm, vl::Real, vh::Real)
+function eigen(A::RealHermSymComplexHerm, vl::Real, vh::Real)
     T = eltype(A)
     S = eigtype(T)
     eigfact!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh)
@@ -634,7 +638,7 @@ eigvals!(A::HermOrSym{T,S}, B::HermOrSym{T,S}) where {T<:BlasReal,S<:StridedMatr
 eigvals!(A::Hermitian{T,S}, B::Hermitian{T,S}) where {T<:BlasComplex,S<:StridedMatrix} =
     LAPACK.sygvd!(1, 'N', A.uplo, A.data, B.uplo == A.uplo ? B.data : copy(B.data'))[1]
 
-eigvecs(A::HermOrSym) = eigvecs(eigfact(A))
+eigvecs(A::HermOrSym) = eigvecs(eigen(A))
 
 function svdvals!(A::RealHermSymComplexHerm)
     vals = eigvals!(A)
@@ -656,7 +660,7 @@ function sympow(A::Symmetric, p::Integer)
 end
 function ^(A::Symmetric{<:Real}, p::Real)
     isinteger(p) && return integerpow(A, p)
-    F = eigfact(A)
+    F = eigen(A)
     if all(λ -> λ ≥ 0, F.values)
         return Symmetric((F.vectors * Diagonal((F.values).^p)) * F.vectors')
     else
@@ -680,7 +684,7 @@ function ^(A::Hermitian, p::Integer)
 end
 function ^(A::Hermitian{T}, p::Real) where T
     isinteger(p) && return integerpow(A, p)
-    F = eigfact(A)
+    F = eigen(A)
     if all(λ -> λ ≥ 0, F.values)
         retmat = (F.vectors * Diagonal((F.values).^p)) * F.vectors'
         if T <: Real
@@ -699,12 +703,12 @@ end
 for func in (:exp, :cos, :sin, :tan, :cosh, :sinh, :tanh, :atan, :asinh, :atanh)
     @eval begin
         function ($func)(A::HermOrSym{<:Real})
-            F = eigfact(A)
+            F = eigen(A)
             return Symmetric((F.vectors * Diagonal(($func).(F.values))) * F.vectors')
         end
         function ($func)(A::Hermitian{<:Complex})
             n = checksquare(A)
-            F = eigfact(A)
+            F = eigen(A)
             retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
             for i = 1:n
                 retmat[i,i] = real(retmat[i,i])
@@ -717,7 +721,7 @@ end
 for func in (:acos, :asin)
     @eval begin
         function ($func)(A::HermOrSym{<:Real})
-            F = eigfact(A)
+            F = eigen(A)
             if all(λ -> -1 ≤ λ ≤ 1, F.values)
                 retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
             else
@@ -727,7 +731,7 @@ for func in (:acos, :asin)
         end
         function ($func)(A::Hermitian{<:Complex})
             n = checksquare(A)
-            F = eigfact(A)
+            F = eigen(A)
             if all(λ -> -1 ≤ λ ≤ 1, F.values)
                 retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
                 for i = 1:n
@@ -742,7 +746,7 @@ for func in (:acos, :asin)
 end
 
 function acosh(A::HermOrSym{<:Real})
-    F = eigfact(A)
+    F = eigen(A)
     if all(λ -> λ ≥ 1, F.values)
         retmat = (F.vectors * Diagonal(acosh.(F.values))) * F.vectors'
     else
@@ -752,7 +756,7 @@ function acosh(A::HermOrSym{<:Real})
 end
 function acosh(A::Hermitian{<:Complex})
     n = checksquare(A)
-    F = eigfact(A)
+    F = eigen(A)
     if all(λ -> λ ≥ 1, F.values)
         retmat = (F.vectors * Diagonal(acosh.(F.values))) * F.vectors'
         for i = 1:n
@@ -766,7 +770,7 @@ end
 
 function sincos(A::HermOrSym{<:Real})
     n = checksquare(A)
-    F = eigfact(A)
+    F = eigen(A)
     S, C = Diagonal(similar(A, (n,))), Diagonal(similar(A, (n,)))
     for i in 1:n
         S.diag[i], C.diag[i] = sincos(F.values[i])
@@ -775,7 +779,7 @@ function sincos(A::HermOrSym{<:Real})
 end
 function sincos(A::Hermitian{<:Complex})
     n = checksquare(A)
-    F = eigfact(A)
+    F = eigen(A)
     S, C = Diagonal(similar(A, (n,))), Diagonal(similar(A, (n,)))
     for i in 1:n
         S.diag[i], C.diag[i] = sincos(F.values[i])
@@ -792,7 +796,7 @@ end
 for func in (:log, :sqrt)
     @eval begin
         function ($func)(A::HermOrSym{<:Real})
-            F = eigfact(A)
+            F = eigen(A)
             if all(λ -> λ ≥ 0, F.values)
                 retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
             else
@@ -803,7 +807,7 @@ for func in (:log, :sqrt)
 
         function ($func)(A::Hermitian{<:Complex})
             n = checksquare(A)
-            F = eigfact(A)
+            F = eigen(A)
             if all(λ -> λ ≥ 0, F.values)
                 retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
                 for i = 1:n
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index d77082425d394..9b0755cb2626d 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -2432,7 +2432,7 @@ function logabsdet(A::Union{UpperTriangular{T},LowerTriangular{T}}) where T
     return abs_det, sgn
 end
 
-eigfact(A::AbstractTriangular) = Eigen(eigvals(A), eigvecs(A))
+eigen(A::AbstractTriangular) = Eigen(eigvals(A), eigvecs(A))
 
 # Generic singular systems
 for func in (:svd, :svd!, :svdvals)
diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl
index 2be5dbd7ad9d5..6e5bd0baa3051 100644
--- a/stdlib/LinearAlgebra/src/tridiag.jl
+++ b/stdlib/LinearAlgebra/src/tridiag.jl
@@ -184,16 +184,16 @@ end
 (\)(T::SymTridiagonal, B::StridedVecOrMat) = ldlt(T)\B
 
 eigfact!(A::SymTridiagonal{<:BlasReal}) = Eigen(LAPACK.stegr!('V', A.dv, A.ev)...)
-eigfact(A::SymTridiagonal{T}) where T = eigfact!(copy_oftype(A, eigtype(T)))
+eigen(A::SymTridiagonal{T}) where T = eigfact!(copy_oftype(A, eigtype(T)))
 
 eigfact!(A::SymTridiagonal{<:BlasReal}, irange::UnitRange) =
     Eigen(LAPACK.stegr!('V', 'I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop)...)
-eigfact(A::SymTridiagonal{T}, irange::UnitRange) where T =
+eigen(A::SymTridiagonal{T}, irange::UnitRange) where T =
     eigfact!(copy_oftype(A, eigtype(T)), irange)
 
 eigfact!(A::SymTridiagonal{<:BlasReal}, vl::Real, vu::Real) =
     Eigen(LAPACK.stegr!('V', 'V', A.dv, A.ev, vl, vu, 0, 0)...)
-eigfact(A::SymTridiagonal{T}, vl::Real, vu::Real) where T =
+eigen(A::SymTridiagonal{T}, vl::Real, vu::Real) where T =
     eigfact!(copy_oftype(A, eigtype(T)), vl, vu)
 
 eigvals!(A::SymTridiagonal{<:BlasReal}) = LAPACK.stev!('N', A.dv, A.ev)[1]
@@ -214,7 +214,7 @@ eigmax(A::SymTridiagonal) = eigvals(A, size(A, 1):size(A, 1))[1]
 eigmin(A::SymTridiagonal) = eigvals(A, 1:1)[1]
 
 #Compute selected eigenvectors only corresponding to particular eigenvalues
-eigvecs(A::SymTridiagonal) = eigfact(A).vectors
+eigvecs(A::SymTridiagonal) = eigen(A).vectors
 
 """
     eigvecs(A::SymTridiagonal[, eigvals]) -> Matrix
diff --git a/stdlib/LinearAlgebra/test/dense.jl b/stdlib/LinearAlgebra/test/dense.jl
index 2628165945936..2e809075c4688 100644
--- a/stdlib/LinearAlgebra/test/dense.jl
+++ b/stdlib/LinearAlgebra/test/dense.jl
@@ -581,19 +581,19 @@ end
             @test all(z -> (0 < real(z) < π ||
                             abs(real(z)) < abstol && imag(z) >= 0 ||
                             abs(real(z) - π) < abstol && imag(z) <= 0),
-                      eigfact(acos(A)).values)
+                      eigen(acos(A)).values)
             @test all(z -> (-π/2 < real(z) < π/2 ||
                             abs(real(z) + π/2) < abstol && imag(z) >= 0 ||
                             abs(real(z) - π/2) < abstol && imag(z) <= 0),
-                      eigfact(asin(A)).values)
+                      eigen(asin(A)).values)
             @test all(z -> (-π < imag(z) < π && real(z) > 0 ||
                             0 <= imag(z) < π && abs(real(z)) < abstol ||
                             abs(imag(z) - π) < abstol && real(z) >= 0),
-                      eigfact(acosh(A)).values)
+                      eigen(acosh(A)).values)
             @test all(z -> (-π/2 < imag(z) < π/2 ||
                             abs(imag(z) + π/2) < abstol && real(z) <= 0 ||
                             abs(imag(z) - π/2) < abstol && real(z) <= 0),
-                      eigfact(asinh(A)).values)
+                      eigen(asinh(A)).values)
         end
     end
 end
diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl
index 651e14b819115..b3467621fcc10 100644
--- a/stdlib/LinearAlgebra/test/diagonal.jl
+++ b/stdlib/LinearAlgebra/test/diagonal.jl
@@ -207,7 +207,7 @@ srand(1)
     @test factorize(D) == D
 
     @testset "Eigensystem" begin
-        eigD = eigfact(D)
+        eigD = eigen(D)
         @test Diagonal(eigD.values) ≈ D
         @test eigD.vectors == Matrix(I, size(D))
     end
diff --git a/stdlib/LinearAlgebra/test/eigen.jl b/stdlib/LinearAlgebra/test/eigen.jl
index eec694ed7f475..c4e3c087e4b3d 100644
--- a/stdlib/LinearAlgebra/test/eigen.jl
+++ b/stdlib/LinearAlgebra/test/eigen.jl
@@ -37,7 +37,7 @@ aimg  = randn(n,n)/2
             for i in 1:size(a,2)
                 @test a*v[:,i] ≈ d[i]*v[:,i]
             end
-            f = eigfact(a)
+            f = eigen(a)
             @test det(a) ≈ det(f)
             @test inv(a) ≈ inv(f)
             @test isposdef(a) == isposdef(f)
@@ -45,7 +45,7 @@ aimg  = randn(n,n)/2
             @test eigvecs(f) === f.vectors
             @test Array(f) ≈ a
 
-            num_fact = eigfact(one(eltya))
+            num_fact = eigen(one(eltya))
             @test num_fact.values[1] == one(eltya)
             h = asym
             @test minimum(eigvals(h)) ≈ eigmin(h)
@@ -61,7 +61,7 @@ aimg  = randn(n,n)/2
                 asym_sg = view(asym, 1:n1, 1:n1)
                 a_sg = view(a, 1:n, n1+1:n2)
             end
-            f = eigfact(asym_sg, a_sg'a_sg)
+            f = eigen(asym_sg, a_sg'a_sg)
             @test asym_sg*f.vectors ≈ (a_sg'a_sg*f.vectors) * Diagonal(f.values)
             @test f.values ≈ eigvals(asym_sg, a_sg'a_sg)
             @test prod(f.values) ≈ prod(eigvals(asym_sg/(a_sg'a_sg))) atol=200ε
@@ -82,7 +82,7 @@ aimg  = randn(n,n)/2
                 a1_nsg = view(a, 1:n1, 1:n1)
                 a2_nsg = view(a, n1+1:n2, n1+1:n2)
             end
-            f = eigfact(a1_nsg, a2_nsg)
+            f = eigen(a1_nsg, a2_nsg)
             @test a1_nsg*f.vectors ≈ (a2_nsg*f.vectors) * Diagonal(f.values)
             @test f.values ≈ eigvals(a1_nsg, a2_nsg)
             @test prod(f.values) ≈ prod(eigvals(a1_nsg/a2_nsg)) atol=50000ε
@@ -109,7 +109,7 @@ 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))
-        f = eigfact(a)
+        f = eigen(a)
         @test a ≈ f.vectors * Diagonal(f.values) / f.vectors
     end
 end
@@ -124,8 +124,8 @@ end
 
 @testset "text/plain (REPL) printing of Eigen and GeneralizedEigen" begin
     A, B = randn(5,5), randn(5,5)
-    e    = eigfact(A)
-    ge   = eigfact(A, B)
+    e    = eigen(A)
+    ge   = eigen(A, B)
     valsstring = sprint((t, s) -> show(t, "text/plain", s), e.values)
     vecsstring = sprint((t, s) -> show(t, "text/plain", s), e.vectors)
     factstring = sprint((t, s) -> show(t, "text/plain", s), e)
diff --git a/stdlib/LinearAlgebra/test/lapack.jl b/stdlib/LinearAlgebra/test/lapack.jl
index 642827da3082a..9289a0b3efbc9 100644
--- a/stdlib/LinearAlgebra/test/lapack.jl
+++ b/stdlib/LinearAlgebra/test/lapack.jl
@@ -231,7 +231,7 @@ end
     @testset for elty in (ComplexF32, ComplexF64)
         A = rand(elty,10,10)
         Aw, Avl, Avr = LAPACK.geev!('N','V',copy(A))
-        fA = eigfact(A)
+        fA = eigen(A)
         @test fA.values  ≈ Aw
         @test fA.vectors ≈ Avr
     end
@@ -660,7 +660,7 @@ end
 
 # Issue 14065 (and 14220)
 let A = [NaN NaN; NaN NaN]
-    @test_throws ArgumentError eigfact(A)
+    @test_throws ArgumentError eigen(A)
 end
 
 end # module TestLAPACK
diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl
index 96eab7b13fc88..4b4353f6e4d36 100644
--- a/stdlib/LinearAlgebra/test/lu.jl
+++ b/stdlib/LinearAlgebra/test/lu.jl
@@ -50,10 +50,10 @@ dimg  = randn(n)/2
                                        -2.0      4.0       1.0    -eps(real(one(eltya)));
                                        -eps(real(one(eltya)))/4  eps(real(one(eltya)))/2  -1.0     0;
                                        -0.5     -0.5       0.1     1.0])
-            F = eigfact(A, permute=false, scale=false)
+            F = eigen(A, permute=false, scale=false)
             eig(A, permute=false, scale=false)
             @test F.vectors*Diagonal(F.values)/F.vectors ≈ A
-            F = eigfact(A)
+            F = eigen(A)
             # @test norm(F.vectors*Diagonal(F.values)/F.vectors - A) > 0.01
         end
     end
diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl
index 5ff60262510a6..54172dfcb2af7 100644
--- a/stdlib/LinearAlgebra/test/symmetric.jl
+++ b/stdlib/LinearAlgebra/test/symmetric.jl
@@ -222,14 +222,14 @@ end
                         @test asym*v[:,1] ≈ d[1]*v[:,1]
                         @test v*Diagonal(d)*transpose(v) ≈ asym
                         @test isequal(eigvals(asym[1]), eigvals(asym[1:1,1:1]))
-                        @test abs.(eigfact(Symmetric(asym), 1:2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
+                        @test abs.(eigen(Symmetric(asym), 1:2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
                         eig(Symmetric(asym), 1:2) # same result, but checks that method works
-                        @test abs.(eigfact(Symmetric(asym), d[1] - 1, (d[2] + d[3])/2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
+                        @test abs.(eigen(Symmetric(asym), d[1] - 1, (d[2] + d[3])/2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
                         eig(Symmetric(asym), d[1] - 1, (d[2] + d[3])/2) # same result, but checks that method works
                         @test eigvals(Symmetric(asym), 1:2) ≈ d[1:2]
                         @test eigvals(Symmetric(asym), d[1] - 1, (d[2] + d[3])/2) ≈ d[1:2]
-                        # eigfact doesn't support Symmetric{Complex}
-                        @test Matrix(eigfact(asym)) ≈ asym
+                        # eigen doesn't support Symmetric{Complex}
+                        @test Matrix(eigen(asym)) ≈ asym
                         @test eigvecs(Symmetric(asym)) ≈ eigvecs(asym)
                     end
 
@@ -237,13 +237,13 @@ end
                     @test aherm*v[:,1] ≈ d[1]*v[:,1]
                     @test v*Diagonal(d)*v' ≈ aherm
                     @test isequal(eigvals(aherm[1]), eigvals(aherm[1:1,1:1]))
-                    @test abs.(eigfact(Hermitian(aherm), 1:2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
+                    @test abs.(eigen(Hermitian(aherm), 1:2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
                     eig(Hermitian(aherm), 1:2) # same result, but checks that method works
-                    @test abs.(eigfact(Hermitian(aherm), d[1] - 1, (d[2] + d[3])/2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
+                    @test abs.(eigen(Hermitian(aherm), d[1] - 1, (d[2] + d[3])/2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
                     eig(Hermitian(aherm), d[1] - 1, (d[2] + d[3])/2) # same result, but checks that method works
                     @test eigvals(Hermitian(aherm), 1:2) ≈ d[1:2]
                     @test eigvals(Hermitian(aherm), d[1] - 1, (d[2] + d[3])/2) ≈ d[1:2]
-                    @test Matrix(eigfact(aherm)) ≈ aherm
+                    @test Matrix(eigen(aherm)) ≈ aherm
                     @test eigvecs(Hermitian(aherm)) ≈ eigvecs(aherm)
 
                     # relation to svdvals
@@ -365,7 +365,7 @@ end
 end
 
 @testset "Issues #8057 and #8058. f=$f, A=$A" for f in
-        (eigfact, eigvals, eig),
+        (eigen, eigvals, eig),
             A in (Symmetric([0 1; 1 0]), Hermitian([0 im; -im 0]))
     @test_throws ArgumentError f(A, 3, 2)
     @test_throws ArgumentError f(A, 1:4)
diff --git a/stdlib/LinearAlgebra/test/tridiag.jl b/stdlib/LinearAlgebra/test/tridiag.jl
index 679ed2dfd20a0..903d09a02d27f 100644
--- a/stdlib/LinearAlgebra/test/tridiag.jl
+++ b/stdlib/LinearAlgebra/test/tridiag.jl
@@ -253,14 +253,14 @@ end
                         test_approx_eq_vecs(v, evecs)
                     end
                     @testset "stegr! call with index range" begin
-                        F = eigfact(SymTridiagonal(b, a),1:2)
-                        fF = eigfact(Symmetric(Array(SymTridiagonal(b, a))),1:2)
+                        F = eigen(SymTridiagonal(b, a),1:2)
+                        fF = eigen(Symmetric(Array(SymTridiagonal(b, a))),1:2)
                         test_approx_eq_modphase(F.vectors, fF.vectors)
                         @test F.values ≈ fF.values
                     end
                     @testset "stegr! call with value range" begin
-                        F = eigfact(SymTridiagonal(b, a),0.0,1.0)
-                        fF = eigfact(Symmetric(Array(SymTridiagonal(b, a))),0.0,1.0)
+                        F = eigen(SymTridiagonal(b, a),0.0,1.0)
+                        fF = eigen(Symmetric(Array(SymTridiagonal(b, a))),0.0,1.0)
                         test_approx_eq_modphase(F.vectors, fF.vectors)
                         @test F.values ≈ fF.values
                     end