diff --git a/README.md b/README.md
index ec28fbbf7..ee999d622 100644
--- a/README.md
+++ b/README.md
@@ -241,26 +241,26 @@ end
 
 This solves the confined anharmonic oscillator, `[-𝒟² + V(x)] u = λu`, where `u(±10) = 0`, `V(x) = ω*x² + x⁴`, and `ω = 25`.
 ```julia
-    n = 3000
-    ω = 25.0
-    d = Segment(-10..10)
-    S = Ultraspherical(0.5, d)
-    NS = NormalizedPolynomialSpace(S)
-    V = Fun(x->ω*x^2+x^4, S)
-    L = -Derivative(S, 2) + V
-    C = Conversion(domainspace(L), rangespace(L))
-    B = Dirichlet(S)
-    QS = QuotientSpace(B)
-    Q = Conversion(QS, S)
-    D1 = Conversion(S, NS)
-    D2 = Conversion(NS, S)
-    R = D1*Q
-    P = cache(PartialInverseOperator(C, (0, ApproxFun.bandwidth(L, 1) + ApproxFun.bandwidth(R, 1) + ApproxFun.bandwidth(C, 2))))
-    A = R'D1*P*L*D2*R
-    B = R'R
-    SA = Symmetric(A[1:n,1:n], :L)
-    SB = Symmetric(B[1:n,1:n], :L)
-    λ = eigvals(SA, SB)[1:round(Int, 3n/5)]
+n = 3000
+ω = 25.0
+d = Segment(-10..10)
+S = Ultraspherical(0.5, d)
+NS = NormalizedPolynomialSpace(S)
+V = Fun(x->ω*x^2+x^4, S)
+L = -Derivative(S, 2) + V
+C = Conversion(domainspace(L), rangespace(L))
+B = Dirichlet(S)
+QS = QuotientSpace(B)
+Q = Conversion(QS, S)
+D1 = Conversion(S, NS)
+D2 = Conversion(NS, S)
+R = D1*Q
+P = cache(PartialInverseOperator(C, (0, bandwidth(L, 1) + bandwidth(R, 1) + bandwidth(C, 2))))
+A = R'D1*P*L*D2*R
+B = R'R
+SA = Symmetric(A[1:n,1:n], :L)
+SB = Symmetric(B[1:n,1:n], :L)
+λ = eigvals(SA, SB)[1:round(Int, 3n/5)]
 ```
 
 
diff --git a/src/ApproxFun.jl b/src/ApproxFun.jl
index 4f393329c..28cfe8b3d 100644
--- a/src/ApproxFun.jl
+++ b/src/ApproxFun.jl
@@ -27,14 +27,14 @@ import Base: values, convert, getindex, setindex!, *, +, -, ==, <, <=, >, |, !,
                 getproperty, findfirst, unsafe_getindex, fld, cld, div, real, imag,
                 @_inline_meta, eachindex, firstindex, lastindex, keys, isreal, OneTo,
                 Array, Vector, Matrix, view, ones, @propagate_inbounds, print_array,
-                split
+                split, checkbounds
 
 import Base.Broadcast: BroadcastStyle, Broadcasted, AbstractArrayStyle, broadcastable,
                         DefaultArrayStyle, broadcasted
 
 import Statistics: mean
 
-import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, det, eigvals, dot, cross,
+import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, lmul!, det, eigvals, dot, cross,
                         qr, qr!, rank, isdiag, istril, istriu, issymmetric, ishermitian,
                         Tridiagonal, diagm, diagm_container, factorize, nullspace,
                         Hermitian, Symmetric, adjoint, transpose, char_uplo
diff --git a/src/LinearAlgebra/LinearAlgebra.jl b/src/LinearAlgebra/LinearAlgebra.jl
index ef6992036..539a5b28d 100644
--- a/src/LinearAlgebra/LinearAlgebra.jl
+++ b/src/LinearAlgebra/LinearAlgebra.jl
@@ -19,3 +19,7 @@ include("RaggedMatrix.jl")
 
 
 include("rowvector.jl")
