ArrowheadMatrix -> BBBArrowheadMatrix (#61)

Author: dlfivefifty

src/arrowhead.jl

@@ -1,52 +1,52 @@
 
 """
-ArrowheadMatrix
+BBBArrowheadMatrix
 
     A   B   B
     C   D   D  D
     …   …   …   …   …
     C   D   D
 """
-struct ArrowheadMatrix{T, AA<:AbstractMatrix{T}, BB, CC, DD} <: AbstractBandedBlockBandedMatrix{T}
+struct BBBArrowheadMatrix{T, AA<:AbstractMatrix{T}, BB, CC, DD} <: AbstractBandedBlockBandedMatrix{T}
     A::AA
     B::BB # first row blocks
     C::CC # first col blocks
     D::DD # these are interlaces
 
-    ArrowheadMatrix{T, AA, BB, CC, DD}(A, B, C, D) where {T,AA,BB,CC,DD} = new{T,AA,BB,CC,DD}(A, B, C, D)
+    BBBArrowheadMatrix{T, AA, BB, CC, DD}(A, B, C, D) where {T,AA,BB,CC,DD} = new{T,AA,BB,CC,DD}(A, B, C, D)
 end
 
-ArrowheadMatrix{T}(A::AbstractMatrix{T}, B, C, D) where T = ArrowheadMatrix{T, typeof(A), typeof(B), typeof(C), typeof(D)}(A, B, C, D)
-ArrowheadMatrix{T}(A, B, C, D) where T = ArrowheadMatrix(convert(AbstractMatrix{T}, A), B, C, D)
+BBBArrowheadMatrix{T}(A::AbstractMatrix{T}, B, C, D) where T = BBBArrowheadMatrix{T, typeof(A), typeof(B), typeof(C), typeof(D)}(A, B, C, D)
+BBBArrowheadMatrix{T}(A, B, C, D) where T = BBBArrowheadMatrix(convert(AbstractMatrix{T}, A), B, C, D)
 
-const ArrowheadMatrices = Union{ArrowheadMatrix,Symmetric{<:Any,<:ArrowheadMatrix},Hermitian{<:Any,<:ArrowheadMatrix},
-                                UpperOrLowerTriangular{<:Any,<:ArrowheadMatrix}}
+const ArrowheadMatrices = Union{BBBArrowheadMatrix,Symmetric{<:Any,<:BBBArrowheadMatrix},Hermitian{<:Any,<:BBBArrowheadMatrix},
+                                UpperOrLowerTriangular{<:Any,<:BBBArrowheadMatrix}}
 
 
 
 subblockbandwidths(A::ArrowheadMatrices) = (1,1)
 
-BlockArrays._show_typeof(io::IO, B::ArrowheadMatrix{T}) where T = print(io, "ArrowheadMatrix{$T}")
+BlockArrays._show_typeof(io::IO, B::BBBArrowheadMatrix{T}) where T = print(io, "BBBArrowheadMatrix{$T}")
 
-function blockbandwidths(A::ArrowheadMatrix)
+function blockbandwidths(A::BBBArrowheadMatrix)
     l,u = bandwidths(A.D[1])
     max(l,length(A.C)),max(u,length(A.B))
 end
 
-function axes(L::ArrowheadMatrix)
+function axes(L::BBBArrowheadMatrix)
     ξ,n = size(L.A)
     m = length(L.D)
     μ,ν = size(L.D[1])
     blockedrange(Vcat(ξ, Fill(m,μ))), blockedrange(Vcat(n, Fill(m,ν)))
 end
 
-copy(A::ArrowheadMatrix) = ArrowheadMatrix(copy(A.A), map(copy, A.B), map(copy, A.C), map(copy, A.D))
+copy(A::BBBArrowheadMatrix) = BBBArrowheadMatrix(copy(A.A), map(copy, A.B), map(copy, A.C), map(copy, A.D))
 
 for adj in (:adjoint, :transpose)
-    @eval $adj(A::ArrowheadMatrix) = ArrowheadMatrix($adj(A.A), map($adj, A.C), map($adj, A.B), map($adj, A.D))
+    @eval $adj(A::BBBArrowheadMatrix) = BBBArrowheadMatrix($adj(A.A), map($adj, A.C), map($adj, A.B), map($adj, A.D))
 end
 
-function getindex(L::ArrowheadMatrix{T}, Kk::BlockIndex{1}, Jj::BlockIndex{1})::T where T
+function getindex(L::BBBArrowheadMatrix{T}, Kk::BlockIndex{1}, Jj::BlockIndex{1})::T where T
     K,k = block(Kk),blockindex(Kk)
     J,j = block(Jj),blockindex(Jj)
     J == K == Block(1) && return L.A[k,j]
@@ -60,12 +60,12 @@ function getindex(L::ArrowheadMatrix{T}, Kk::BlockIndex{1}, Jj::BlockIndex{1})::
     return L.D[k][Int(K)-1, Int(J)-1]
 end
 
-function getindex(L::ArrowheadMatrix, k::Int, j::Int)
+function getindex(L::BBBArrowheadMatrix, k::Int, j::Int)
     ax,bx = axes(L)
     L[findblockindex(ax, k), findblockindex(bx, j)]
 end
 
-function ArrowheadMatrix(A, B, C, D)
+function BBBArrowheadMatrix(A, B, C, D)
     ξ,n = size(A)
     m = length(D)
     μ,ν = size(D[1])
@@ -95,7 +95,7 @@ function ArrowheadMatrix(A, B, C, D)
     end
     T = promote_type(eltype(A), mapreduce(eltype, promote_type, B; init=eltype(A)),
                      mapreduce(eltype, promote_type, C; init=eltype(A)), mapreduce(eltype, promote_type, D; init=eltype(A)))
-    ArrowheadMatrix{T}(A, B, C, D)
+    BBBArrowheadMatrix{T}(A, B, C, D)
 end
 
 
@@ -113,24 +113,24 @@ arrowheadlayout(::BandedLazyLayouts) = LazyArrowheadLayout()
 arrowheadlayout(::DiagonalLayout{<:AbstractLazyLayout}) = LazyArrowheadLayout()
 symmetriclayout(lay::ArrowheadLayouts) = SymmetricLayout{typeof(lay)}()
 
-MemoryLayout(::Type{<:ArrowheadMatrix{<:Any,<:Any,<:Any,<:Any,<:AbstractVector{D}}}) where D = arrowheadlayout(MemoryLayout(D))
+MemoryLayout(::Type{<:BBBArrowheadMatrix{<:Any,<:Any,<:Any,<:Any,<:AbstractVector{D}}}) where D = arrowheadlayout(MemoryLayout(D))
 # Hack since ∞ diagonal fill is only DiagonalLayout{FillLayout}
-MemoryLayout(::Type{<:ArrowheadMatrix{<:Any,<:Any,<:Any,<:Any,<:AbstractVector{<:Diagonal{<:Any,<:AbstractFill{<:Any,1,<:Tuple{OneToInf}}}}}}) = LazyArrowheadLayout()
+MemoryLayout(::Type{<:BBBArrowheadMatrix{<:Any,<:Any,<:Any,<:Any,<:AbstractVector{<:Diagonal{<:Any,<:AbstractFill{<:Any,1,<:Tuple{OneToInf}}}}}}) = LazyArrowheadLayout()
 
 sublayout(::ArrowheadLayouts,
           ::Type{<:NTuple{2,BlockSlice{<:BlockRange{1, Tuple{OneTo{Int}}}}}}) = ArrowheadLayout()
 function sub_materialize(::ArrowheadLayout, V::AbstractMatrix)
     KR,JR = parentindices(V)
     P = parent(V)
     M,N =  KR.block[end],JR.block[end]
-    ArrowheadMatrix(P.A, P.B, P.C,
+    BBBArrowheadMatrix(P.A, P.B, P.C,
                     layout_getindex.(P.D, Ref(oneto(Int(M)-1)), Ref(oneto(Int(N)-1))))
 end
 
 symmetric(A) = Symmetric(A)
 symmetric(A::Union{SymTridiagonal,Symmetric,Diagonal}) = A
 
-function getproperty(F::Symmetric{<:Any,<:ArrowheadMatrix}, d::Symbol)
+function getproperty(F::Symmetric{<:Any,<:BBBArrowheadMatrix}, d::Symbol)
     P = getfield(F, :data)
     if d == :A
         return symmetric(P.A)
@@ -242,19 +242,19 @@ end
 # Cholesky
 ####
 
-function reversecholcopy(S::Symmetric{<:Any,<:ArrowheadMatrix})
+function reversecholcopy(S::Symmetric{<:Any,<:BBBArrowheadMatrix})
     T = LinearAlgebra.choltype(S)
     A = parent(S)
-    Symmetric(ArrowheadMatrix(LinearAlgebra.copymutable_oftype(A.A, T),
+    Symmetric(BBBArrowheadMatrix(LinearAlgebra.copymutable_oftype(A.A, T),
     LinearAlgebra.copymutable_oftype.(A.B, T), LinearAlgebra.copymutable_oftype.(A.C, T),
     LinearAlgebra.copymutable_oftype.(A.D, T)))
 end
 
 
 
-function MatrixFactorizations._reverse_chol!(A::ArrowheadMatrix, ::Type{UpperTriangular})
-
-    for B in A.D
+function MatrixFactorizations._reverse_chol!(A::BBBArrowheadMatrix, ::Type{UpperTriangular})
+    
+    Threads.@threads for B in A.D
         reversecholesky!(Symmetric(B))
     end
 
@@ -384,16 +384,16 @@ tupleop(op, A::Tuple, B::Tuple) = (op(first(A), first(B)), tupleop(op, tail(A),
 
 
 for op in (:+, :-)
-    @eval $op(A::ArrowheadMatrix, B::ArrowheadMatrix) = ArrowheadMatrix($op(A.A, B.A), tupleop($op, A.B, B.B), tupleop($op, A.C, B.C), $op(A.D, B.D))
+    @eval $op(A::BBBArrowheadMatrix, B::BBBArrowheadMatrix) = BBBArrowheadMatrix($op(A.A, B.A), tupleop($op, A.B, B.B), tupleop($op, A.C, B.C), $op(A.D, B.D))
 end
--(A::ArrowheadMatrix) = ArrowheadMatrix(-A.A, map(-, A.B), map(-, A.C), -A.D)
+-(A::BBBArrowheadMatrix) = BBBArrowheadMatrix(-A.A, map(-, A.B), map(-, A.C), -A.D)
 
 for op in (:*, :\)
-    @eval $op(c::Number, A::ArrowheadMatrix) = ArrowheadMatrix($op(c, A.A), broadcast($op, c, A.B), broadcast($op, c, A.C), broadcast($op, c, A.D))
+    @eval $op(c::Number, A::BBBArrowheadMatrix) = BBBArrowheadMatrix($op(c, A.A), broadcast($op, c, A.B), broadcast($op, c, A.C), broadcast($op, c, A.D))
 end
 
 for op in (:*, :/)
-    @eval $op(A::ArrowheadMatrix, c::Number) = ArrowheadMatrix($op(A.A, c), broadcast($op, A.B, c), broadcast($op, A.C, c), broadcast($op, A.D, c))
+    @eval $op(A::BBBArrowheadMatrix, c::Number) = BBBArrowheadMatrix($op(A.A, c), broadcast($op, A.B, c), broadcast($op, A.C, c), broadcast($op, A.D, c))
 end
 
 ###
@@ -416,7 +416,7 @@ for (UNIT, Tri) in (('U',UnitUpperTriangular), ('N', UpperTriangular))
         A,B,D = P.A,P.B,P.D
         m = length(D)
 
-        for k = 1:m
+        Threads.@threads for k = 1:m
             ArrayLayouts.ldiv!($Tri(D[k]), view(dest, n+k:m:length(dest)))
         end
 
@@ -467,7 +467,7 @@ end
 
 
 for Tri in (:UpperTriangular, :UnitUpperTriangular)
-    @eval function getproperty(F::$Tri{<:Any,<:ArrowheadMatrix}, d::Symbol)
+    @eval function getproperty(F::$Tri{<:Any,<:BBBArrowheadMatrix}, d::Symbol)
         P = getfield(F, :data)
         if d == :A
             return $Tri(P.A)
@@ -484,7 +484,7 @@ for Tri in (:UpperTriangular, :UnitUpperTriangular)
 end
 
 for Tri in (:LowerTriangular, :UnitLowerTriangular)
-    @eval function getproperty(F::$Tri{<:Any,<:ArrowheadMatrix}, d::Symbol)
+    @eval function getproperty(F::$Tri{<:Any,<:BBBArrowheadMatrix}, d::Symbol)
         P = getfield(F, :data)
         if d == :A
             return $Tri(P.A)

src/continuouspolynomial.jl

@@ -228,7 +228,7 @@ function \(P::ContinuousPolynomial{0}, C::ContinuousPolynomial{1})
     @assert P.points == C.points
     v = one(T) ./ (3:2:∞)
     N = length(P.points)
-    ArrowheadMatrix(_BandedMatrix(Ones{T}(2, N)/2, oneto(N-1), 0, 1),
+    BBBArrowheadMatrix(_BandedMatrix(Ones{T}(2, N)/2, oneto(N-1), 0, 1),
         (_BandedMatrix(Fill(v[1], 1, N-1), oneto(N-1), 0, 0),),
         (_BandedMatrix(Vcat(Ones{T}(1, N)/2, -Ones{T}(1, N)/2), oneto(N-1), 0, 1),),
         Fill(_BandedMatrix(Hcat(v, Zeros{T}(∞), -v)', axes(v,1), 1, 1), N-1))
@@ -243,7 +243,7 @@ end
 
     # The basis for D is defined in each element a..b as diff(legendre(a..b)) with delta functions. But note that
     # the delta functions don't change! We need to kill off the deltas in conversion.
-    ArrowheadMatrix(_BandedMatrix(Vcat(Fill(-one(T),1,N-1), Fill(one(T),1,N-1)), N-2, 0, 1),
+    BBBArrowheadMatrix(_BandedMatrix(Vcat(Fill(-one(T),1,N-1), Fill(one(T),1,N-1)), N-2, 0, 1),
         (_BandedMatrix(Ones{T}(2,N-1), N-2, 0, 1),),
         (SquareEye{T}(N-1),),
         Fill(R[:,2:end], N-1))
@@ -261,7 +261,7 @@ function grammatrix(A::ContinuousPolynomial{0,T}) where T
     N = length(r) - 1
     hs = diff(r)
     M = grammatrix(Legendre{T}())
-    ArrowheadMatrix{T}(Diagonal(hs), (), (), [Diagonal(M.diag[2:end] * h/2) for h in hs])
+    BBBArrowheadMatrix{T}(Diagonal(hs), (), (), [Diagonal(M.diag[2:end] * h/2) for h in hs])
 end
 
 function grammatrix(A::ContinuousPolynomial{0,T, <:AbstractRange}) where T
@@ -283,7 +283,7 @@ function grammatrix(C::ContinuousPolynomial{1, T, <:AbstractRange}) where T
     a21 = _BandedMatrix(Fill(h/6, 2, N), N+1, 1, 0)
     a31 = _BandedMatrix(Vcat(Fill(-h/30, 1, N), Fill(h/30, 1, N)), N+1, 1, 0)
 
-    Symmetric(ArrowheadMatrix(a11, (a21, a31), (),
+    Symmetric(BBBArrowheadMatrix(a11, (a21, a31), (),
                 Fill(_BandedMatrix(Vcat((h*a/2)',
                 Zeros{T}(1,∞),
                 (h*b/2)'), ∞, 0, 2), N)))
@@ -320,7 +320,7 @@ function diff(C::ContinuousPolynomial{1,T}; dims=1) where T
     r = C.points
     N = length(r)
     s = one(T) ./ (r[2:end]-r[1:end-1])
-    ContinuousPolynomial{0}(r) * ArrowheadMatrix(_BandedMatrix(Vcat([0; s]', [-s; 0]'), length(s), 0, 1), (), (),
+    ContinuousPolynomial{0}(r) * BBBArrowheadMatrix(_BandedMatrix(Vcat([0; s]', [-s; 0]'), length(s), 0, 1), (), (),
                                                  (-2s) .* Ref(Eye{T}(∞)))
 end
 
@@ -329,7 +329,7 @@ function diff(C::ContinuousPolynomial{1,T,<:AbstractRange}; dims=1) where T
     r = C.points
     N = length(r)
     s = 1 ./ step(r)
-    ContinuousPolynomial{0}(r) * ArrowheadMatrix(_BandedMatrix(Vcat(Fill(s, 1, N), Fill(-s, 1, N)), N-1, 0, 1), (), (),
+    ContinuousPolynomial{0}(r) * BBBArrowheadMatrix(_BandedMatrix(Vcat(Fill(s, 1, N), Fill(-s, 1, N)), N-1, 0, 1), (), (),
                                                  Fill(-2s*Eye{T}(∞), N-1))
 end
 
@@ -338,7 +338,7 @@ function diff(P::ContinuousPolynomial{0,T,<:AbstractRange}; dims=1) where T
     N = length(r)
     s = step(r)
 
-    ContinuousPolynomial{-1}(r) * ArrowheadMatrix(_BandedMatrix(Vcat(Fill(one(T),1,N-1), Fill(-one(T),1,N-1)), N-2, 0, 1),
+    ContinuousPolynomial{-1}(r) * BBBArrowheadMatrix(_BandedMatrix(Vcat(Fill(one(T),1,N-1), Fill(-one(T),1,N-1)), N-2, 0, 1),
         (),
         (),
         Fill(Eye{T}(∞) * (2/s), N-1))
@@ -352,7 +352,7 @@ function weaklaplacian(C::ContinuousPolynomial{1,T,<:AbstractRange}) where T
     si = inv(s)
     t1 = Vcat(-si, Fill(-2si, N-2), -si)
     t2 = Fill(si, N-1)
-    Symmetric(ArrowheadMatrix(LazyBandedMatrices.Bidiagonal(t1, t2, :U), (), (),
+    Symmetric(BBBArrowheadMatrix(LazyBandedMatrices.Bidiagonal(t1, t2, :U), (), (),
         Fill(Diagonal(convert(T,-4) ./ (s*(convert(T,3):2:∞))), N-1)))
 end
 

src/dirichlet.jl

@@ -36,15 +36,15 @@ function \(C::ContinuousPolynomial{1}, Q::DirichletPolynomial)
     T = promote_type(eltype(Q), eltype(C))
     @assert C.points == Q.points
     m = length(C.points)
-    ArrowheadMatrix(layout_getindex(Eye{T}(m), :, 2:m-1), (), (), Fill(Eye{T}(∞), m-1))
+    BBBArrowheadMatrix(layout_getindex(Eye{T}(m), :, 2:m-1), (), (), Fill(Eye{T}(∞), m-1))
 end
 
 # \ understood as pseudo invers
 function \(Q::DirichletPolynomial, C::ContinuousPolynomial{1})
     T = promote_type(eltype(Q), eltype(C))
     @assert C.points == Q.points
     m = length(C.points)
-    ArrowheadMatrix(layout_getindex(Eye{T}(m), 2:m-1, :), (), (), Fill(Eye{T}(∞), m-1))
+    BBBArrowheadMatrix(layout_getindex(Eye{T}(m), 2:m-1, :), (), (), Fill(Eye{T}(∞), m-1))
 end
 
 function \(P::ContinuousPolynomial, Q::DirichletPolynomial)
@@ -173,7 +173,7 @@ function weaklaplacian(C::DirichletPolynomial{T,<:AbstractRange}) where T
     si = inv(s)
     t1 = Fill(-2si, N-2)
     t2 = Fill(si, N-3)
-    Symmetric(ArrowheadMatrix(LazyBandedMatrices.Bidiagonal(t1, t2, :U), (), (),
+    Symmetric(BBBArrowheadMatrix(LazyBandedMatrices.Bidiagonal(t1, t2, :U), (), (),
         Fill(Diagonal(convert(T, -4) ./ (s*(convert(T,3):2:∞))), N-1)))
 end
 
@@ -189,7 +189,7 @@ function grammatrix(Q::DirichletPolynomial{T, <:AbstractRange}) where T
     a21 = _BandedMatrix(Fill(h/6, 2, N), N-1, 0, 1)
     a31 = _BandedMatrix(Vcat(Fill(-h/30, 1, N), Fill(h/30, 1, N)), N-1, 0, 1)
 
-    Symmetric(ArrowheadMatrix(a11, (a21, a31), (),
+    Symmetric(BBBArrowheadMatrix(a11, (a21, a31), (),
                 Fill(_BandedMatrix(Vcat((h*a/2)',
                 Zeros{T}(1,∞),
                 (h*b/2)'), ∞, 0, 2), N)))