Feat symmetric eigen

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

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

src/LinearAlgebra/LinearAlgebra.jl

@@ -19,3 +19,7 @@ include("RaggedMatrix.jl")
 
 
 include("rowvector.jl")
+
+include("deflationcheck.jl")
+include("SymBandedPlusBulge.jl")
+include("SymBlockArrowHead.jl")

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