+
+include("deflationcheck.jl")
+include("SymBandedPlusBulge.jl")
+include("SymBlockArrowHead.jl")
diff --git a/src/LinearAlgebra/SymBandedPlusBulge.jl b/src/LinearAlgebra/SymBandedPlusBulge.jl
new file mode 100644
index 000000000..b6cf65749
--- /dev/null
+++ b/src/LinearAlgebra/SymBandedPlusBulge.jl
@@ -0,0 +1,178 @@
+"""
+    Represent a symmetric banded matrix plus a bulge:
+    [ □ ◺
+      ◹ □ ◺
+        ◹ □ □ ◺
+          □ □ □
+          ◹ □ □ ◺
+              ◹ □ ◺
+                ◹ □ ]
+"""
+mutable struct SymBandedPlusBulge{T,M<:BandedMatrix{T}} <: AbstractMatrix{T}
+    A::M
+    b::Int
+    bulge::Int
+end
+
+SymBandedPlusBulge(A::AbstractMatrix, b::Int, bulge::Int) = SymBandedPlusBulge(BandedMatrix(A, (2b, 2b)), b, bulge)
+
+size(A::SymBandedPlusBulge) = size(A.A)
+
+function getindex(A::SymBandedPlusBulge{T}, i::Integer, j::Integer) where T
+    b = A.b
+    bulge = A.bulge
+    col1 = bulge - 2b
+    col2 = bulge - b
+    row1 = col1 - b
+    row2 = col2 - 2b
+    AA = A.A
+    if -b ≤ i-j ≤ b
+        if j ≥ i
+            return AA[i,j]
+        else
+            return AA[j,i]
+        end
+    elseif 1 ≤ j-col1 ≤ b && 1 ≤ i-row1 ≤ b && j-col1 > i-row1
+        return AA[i,j]
+    elseif 1 ≤ j-col2 ≤ b && 1 ≤ i-row2-(j-col2-1) ≤ b
+        return AA[i,j]
+    elseif 1 ≤ i-col1 ≤ b && 1 ≤ j-row1 ≤ b && i-col1 > j-row1
+        return AA[j,i]
+    elseif 1 ≤ i-col2 ≤ b && 1 ≤ j-row2-(i-col2-1) ≤ b
+        return AA[j,i]
+    else
+        return zero(T)
+    end
+end
+
+setindex!(A::SymBandedPlusBulge{T}, v, i::Integer, j::Integer) where T = (b = A.b; -2b ≤ i-j ≤ 2b && setindex!(A.A, v, i, j))
+
+function computeHouseholder!(A::SymBandedPlusBulge{T}, w::Vector{T}, col::Int) where T
+    b = A.b
+    corr = zero(T)
+    fill!(w, corr)
+    row = col - b
+    @inbounds for i = max(row-b, 1):row
+        corr += abs2(A[i, col])
+        w[i] = A[i, col]
+    end
+    normw = sqrt(corr)
+    corr = copysign(normw, A[row, col])
+    w[row] += corr
+    normw = sqrt(normw^2+corr*muladd(T(2),w[row],-corr))
+    if normw == zero(normw)
+        return w
+    else
+        return LinearAlgebra.__normalize!(w, normw)
+    end
+end
+
+function applyHouseholder!(w::Vector{T}, wQ::Vector{T}, Q::Matrix{T}, bulge::Int, b::Int) where T
+    sw = max(bulge - 2b, 1)
+    fw = min(bulge - b, size(Q, 1))
+    s = 1
+    f = size(Q, 1)
+
+    # wQ = (w^⊤ Q)^⊤
+    fill!(wQ, zero(T))
+    @inbounds for j = s:f
+        wQj = zero(T)
+        @simd for i = sw:fw
+            wQj += w[i]*Q[i,j]
+        end
+        wQ[j] = wQj
+    end
+
+    # H Q = Q - 2/dot(w, w) (w (w^⊤ Q))
+    # twodivnrmw = 2/dot(w, w)
+    twodivnrmw = T(2)
+    @inbounds for j = s:f
+        wQj = wQ[j]
+        @simd for i = sw:fw
+            Q[i,j] -= twodivnrmw*w[i]*wQj
+        end
+    end
+    #Q .= Q - 2*w*(w'Q)
+
+    Q
+end
+
+function similarity!(w::Vector{T}, v::Vector{T}, Aw::Vector{T}, A::SymBandedPlusBulge{T}) where T
+    bulge = A.bulge
+    b = A.b
+    sw = max(bulge - 2b, 1)
+    fw = min(bulge - b, size(A, 1))
+    sv = max(bulge - 3b, 1)
+    fv = min(bulge, size(A, 1))
+    s = min(sv, sw)
+    f = max(fv, fw)
+
+    # Aw = A*w
+    fill!(Aw, zero(T))
+    @inbounds for j = sw:fw
+        wj = w[j]
+        @simd for i = sv:fv
+            Aw[i] += A[i,j]*wj
+        end
+    end
+
+    # v = Aw - dot(w, Aw)/dot(w, w)*w
+    fill!(v, zero(T))
+    # cst = -dot(w, Aw)/dot(w, w)
+    cst = -dot(w, Aw)
+    @inbounds for i = sv:fv
+        v[i] = muladd(cst, w[i], Aw[i])
+    end
+
+    # H A H^⊤ = A - 2/dot(w, w) ( vw^⊤ + wv^⊤ )
+    # twodivnrmw = 2/dot(w, w)
+    twodivnrmw = T(2)
+    @inbounds for j = s:f
+        vj = v[j]
+        wj = w[j]
+        @simd for i = s:f
+            A[i,j] -= twodivnrmw*(v[i]*wj+w[i]*vj)
+        end
+    end
+
+    A.bulge -= 1
+    A
+end
+
+function chasebulge!(w::Vector{T}, v::Vector{T}, Aw::Vector{T}, A::SymBandedPlusBulge{T}) where T
+    bulge = A.bulge
+    b = A.b
+    for k = bulge:-1:b+1
+        computeHouseholder!(A, w, k)
+        similarity!(w, v, Aw, A)
+    end
+
+    A
+end
+
+function chasebulge!(A::SymBandedPlusBulge{T}) where T
+    w = zeros(T, size(A, 2))
+    v = zero(w)
+    Aw = zero(w)
+    chasebulge!(w, v, Aw, A)
+end
+
+function chasebulge!(w::Vector{T}, v::Vector{T}, Aw::Vector{T}, A::SymBandedPlusBulge{T}, wQ::Vector{T}, Q::Matrix{T}) where T
+    bulge = A.bulge
+    b = A.b
+    for k = bulge:-1:b+1
+        computeHouseholder!(A, w, k)
+        applyHouseholder!(w, wQ, Q, k, b)
+        similarity!(w, v, Aw, A)
+    end
+
+    A, Q
+end
+
+function chasebulge!(A::SymBandedPlusBulge{T}, Q::Matrix{T}) where T
+    w = zeros(T, size(A, 2))
+    v = zero(w)
+    Aw = zero(w)
+    wQ = zero(w)
+    chasebulge!(w, v, Aw, A, wQ, Q)
+end
diff --git a/src/LinearAlgebra/SymBlockArrowHead.jl b/src/LinearAlgebra/SymBlockArrowHead.jl
new file mode 100644
index 000000000..b1555c7d2
--- /dev/null
+++ b/src/LinearAlgebra/SymBlockArrowHead.jl
@@ -0,0 +1,252 @@
+"""
+    Represent a symmetric block arrowhead matrix:
+    [A   B
+     B^⊤ C]
+"""
+struct SymBlockArrowHead{T,M<:AbstractMatrix{T}} <: AbstractMatrix{T}
+    A::M
+    B::Matrix{T}
+    C::Matrix{T}
+    function SymBlockArrowHead(A::M, B::Matrix{T}, C::Matrix{T}) where {T, M<:AbstractMatrix{T}}
+        @assert size(A, 1) == size(A, 2) == size(B, 1)
+        @assert size(B, 2) == size(C, 1) == size(C, 2)
+        new{T, M}(A, B, C)
+    end
+end
+
+size(A::SymBlockArrowHead) = (s = sum(size(A.B)); (s, s))
+
+function getindex(A::SymBlockArrowHead{T}, i::Integer, j::Integer) where T
+    @boundscheck checkbounds(A, i, j)
+    if i ≤ size(A.A, 1) && j ≤ size(A.A, 2)
+        A.A[i,j]
+    elseif i ≤ size(A.A, 1) && j > size(A.A, 2)
+        return A.B[i, j - size(A.A, 2)]
+    elseif j ≤ size(A.A, 2) && i > size(A.A, 1)
+        return A.B[j, i - size(A.A, 1)]
+    elseif size(A.A, 1) < i ≤ size(A, 1) && size(A.A, 2) < j ≤ size(A, 2)
+        return A.C[i - size(A.A, 1), j - size(A.A, 2)]
+    else
+        return zero(T)
+    end
+end
+
+function setindex!(A::SymBlockArrowHead{T}, v, i::Integer, j::Integer) where T
+    if i ≤ size(A.A, 1) && j ≤ size(A.A, 2)
+        A.A[i,j] = v
+    elseif i ≤ size(A.A, 1) && j > size(A.A, 2)
+        return A.B[i, j - size(A.A, 2)] = v
+    elseif j ≤ size(A.A, 2) && i > size(A.A, 1)
+        return A.B[j, i - size(A.A, 1)] = v
+    elseif size(A.A, 1) < i ≤ size(A, 1) && size(A.A, 2) < j ≤ size(A, 2)
+        return A.C[i - size(A.A, 1), j - size(A.A, 2)] = v
+    else
+        throw(error("Index out of reach."))
+    end
+end
+
+function computeHouseholder!(A::SymBlockArrowHead{T,SymBandedPlusBulge{T,V}}, w::Vector{T}, row::Int, col::Int) where {T,V}
+    b = A.A.b
+    corr = zero(T)
+    fill!(w, corr)
+    @inbounds for i = max(row-b, 1):row
+        corr += abs2(A[i, col])
+        w[i] = A[i, col]
+    end
+    normw = sqrt(corr)
+    corr = copysign(normw, A[row, col])
+    w[row] += corr
+    normw = sqrt(normw^2+corr*muladd(T(2),w[row],-corr))
+    if normw == zero(normw)
+        return w
+    else
+        return LinearAlgebra.__normalize!(w, normw)
+    end
+end
+
+function applyHouseholder!(w::Vector{T}, wQ::Vector{T}, Q::Matrix{T}, bulge::Int, b::Int, row::Int) where T
+    sw = max(bulge - b, 1)
+    fw = min(bulge, size(Q, 1))
+    s = 1
+    f = size(Q, 1)
+
+    # wQ = (w^⊤ Q)^⊤
+    fill!(wQ, zero(T))
+    @inbounds for j = s:f
+        wQj = zero(T)
+        @simd for i = sw:fw
+            wQj += w[i]*Q[i,j]
+        end
+        wQ[j] = wQj
+    end
+
+    # H Q = Q - 2/dot(w, w) (w (w^⊤ Q))
+    # twodivnrmw = 2/dot(w, w)
+    twodivnrmw = T(2)
+    @inbounds for j = s:f
+        wQj = wQ[j]
+        @simd for i = sw:fw
+            Q[i,j] -= twodivnrmw*w[i]*wQj
+        end
+    end
+    #Q .= Q - 2*w*(w'Q)
+
+    Q
+end
+
+function similarity!(w::Vector{T}, v::Vector{T}, Aw::Vector{T}, A::SymBlockArrowHead{T,SymBandedPlusBulge{T,V}}, row::Int) where {T,V}
+    bulge = A.A.bulge
+    b = A.A.b
+    sw = max(row-b, 1)
+    fw = min(row, size(A, 2))
+    sv = max(bulge - 3b, 1)
+    fv = min(bulge, size(A, 1))
+    sv1 = size(A, 2) - b + 1
+    fv1 = size(A, 2)
+    s = min(sv, sw)
+    f = max(fv, fw)
+
+    # Aw = A*w
+    fill!(Aw, zero(T))
+    @inbounds for j = sw:fw
+        wj = w[j]
+        @simd for i = sv:fv
+            Aw[i] += A[i,j]*wj
+        end
+        @simd for i = sv1:fv1
+            Aw[i] += A[i,j]*wj
+        end
+    end
+
+    # v = Aw - dot(w, Aw)/dot(w, w)*w
+    fill!(v, zero(T))
+    # cst = -dot(w, Aw)/dot(w, w)
+    cst = -dot(w, Aw)
+    @inbounds @simd for i = sv:fv
+        v[i] = muladd(cst, w[i], Aw[i])
+    end
+    @inbounds @simd for i = sv1:fv1
+        v[i] = muladd(cst, w[i], Aw[i])
+    end
+
+    # H A H^⊤ = A - 2/dot(w, w) ( vw^⊤ + wv^⊤ )
+    # twodivnrmw = 2/dot(w, w)
+    twodivnrmw = T(2)
+    @inbounds for j = s:f
+        vj = v[j]
+        wj = w[j]
+        @simd for i = s:f
+            A[i,j] -= twodivnrmw*(v[i]*wj+w[i]*vj)
+        end
+    end
+    @inbounds for j = sv1:fv1
+        vj = v[j]
+        wj = w[j]
+        @simd for i = s:f
+            A[i,j] -= twodivnrmw*(v[i]*wj+w[i]*vj)
+        end
+    end
+
+    A
+end
+
+function symblockarrowhead2symbanded!(w::Vector{T}, v::Vector{T}, Aw::Vector{T}, A::SymBlockArrowHead{T,SymBandedPlusBulge{T,V}}) where {T,V}
+    n = size(A, 2)
+    b = A.A.b
+
+    for j = 1:b-1
+        computeHouseholder!(A, w, b+1-j, n+1-j)
+        similarity!(w, v, Aw, A, b+1-j)
+    end
+
+    for m = 2:floor(Int, n÷b)-1
+        A.A.bulge = m*A.A.b
+        for j = 1:b
+            computeHouseholder!(A, w, m*b+1-j, n+1-j)
+            similarity!(w, v, Aw, A, m*b+1-j)
+        end
+        chasebulge!(w, v, Aw, A.A)
+    end
+
+    A
+end
+
+function symblockarrowhead2symbanded!(A::SymBlockArrowHead{T,SymBandedPlusBulge{T,V}}) where {T,V}
+    w = zeros(T, size(A, 2))
+    v = zero(w)
+    Aw = zero(w)
+    symblockarrowhead2symbanded!(w, v, Aw, A)
+end
+
+function symblockarrowhead2symbanded!(w::Vector{T}, v::Vector{T}, Aw::Vector{T}, A::SymBlockArrowHead{T,SymBandedPlusBulge{T,V}}, wQ::Vector{T}, Q::Matrix{T}) where {T,V}
+    n = size(A, 2)
+    b = A.A.b
+
+    for j = 1:b-1
+        computeHouseholder!(A, w, b+1-j, n+1-j)
+        applyHouseholder!(w, wQ, Q, b+1-j, b, n+1-j)
+        similarity!(w, v, Aw, A, b+1-j)
+    end
+
+    for m = 2:floor(Int, n÷b)-1
+        A.A.bulge = m*A.A.b
+        for j = 1:b
+            computeHouseholder!(A, w, m*b+1-j, n+1-j)
+            applyHouseholder!(w, wQ, Q, m*b+1-j, b, n+1-j)
+            similarity!(w, v, Aw, A, m*b+1-j)
+        end
+        chasebulge!(w, v, Aw, A.A, wQ, Q)
+    end
+
+    A
+end
+
+function symblockarrowhead2symbanded!(A::SymBlockArrowHead{T,SymBandedPlusBulge{T,V}}, Q::Matrix{T}) where {T,V}
+    w = zeros(T, size(A, 2))
+    v = zero(w)
+    Aw = zero(w)
+    wQ = zero(w)
+    symblockarrowhead2symbanded!(w, v, Aw, A, wQ, Q)
+end
+
+function structured_cholesky!(A::SymBlockArrowHead{T,Diagonal{T,V}}) where {T,V}
+    D = A.A
+    d = D.diag
+    d .= sqrt.(d)
+    B = A.B
+    ldiv!(D, B)
+    A.C .= cholesky(A.C-B'B).factors
+    UpperTriangular(A)
+end
+
+function structured_invcholesky!(A::SymBlockArrowHead{T, Diagonal{T,V}}) where {T,V}
+    D = A.A
+    d = D.diag
+    d .= sqrt.(d)
+    B = A.B
+    ldiv!(D, B)
+    A.C .= cholesky(A.C-B'B).factors
+    # A, B, C are the old D^{1/2}, E, F.
+    d .= inv.(d)
+    A.C .= inv(UpperTriangular(A.C)).data
+    B .= -lmul!(D, B*A.C)
+    UpperTriangular(A)
+end
+
+function congruence!(L::SymBlockArrowHead{T,M}, C::SymBlockArrowHead{T, Diagonal{T,V}}) where {T, M, V}
+    R = structured_invcholesky!(C)
+    A = L.A
+    B = L.B
+    CC = L.C
+    D = R.data.A
+    E = R.data.B
+    F = UpperTriangular(R.data.C)
+    DAD = D*A*D
+    DAEpDBF = D*(A*E + B*F)
+    EtBF = E'B*F
+    EtAEpEtBFpFtBtEpFtCF = E'A*E + F'CC*F + (EtBF+EtBF')
+    A .= DAD
+    B .= DAEpDBF
+    CC .= EtAEpEtBFpFtBtEpFtCF
+    L, C
+end
diff --git a/src/LinearAlgebra/deflationcheck.jl b/src/LinearAlgebra/deflationcheck.jl
new file mode 100644
index 000000000..2b988706a
--- /dev/null
+++ b/src/LinearAlgebra/deflationcheck.jl
@@ -0,0 +1,45 @@
+function deflationcheck(M::Matrix, Λ::Vector, tol::Real, verbose::Bool)
+    b, k = size(M)
+    @assert k == length(Λ)
+    maxΛ = maximum(abs, Λ)
+    minΛ = minimum(abs, Λ)
+    j = 1
+    maxΛ2j = (abs(Λ[j]) + minΛ)*(abs(Λ[j]) + maxΛ)
+    FroM2 = colnorm2(M, j)
+    while FroM2 ≤ tol*maxΛ2j && j < k
+        j += 1
+        maxΛ2j = max(maxΛ2j, (abs(Λ[j]) + minΛ)*(abs(Λ[j]) + maxΛ))
+        FroM2 += colnorm2(M, j)
+        verbose && println("    j = ", j, ", ||M_{1:b,1:j}||_F^2 = ", FroM2, ", ϵ (|λ_j| + |λ_min|)(|λ_j| + |λ_max|)  = ", tol*maxΛ2j, ".")
+    end
+    j-1
+end
+
+function deflationcheck(M::Matrix, N::Matrix, Λ::Vector, tol::Real, verbose::Bool)
+    b, k = size(M)
+    β, κ = size(N)
+    @assert b == β && k == κ
+    @assert k == length(Λ)
+    maxΛ = maximum(abs, Λ)
+    minΛ = minimum(abs, Λ)
+    j = 1
+    maxΛj = abs2(Λ[j])
+    maxΛ2j = (abs(Λ[j]) + minΛ)*(abs(Λ[j]) + maxΛ)
+    FroMN2 = colnorm2(M, j) + maxΛj*colnorm2(N, j)
+    while FroMN2 ≤ tol*maxΛ2j && j < k
+        j += 1
+        maxΛj = max(maxΛj, abs2(Λ[j]))
+        maxΛ2j = max(maxΛ2j, (abs(Λ[j]) + minΛ)*(abs(Λ[j]) + maxΛ))
+        FroMN2 += colnorm2(M, j) + maxΛj*colnorm2(N, j)
+        verbose && println("    j = ", j, ", ||M_{1:b,1:j}||_F^2 + |λ_j|^2||N_{1:b,1:j}||_F^2 = ", FroMN2, ", ϵ (|λ_j| + |λ_min|)(|λ_j| + |λ_max|)  = ", tol*maxΛ2j, ".")
+    end
+    j-1
+end
+
+function colnorm2(M::Matrix, J::Int)
+    ret = zero(real(eltype(M)))
+    @inbounds for i = 1:size(M, 1)
+        ret += abs2(M[i, J])
+    end
+    ret
+end
diff --git a/src/Operators/Operator.jl b/src/Operators/Operator.jl
index 6abd4a0a4..960eb2054 100644
--- a/src/Operators/Operator.jl
+++ b/src/Operators/Operator.jl
@@ -1,5 +1,5 @@
 export Operator
-export bandwidths, bandrange, \, periodic
+export bandwidth, bandwidths, bandrange, \, periodic
 export neumann
 export ldirichlet,rdirichlet,lneumann,rneumann
 export ldiffbc,rdiffbc,diffbcs
@@ -608,6 +608,7 @@ include("systems.jl")
 include("qr.jl")
 include("nullspace.jl")
 
+include("eigen.jl")
 
 
 
diff --git a/src/Operators/banded/DiagonalOperator.jl b/src/Operators/banded/DiagonalOperator.jl
deleted file mode 100644
index e69de29bb..000000000
diff --git a/src/Operators/eigen.jl b/src/Operators/eigen.jl
new file mode 100644
index 000000000..5ac96c2e8
--- /dev/null
+++ b/src/Operators/eigen.jl
@@ -0,0 +1,294 @@
+const sqrt2m1 = sqrt(2)-1
+
+function symmetric_eigen(L::Operator{T}, n::Int; verbose::Bool = false, tol::Real = eps(T), CONST::Real = sqrt2m1) where T
+    # ndef is the number of deflated eigenvalues
+    # ALLEIGS is a container for the converged eigenvalues
+    # b is the bandwidth of the symmetric banded operator
+    # (k, k) are the dimensions of the deflation window
+    ds = domainspace(L)
+    rs = rangespace(L)
+    ndef = 0
+    ALLEIGS = T[]
+    ALLQ = EmbeddedBlockOperator{T,Matrix{T},typeof(ds),typeof(rs)}[]
+    b = maximum(bandwidths(L))
+    k = 2b
+
+    # We start by determining a good deflation window.
+
+    Lk = Symmetric(L[1:k,1:k])
+    Lb = L[k+1:k+b,1:k]
+    Λ, Q = eigen(Lk)
+    M = Lb*Q
+
+    j = 0
+    while j < 1
+        k *= 2
+        Lk = Symmetric(L[1:k,1:k])
+        Lb = L[k+1:k+b,1:k]
+        Λ, Q = eigen(Lk)
+        M = Lb*Q
+        j = deflationcheck(M, Λ, tol^2, verbose)÷b*b
+        verbose && println("This is the size of the deflation window: ",k," and this is the number of deflated eigenvalues: ",j)
+    end
+
+    SBA = SymBlockArrowHead(SymBandedPlusBulge(Diagonal(Λ[j+1:k]), b, b), Matrix(M[:,j+1:k]'), Matrix(L[k+1:k+b,k+1:k+b]))
+
+    push!(ALLEIGS, Λ[1:j]...)
+    push!(ALLQ, EmbeddedBlockOperator(Q, ndef, ds, rs))
+    ndef += j
+
+    w = zeros(T, k+b)
+    v = zero(w)
+    Lw = zero(w)
+    wQ = zero(w)
+
+    while ndef < n
+        # With the dimensions of the deflation window, (k, k), found & fixed,
+        # we now create a SymBlockArrowHead matrix upon which we are able to
+        # perform similarity transformations to return to banded form.
+
+        Q = Matrix{T}(I, k, k)
+        symblockarrowhead2symbanded!(w, v, Lw, SBA, wQ, Q)
+        push!(ALLQ, EmbeddedBlockOperator(Matrix(Q'), ndef, ds, rs))
+
+        replenishLk!(Lk, SBA, L, ndef)
+        replenishLb!(Lb, L, ndef, k, b)
+
+        if j < CONST*k
+            verbose && println("Failed to deflate! Doubling the dimensions of the deflation window.")
+            Lk = [Lk L[ndef+1:ndef+k,ndef+k+1:ndef+2k]; L[ndef+k+1:ndef+2k,ndef+1:ndef+2k]]
+            k *= 2
+            Lk = SymBandedMatrix(Lk, b, 1:k)
+            Lb = L[ndef+k+1:ndef+k+b, ndef+1:ndef+k]
+            M = zeros(T, b, k)
+            w = zeros(T, k+b)
+            v = zero(w)
+            Lw = zero(w)
+            wQ = zero(w)
+        end
+
+        Λ, Q = eigen(Lk)
+
+        mul!(fill!(M, zero(T)), Lb, Q)
+        j = deflationcheck(M, Λ, tol^2, verbose)÷b*b
+        verbose && println("This is the size of the deflation window: ",k," and this is the number of deflated eigenvalues: ",j)
+
+        SBA = SymBlockArrowHead(SymBandedPlusBulge(Diagonal(Λ[j+1:k]), b, b), Matrix(M[:,j+1:k]'), Matrix(L[ndef+k+1:ndef+k+b,ndef+k+1:ndef+k+b]))
+
+        push!(ALLEIGS, Λ[1:j]...)
+        push!(ALLQ, EmbeddedBlockOperator(Q, ndef, ds, rs))
+        ndef += j
+    end
+
+    ALLEIGS, ALLQ
+end
+
+function symmetric_eigen(L::Operator{T}, C::Operator{T}, n::Int; verbose::Bool = false, tol::Real = eps(T), CONST::Real = sqrt2m1) where T
+    # ndef is the number of deflated eigenvalues
+    # ALLEIGS is a container for the converged eigenvalues
+    # b is the bandwidth of the symmetric banded operator
+    # (k, k) are the dimensions of the deflation window
+    ds = domainspace(L)
+    rs = rangespace(L)
+    ndef = 0
+    ALLEIGS = T[]
+    ALLV = EmbeddedBlockOperator{T,Matrix{T},typeof(ds),typeof(rs)}[]
+    b = max(maximum(bandwidths(L)), maximum(bandwidths(C)))
+    k = 2b
+
+    # We start by determining a good deflation window.
+
+    Lk = Symmetric(L[1:k,1:k])
+    Lb = L[k+1:k+b,1:k]
+    Ck = Symmetric(C[1:k,1:k])
+    Cb = C[k+1:k+b,1:k]
+    Λ, V = eigen(Lk, Ck)
+    M = Lb*V
+    N = Cb*V
+
+    j = 0
+    while j < 2b
+        k *= 2
+        Lk = Symmetric(L[1:k,1:k])
+        Lb = L[k+1:k+b,1:k]
+        Ck = Symmetric(C[1:k,1:k])
+        Cb = C[k+1:k+b,1:k]
+        Λ, V = eigen(Lk, Ck)
+        M = Lb*V
+        N = Cb*V
+        j = deflationcheck(M, N, Λ, tol^2, verbose)÷b*b
+        verbose && println("This is the size of the deflation window: ",k," and this is the number of deflated eigenvalues: ",j)
+    end
+
+    SBAL = SymBlockArrowHead(SymBandedPlusBulge(Diagonal(Λ[j+1:k]), b, b), Matrix(M[:,j+1:k]'), Matrix(L[k+1:k+b,k+1:k+b]))
+    SBAC = SymBlockArrowHead(Diagonal(ones(T, k-j)), Matrix(N[:,j+1:k]'), Matrix(C[k+1:k+b,k+1:k+b]))
+
+    push!(ALLEIGS, Λ[1:j]...)
+    push!(ALLV, EmbeddedBlockOperator(V, ndef, ds, rs))
+    ndef += j
+
+    w = zeros(T, k+b)
+    v = zero(w)
+    Lw = zero(w)
+    wQ = zero(w)
+
+    while ndef < n
+        # With the dimensions of the deflation window, (k, k), found & fixed,
+        # we now create a SymBlockArrowHead matrix upon which we are able to
+        # perform similarity transformations to return to banded form.
+
+        congruence!(SBAL, SBAC)
+        push!(ALLV, EmbeddedBlockOperator(Matrix(UpperTriangular(SBAC)), ndef, ds, rs))
+
+        Q = Matrix{T}(I, k, k)
+        symblockarrowhead2symbanded!(w, v, Lw, SBAL, wQ, Q)
+        push!(ALLV, EmbeddedBlockOperator(Matrix(Q'), ndef, ds, rs))
+
+        replenishLk!(Lk, SBAL, SBAC, L, ndef)
+        replenishLb!(Lb, L, ndef, k, b)
+        replenishCk!(Ck, SBAC, C, b, ndef)
+        replenishLb!(Cb, C, ndef, k, b)
+
+        if j < CONST*k
+            verbose && println("Failed to deflate! Doubling the dimensions of the deflation window.")
+
+            Lk = [Lk L[ndef+1:ndef+k,ndef+k+1:ndef+2k]; L[ndef+k+1:ndef+2k,ndef+1:ndef+2k]]
+            Ck = [Ck C[ndef+1:ndef+k,ndef+k+1:ndef+2k]; C[ndef+k+1:ndef+2k,ndef+1:ndef+2k]]
+            k *= 2
+            Lk = SymBandedMatrix(Lk, b, 1:k)
+            Ck = SymBandedMatrix(Ck, b, 1:k)
+            Lb = L[ndef+k+1:ndef+k+b, ndef+1:ndef+k]
+            Cb = C[ndef+k+1:ndef+k+b, ndef+1:ndef+k]
+            M = zeros(T, b, k)
+            N = zeros(T, b, k)
+            w = zeros(T, k+b)
+            v = zero(w)
+            Lw = zero(w)
+            wQ = zero(w)
+        end
+
+        Λ, V = eigen(Lk, Ck)
+        mul!(fill!(M, zero(T)), Lb, V)
+        mul!(fill!(N, zero(T)), Cb, V)
+        j = deflationcheck(M, N, Λ, tol^2, verbose)÷b*b
+        verbose && println("This is the size of the deflation window: ",k," and this is the number of deflated eigenvalues: ",j)
+
+        SBAL = SymBlockArrowHead(SymBandedPlusBulge(Diagonal(Λ[j+1:k]), b, b), Matrix(M[:,j+1:k]'), Matrix(L[ndef+k+1:ndef+k+b,ndef+k+1:ndef+k+b]))
+        SBAC = SymBlockArrowHead(Diagonal(ones(T, k-j)), Matrix(N[:,j+1:k]'), Matrix(C[ndef+k+1:ndef+k+b,ndef+k+1:ndef+k+b]))
+
+        push!(ALLEIGS, Λ[1:j]...)
+        push!(ALLV, EmbeddedBlockOperator(V, ndef, ds, rs))
+        ndef += j
+    end
+
+    ALLEIGS, ALLV
+end
+
+
+
+symmetric_eigvals(L::Operator{T}, n::Int; kwargs...) where T = symmetric_eigen(L, n; kwargs...)[1][1:n]
+symmetric_eigvals(L::Operator{T}, C::Operator{T}, n::Int; kwargs...) where T = symmetric_eigen(L, C, n; kwargs...)[1][1:n]
+
+
+
+
+
+# Miscellaneous
+
+function SymBandedMatrix(M::AbstractMatrix{T}, b::Int, ir::UnitRange) where T
+    B = BandedMatrix{T}(undef, length(ir), length(ir), b, b)
+    shft = first(ir)-1
+    @inbounds for i in 1:length(ir)
+        for j in colrange(B, i)
+            B[i,j] = M[i+shft, j+shft]
+        end
+    end
+    Symmetric(B)
+end
+
+function replenishLk!(Lk::Symmetric{T,MB}, SBA::SymBlockArrowHead, L::Operator, ndef::Int) where {T,MB<:BandedMatrix{T}}
+    k = size(Lk, 1)
+    l = size(SBA.A, 1)
+    @inbounds for i = 1:l
+        for j in colrange(Lk, i)
+            Lk.data[i,j] = SBA[i,j]
+        end
+    end
+    @inbounds for i = l+1:k
+        for j in colrange(Lk, i)
+            Lk.data[i,j] = L[i+ndef,j+ndef]
+        end
+    end
+    Lk
+end
+
+function replenishLk!(A::Symmetric{T,MB}, SBAL::SymBlockArrowHead, SBAC::SymBlockArrowHead, M::Operator, ndef::Int) where {T,MB<:BandedMatrix{T}}
+    b = SBAL.A.b
+    m, n = size(SBAL.A)
+    for i = size(SBAL.A, 1)+1:size(A, 1)
+        for j in colrange(A, i)
+            A.data[i,j] = M[i+ndef,j+ndef]
+        end
+    end
+    @inbounds for i = size(SBAL, 1)-2b+1:size(SBAL, 1)-b
+        for j in colrange(A, i)
+            A.data[i,j] = SBAL[i,j]
+        end
+    end
+    @inbounds for i = 1:size(SBAL.A, 1)
+        for j in colrange(A, i)
+            A.data[i,j] = SBAL[i,j]
+        end
+    end
+    @inbounds for i = size(SBAL.A, 1)+1:size(SBAL, 1)
+        for j = i:size(SBAL, 2)
+            A.data[i,j] = SBAL[i,j]
+        end
+    end
+    k = size(SBAC.A, 1)
+    temp1 = UpperTriangular(SBAC.C)
+    temp2 = UpperTriangular(Matrix(M[ndef+k+b+1:ndef+k+2b,ndef+k+1:ndef+k+b]))
+    temp3 = temp2*temp1
+    @inbounds for i = k+b+1:k+2b
+        for j = k+1:k+b
+            A.data[i,j] = A.data[j,i] = temp3[i-k-b,j-k]
+        end
+    end
+    A
+end
+
+function replenishCk!(A::Symmetric{T,MB}, SBA::SymBlockArrowHead, M::Operator, b::Int, ndef::Int)  where {T,MB<:BandedMatrix{T}}
+    m, n = size(SBA.A)
+    @inbounds for i = size(SBA.A, 1)+1:size(A, 1)
+        for j in colrange(A, i)
+            A.data[i,j] = M[i+ndef,j+ndef]
+        end
+    end
+    @inbounds for i = 1:size(SBA.A, 1)+b
+        for j in colrange(A, i)
+            A.data[i,j] = 0
+        end
+        A.data[i,i] = 1
+    end
+    k = size(SBA.A, 1)
+    temp1 = UpperTriangular(SBA.C)
+    temp2 = UpperTriangular(Matrix(M[ndef+k+b+1:ndef+k+2b,ndef+k+1:ndef+k+b]))
+    temp3 = temp2*temp1
+    @inbounds for i = k+b+1:k+2b
+        for j = k+1:k+b
+            A.data[i,j] = A.data[j,i] = temp3[i-k-b,j-k]
+        end
+    end
+    A
+end
+
+function replenishLb!(Lb::BandedMatrix{T}, L::Operator{T}, ndef::Int, k::Int, b::Int) where T
+    ndefpk = ndef+k
+    @inbounds for j = k-b+1:k
+        jpndef = j+ndef
+        for i in 1:(b+j-k)
+            Lb[i,j] = L[i+ndefpk, jpndef]
+        end
+    end
+    Lb
+end
diff --git a/src/Operators/general/EmbeddedBlockOperator.jl b/src/Operators/general/EmbeddedBlockOperator.jl
new file mode 100644
index 000000000..12197d28c
--- /dev/null
+++ b/src/Operators/general/EmbeddedBlockOperator.jl
@@ -0,0 +1,134 @@
+export EmbeddedBlockOperator
+
+"""
+    Represent a square block operator embedded in the identity:
+    [ I
+        X
+          I ]
+    where the first identity is k x k and the second is infinite-dimensional.
+"""
+struct EmbeddedBlockOperator{T,M<:AbstractMatrix{T},DS,RS} <: Operator{T}
+    X::M
+    temp::Vector{T}
+    k::Int
+    domainspace::DS
+    rangespace::RS
+    function EmbeddedBlockOperator{T,M,DS,RS}(X::M, k::Int, ds::DS, rs::RS) where {T,M,DS,RS}
+        m, n = size(X)
+        @assert m == n
+        new{T,M,DS,RS}(X, zeros(T, m), k, ds, rs)
+    end
+end
+
+EmbeddedBlockOperator(X::M, k::Int, ds, rs) where M = EmbeddedBlockOperator{eltype(X), M, typeof(ds), typeof(rs)}(X, k, ds, rs)
+EmbeddedBlockOperator(X::M, k::Int) where M = EmbeddedBlockOperator(X, k, SequenceSpace(), SequenceSpace())
+
+bandwidths(A::EmbeddedBlockOperator) = (size(A.X, 1)-1, size(A.X, 2)-1)
+domainspace(A::EmbeddedBlockOperator) = A.domainspace
+rangespace(A::EmbeddedBlockOperator) = A.rangespace
+
+function getindex(A::EmbeddedBlockOperator{T}, i::Integer, j::Integer) where T
+    k = A.k
+    n = size(A.X, 1)
+    if i ≤ k
+        if j ≤ k
+            if i == j
+                return one(T)
+            else
+                return zero(T)
+            end
+        else
+            return zero(T)
+        end
+    elseif k < i ≤ k+n
+        if k < j ≤ k+n
+            return A.X[i-k,j-k]
+        else
+            return zero(T)
+        end
+    else
+        if i == j
+            return one(T)
+        else
+            return zero(T)
+        end
+    end
+end
+
+function lmul!(A::EmbeddedBlockOperator{T}, x::AbstractVector{T}) where T
+    X = A.X
+    y = A.temp
+    k = A.k
+    n = size(X, 1)
+    @assert length(x) ≥ k + n
+    fill!(y, zero(T))
+    @inbounds for j = 1:n
+        xj = x[k+j]
+        for i = 1:n
+            y[i] += X[i,j]*xj
+        end
+    end
+    @inbounds for i = k+1:k+n x[i] = y[i-k] end
+
+    x
+end
+
+function lmul!(A::AdjointOperator{T,M}, x::AbstractVector{T}) where {T,M<:EmbeddedBlockOperator{T}}
+    X = A.op.X
+    y = A.op.temp
+    k = A.op.k
+    n = size(X, 1)
+    @assert length(x) ≥ k + n
+    fill!(y, zero(T))
+    @inbounds for i = 1:n
+        yi = zero(eltype(y))
+        for j = 1:n
+            yi += conj(X[j,i])*x[k+j]
+        end
+        y[i] += yi
+    end
+    @inbounds for i = k+1:k+n x[i] = y[i-k] end
+
+    x
+end
+
+function lmul!(A::TransposeOperator{T,M}, x::AbstractVector{T}) where {T,M<:EmbeddedBlockOperator{T}}
+    X = A.op.X
+    y = A.op.temp
+    k = A.op.k
+    n = size(X, 1)
+    @assert length(x) ≥ k + n
+    fill!(y, zero(T))
+    @inbounds for i = 1:n
+        yi = zero(eltype(y))
+        for j = 1:n
+            yi += X[j,i]*x[k+j]
+        end
+        y[i] += yi
+    end
+    @inbounds for i = k+1:k+n x[i] = y[i-k] end
+
+    x
+end
+
+
+function lmul!(A::Vector{M}, x::Vector{T}) where {T,M<:EmbeddedBlockOperator{T}}
+    for N = length(A):-1:1
+        lmul!(A[N], x)
+    end
+    x
+end
+
+function lmul!(A::Adjoint{AdjointOperator{T,M},Vector{M}}, x::Vector{T}) where {T,M<:EmbeddedBlockOperator{T}}
+    for N = 1:length(A)
+        lmul!(A[N], x)
+    end
+    x
+end
+
+function lmul!(A::Transpose{TransposeOperator{T,M},Vector{M}}, x::Vector{T}) where {T,M<:EmbeddedBlockOperator{T}}
+    for N = 1:length(A)
+        lmul!(A[N], x)
+    end
+    x
+end
diff --git a/src/Operators/general/general.jl b/src/Operators/general/general.jl
index f4924e025..78f7e66f2 100644
--- a/src/Operators/general/general.jl
+++ b/src/Operators/general/general.jl
@@ -2,6 +2,7 @@ include("algebra.jl")
 include("CachedOperator.jl")
 include("FiniteOperator.jl")
 include("OperatorLayout.jl")
+include("EmbeddedBlockOperator.jl")
 include("PartialInverseOperator.jl")
 include("InterlaceOperator.jl")
 include("OperatorFunction.jl")
diff --git a/test/EigTest.jl b/test/EigTest.jl
index e44263ea2..f828022b8 100644
--- a/test/EigTest.jl
+++ b/test/EigTest.jl
@@ -1,5 +1,4 @@
 using ApproxFun, BandedMatrices, SpecialFunctions, Test
-    import ApproxFun: bandwidth
 
 #
 # Fix for BigFloat stackoverflow