Three-term ldiv for UpperTriangular AlmostBandedMatrix (#615)

* Three-term ldiv for UpperTriangular AlmostBandedMatrix

* Tests for ldiv

* Fix import

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunBase"
 uuid = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
-version = "0.9.19"
+version = "0.9.20"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/Caching/almostbanded.jl

@@ -392,4 +392,10 @@ for ArrTyp in (:AbstractVector, :AbstractMatrix)
         end
         u
     end
+    @eval function ldiv!(v::$ArrTyp{T},
+                U::UpperTriangular{T,<:SubArray{T, 2, <:AlmostBandedMatrix{T}, NTuple{2,UnitRange{Int}}}},
+                u::$ArrTyp{T}) where T
+
+        ldiv!(U, copyto!(v, u))
+    end
 end

src/LinearAlgebra/AlmostBandedMatrix.jl

@@ -35,7 +35,7 @@ end
 size(A::AlmostBandedMatrix) = size(A.bands)
 
 
-function getindex(B::AlmostBandedMatrix,k::Integer,j::Integer)
+Base.@propagate_inbounds function getindex(B::AlmostBandedMatrix,k::Integer,j::Integer)
     if j > k + bandwidth(B.bands,2)
         B.fill[k,j]
     else
@@ -44,7 +44,7 @@ function getindex(B::AlmostBandedMatrix,k::Integer,j::Integer)
 end
 
 # can only change the bands, not the fill
-function setindex!(B::AlmostBandedMatrix,v,k::Integer,j::Integer)
+Base.@propagate_inbounds function setindex!(B::AlmostBandedMatrix,v,k::Integer,j::Integer)
         B.bands[k,j] = v
 end
 

test/runtests.jl

@@ -11,6 +11,7 @@ using FillArrays
 using InfiniteArrays
 using Infinities
 using LinearAlgebra
+using LowRankMatrices
 using Random
 using SpecialFunctions
 using Test
@@ -249,6 +250,11 @@ end
     AInt = convert(AbstractArray{Int}, A)
     @test AInt isa AbstractArray{Int}
     @test AInt == A
+
+    bands, fill = BandedMatrix(0=>Float64[1:4;]), LowRankMatrix(Float64[1:4;], Float64[1:4;])
+    AB = ApproxFunBase.AlmostBandedMatrix(bands, fill)
+    U = UpperTriangular(view(AB, 1:4, 1:4))
+    @test ldiv!(U, Float64[1:4;]) == ldiv!(ones(4), U, Float64[1:4;]) == ldiv!(factorize(Array(U)), Float64[1:4;])
 end
 
 @testset "DiracDelta sampling" begin