bidiag_forwardsub! with padded (#354)

* bidiag_forwardsub! with padded

* avoid method overload error downstream

* Update paddedtests.jl

* increase cov

* Update paddedtests.jl

Project.toml

@@ -1,6 +1,6 @@
 name = "LazyArrays"
 uuid = "5078a376-72f3-5289-bfd5-ec5146d43c02"
-version = "2.2.4"
+version = "2.2.5"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/LazyArrays.jl

@@ -38,7 +38,8 @@ import ArrayLayouts: AbstractQLayout, Dot, Dotu, Ldiv, Lmul, MatMulMatAdd, MatMu
                      hermitianlayout, layout_getindex, layout_replace_in_print_matrix, ldivaxes, materialize,
                      materialize!, mulreduce, reshapedlayout, rowsupport, scalarone, scalarzero, sub_materialize,
                      sublayout, symmetriclayout, symtridiagonallayout, transposelayout, triangulardata,
-                     triangularlayout, tridiagonallayout, zero!, transtype, OnesLayout
+                     triangularlayout, tridiagonallayout, zero!, transtype, OnesLayout,
+                     diagonaldata, subdiagonaldata, supdiagonaldata
 
 import FillArrays: AbstractFill, getindex_value
 

src/padded.jl

@@ -589,3 +589,35 @@ PaddedArray(A::T, n::Vararg{Integer,N}) where {T<:Number,N} = ApplyArray{T,N}(se
 
 
 BroadcastStyle(::Type{<:PaddedArray{<:Any,N}}) where N = LazyArrayStyle{N}()
+
+
+
+function ArrayLayouts._bidiag_forwardsub!(M::Ldiv{<:Any,<:PaddedColumns,<:AbstractMatrix,<:AbstractVector})
+    A, b_in = M.A, M.B
+    dv = diagonaldata(A)
+    ev = subdiagonaldata(A)
+    b = paddeddata(b_in)
+    N = length(b)
+    dvj = dv[1]
+    iszero(dvj) && throw(SingularException(1))
+    b[1] = bj1 = dvj\b[1]
+    @inbounds for j = 2:N
+        bj  = b[j]
+        bj -= ev[j - 1] * bj1
+        dvj = dv[j]
+        iszero(dvj) && throw(SingularException(j))
+        bj   = dvj\bj
+        b[j] = bj1 = bj
+    end
+
+    @inbounds for j = N+1:length(b_in)
+        iszero(bj1) && break
+        bj = -ev[j - 1] * bj1
+        dvj = dv[j]
+        iszero(dvj) && throw(SingularException(j))
+        bj   = dvj\bj
+        b_in[j] = bj1 = bj
+    end
+
+    b_in
+end

test/paddedtests.jl

@@ -420,6 +420,18 @@ paddeddata(a::PaddedPadded) = a
         @test C'b ≈ Matrix(C)'b
         @test b'C ≈ b'Matrix(C)
     end
-end
 
+    @testset "Bidiagonal" begin
+        B = Bidiagonal(1:5, 1:4, :L)
+        b = Vcat(randn(5), Zeros(0))
+        @test ArrayLayouts.ldiv!(B, deepcopy(b)) ≈ B\b
+        c = cache(Zeros(5)); c[1] = 2;
+        @test ArrayLayouts.ldiv!(B, c) ≈ B\[2; zeros(4)]
+
+        c = cache(Zeros(5)); c[1:2] = [1,2];
+        @test_throws SingularException ArrayLayouts.ldiv!(Bidiagonal(0:4, 1:4, :L), c)
+        @test_throws SingularException ArrayLayouts.ldiv!(Bidiagonal(-1:3, 1:4, :L), c)
+        @test_throws SingularException ArrayLayouts.ldiv!(Bidiagonal(-4:0, 1:4, :L), c)
+    end
+end
 end # module