sketch new type

Author: dlfivefifty

src/semiseparable_arrowhead.jl

@@ -0,0 +1,39 @@
+struct SemiseparableBBBArrowheadMatrix{T} <: AbstractBlockBandedMatrix{T}
+    # banded parts
+    A::BandedMatrix{T}
+    B::NTuple{2,BandedMatrix{T}} # first row blocks
+    C::NTuple{4,BandedMatrix{T}} # first col blocks
+    D
+
+    # fill parts
+    Asub::NTuple{2,Vector{T}}
+    Asup::Tuple{2,Matrix{T}} # matrices are m × 2
+
+    Bsub::NTuple{2,Vector{T}}
+    Bsup::NTuple{2,NTuple{2,Vector{T}}}
+
+    Csub::NTuple{2,NTuple{2,Vector{T}}}
+    Csup::NTuple{2,Vector{T}}
+
+    A22sub::NTuple{2,Vector{T}}
+    A32sub::NTuple{2,Vector{T}}
+
+    A32extra::Vector{T}
+    A33extra::Vector{T}
+
+    D::DD # these are interlaces
+
+end
+
+function getindex(L::SemiseparableBBBArrowheadMatrix{T}, Kk::BlockIndex{1}, Jj::BlockIndex{1})::T where T
+    K,k = block(Kk),blockindex(Kk)
+    J,j = block(Jj),blockindex(Jj)
+    # TODO: add getindex
+end
+
+
+function getindex(L::SemiseparableBBBArrowheadMatrix, k::Int, j::Int)
+    ax,bx = axes(L)
+    L[findblockindex(ax, k), findblockindex(bx, j)]
+end
+

test/explore_QL.jl

@@ -1,4 +1,4 @@
-using PiecewiseOrthogonalPolynomials, Plots, BlockArrays
+using PiecewiseOrthogonalPolynomials, Plots, BlockArrays, Test
 using MatrixFactorizations, LinearAlgebra, BlockBandedMatrices
 ###
 # QL
@@ -9,9 +9,9 @@ function my_ql(A::BBBArrowheadMatrix{T}) where T
     m2, n2 = size(A.D[1])
     @assert m == n == l+1
     @assert m2 == n2
-    #results stored in F and tau
+    #results stored in F and τ
     F = BlockedArray(Matrix(A), axes(A))
-    tau = zeros(m+l*m2)
+    τ = zeros(m+l*m2)
     for j in m2:-1:3
         for i in l:-1:1
             upper_entry = F[Block(j-1, j+1)][i, i] #A.D[i][j-2,j]
@@ -28,7 +28,7 @@ function my_ql(A::BBBArrowheadMatrix{T}) where T
             print(dia_entry_new)
             F[m+(j-1)*l+i, m+(j-1)*l+i] = dia_entry_new #update F[Block(j+1, j+1)][i, i]
             F[m+(j-3)*l+i, m+(j-1)*l+i] = v[1]/v[2] #update F[Block(j-1, j+1)][i, i]
-            tau[m+(j-1)*l+i] = coef*v[2]^2
+            τ[m+(j-1)*l+i] = coef*v[2]^2
             #row recombination(householder transformation) for other columns
             current_upper_entry = F[Block(j-1, j-1)][i, i] #A.D[i][j-2,j-2]
             current_lower_entry = F[Block(j+1, j-1)][i, i] #A.D[i][j,j-2]
@@ -53,41 +53,41 @@ function my_ql(A::BBBArrowheadMatrix{T}) where T
     end
 
     #Deal with Block(1,3)
-    #vectors x and Lambda denote a rank 1 semiseperable matrix
-    lambda = 1.0
-    Lambda = []
+    #vectors x and Λ denote a rank 1 semiseperable matrix
+    λ = 1.0
+    Λ = []
     x = [F[Block(1,3)][l+1,l]]
     x_len = abs(x[1])
     for i in l:-1:2 #consider i=1 later
         a = F[Block(1,3)][i,i]
         b = F[Block(1,3)][i,i-1]
         c = F[Block(3,3)][i,i]
-        v_last = c + sign(c) * sqrt(a^2 + lambda^2 * x_len^2 + c^2)
-        v_len = sqrt(a^2 + lambda^2 * x_len^2 + v_last^2)
-        F[m+l+i,m+l+i] = -sign(c) * sqrt(a^2 + lambda^2 * x_len^2 + c^2)
-        pushfirst!(Lambda, lambda / v_last)
-        lambda = -2/v_len^2 * a * b * lambda
+        v_last = c + sign(c) * sqrt(a^2 + λ^2 * x_len^2 + c^2)
+        v_len = sqrt(a^2 + λ^2 * x_len^2 + v_last^2)
+        F[m+l+i,m+l+i] = -sign(c) * sqrt(a^2 + λ^2 * x_len^2 + c^2)
+        pushfirst!(Λ, λ / v_last)
+        λ = -2/v_len^2 * a * b * λ
         F[m+l+i, m+l+i-1] = -2/v_len^2 * v_last * a * b
-        x_first = (1 - 2/v_len^2 * a^2) * b / lambda
+        x_first = (1 - 2/v_len^2 * a^2) * b / λ
         pushfirst!(x, x_first)  
         x_len = sqrt(x_len^2 + x_first^2)  
         #record information of V
         F[i+1, m+l+i] = 0
         F[i, m+l+i] = a / v_last
-        tau[m+l+i] = 2 * v_last^2 / v_len^2
+        τ[m+l+i] = 2 * v_last^2 / v_len^2
     end
     #deal with the last column in Block(1,3)
     a = F[Block(1,3)][1,1]
     c = F[Block(3,3)][1,1]
-    v_last = c + sign(c) * sqrt(a^2 + lambda^2 * x_len^2 + c^2)
-    v_len = sqrt(a^2 + lambda^2 * x_len^2 + v_last^2)
-    pushfirst!(Lambda, lambda / v_last)
-    F[m+l+1,m+l+1] = -sign(c) * sqrt(a^2 + lambda^2 * x_len^2 + c^2)
+    v_last = c + sign(c) * sqrt(a^2 + λ^2 * x_len^2 + c^2)
+    v_len = sqrt(a^2 + λ^2 * x_len^2 + v_last^2)
+    pushfirst!(Λ, λ / v_last)
+    F[m+l+1,m+l+1] = -sign(c) * sqrt(a^2 + λ^2 * x_len^2 + c^2)
     F[2, m+l+1] = 0
     F[1, m+l+1] = a / v_last
-    tau[m+l+1] = 2 * v_last^2 / v_len^2
+    τ[m+l+1] = 2 * v_last^2 / v_len^2
 
-    F, tau, x, Lambda
+    F, τ, x, Λ
 end
 
 
@@ -104,6 +104,11 @@ KR = Block.(Base.OneTo(N))
 Mₙ = M[KR,KR]
 Δₙ = Δ[KR,KR]
 A = Δₙ + 100^2 * Mₙ
-FF,ttau, xx, LLambda = my_ql(A)
-#tau = ql(A).τ
-#f = ql(A).factors
+FF,tτ, xx, LΛ = my_ql(A)
+τ = ql(A).τ
+f = ql(A).factors
+
+@test BlockedArray(tτ, (axes(A,2),))[Block.(3:6)] ≈ BlockedArray(τ, (axes(A,2),))[Block.(3:6)]
+@test f[:,Block.(4:6)] ≈ FF[:,Block.(4:6)]
+
+@test tril(xx * LΛ') + FF[Block(1,3)][2:end,:] ≈ f[Block(1,3)][2:end,:]