Fix lower Cholesky of banded matrices with upper bandwidth > 0

src/symbanded/BandedCholesky.jl

@@ -1,4 +1,11 @@
 
+# LAPACK's lower band storage starts at the diagonal, which `bandeddata(A)` stores
+# in row `bandwidth(A, 2) + 1`, below the superdiagonals.
+function _lower_bandeddata(A)
+    data = bandeddata(A)
+    view(data, bandwidth(A, 2) + 1:size(data, 1), :)
+end
+
 function banded_chol!(::BandedColumns{DenseColumnMajor},
                        A::AbstractMatrix{T}, ::Type{UpperTriangular}) where T<:BlasFloat
     _, info = pbtrf!('U', size(A, 1), bandwidth(A,2), bandeddata(A))
@@ -7,7 +14,7 @@ end
 
 function banded_chol!(::BandedColumns{DenseColumnMajor},
                       A::AbstractMatrix{T}, ::Type{LowerTriangular}) where T<:BlasFloat
-    _, info = pbtrf!('L', size(A, 1), bandwidth(A,1), bandeddata(A))
+    _, info = pbtrf!('L', size(A, 1), bandwidth(A,1), _lower_bandeddata(A))
     LowerTriangular(A), info
 end
 
@@ -74,7 +81,7 @@ function _ldiv!(::BandedColumns{DenseColumnMajor}, ::AbstractStridedLayout,
     if A.uplo == 'U'
         pbtrs!('U', size(A, 1), bandwidth(A.factors,2), bandeddata(A.factors), B)
     else
-        pbtrs!('L', size(A, 1), bandwidth(A.factors,1), bandeddata(A.factors), B)
+        pbtrs!('L', size(A, 1), bandwidth(A.factors,1), _lower_bandeddata(A.factors), B)
     end
     B
 end
@@ -96,9 +103,18 @@ ldiv!(A::Cholesky{T,<:AbstractBandedMatrix}, B::StridedVecOrMat{T}) where T<:Bla
 ArrayLayouts._cholesky(::Union{SymmetricLayout{<:AbstractBandedLayout},HermitianLayout{<:AbstractBandedLayout}}, ax, A, ::Val{false}=Val(false); check::Bool = true) =
     cholesky!(cholcopy(A); check = check)
 
-cholcopy(A::RealHermSymComplexHerm{<:Any,<:BandedMatrix}) =
-    copyto!(similar(A, LinearAlgebra.choltype(A)), A)
+function _banded_cholcopy(A)
+    B = similar(A, LinearAlgebra.choltype(A))
+    # Copying the wrapper may leave some parent storage untouched. Initialize it
+    # for non-isbits types so full-parent copies do not read undefined references.
+    if !isbitstype(eltype(B))
+        fill!(parent(B), zero(eltype(B)))
+    end
+    copyto!(B, A)
+    B
+end
+cholcopy(A::RealHermSymComplexHerm{<:Any,<:BandedMatrix}) = _banded_cholcopy(A)
 cholcopy(A::RealHermSymComplexHerm{<:Any,<:SubArray{<:Number, 2, <:BandedMatrix}}) =
-    copyto!(similar(A, LinearAlgebra.choltype(A)), A)
+    _banded_cholcopy(A)
 # cholesky(A::Hermitian{T,<:BandedMatrix{T}},
 #     ::Val{false}=Val(false); check::Bool = true) where T = cholesky!(cholcopy(A); check = check)

test/test_symbanded.jl

@@ -398,6 +398,57 @@ end
         @test Ac\b ≈ A\b
     end
 
+    @testset "lower Cholesky of a matrix that also stores upper bands" begin
+        for T in (Float64, ComplexF64, BigFloat), u in (1, 2, 3)
+            M = brand(T, 8, 8, 2, 2)
+            A = BandedMatrix(Matrix(M + M') + 10I, (2, u))
+            S = T <: Complex ? Hermitian(A, :L) : Symmetric(A, :L)
+            Ac = cholesky(S)
+
+            @test Ac isa Cholesky{T,<:BandedMatrix{T}}
+            @test Matrix(Ac.L) ≈ cholesky(Matrix(S)).L
+
+            b = rand(T, size(A, 1))
+            @test Ac\b ≈ Matrix(S)\b
+        end
+    end
+
+    @testset "upper Cholesky of a matrix that also stores lower bands" begin
+        for T in (Float64, ComplexF64, BigFloat), l in (1, 2, 3)
+            M = brand(T, 8, 8, 2, 2)
+            A = BandedMatrix(Matrix(M + M') + 10I, (l, 2))
+            S = T <: Complex ? Hermitian(A, :U) : Symmetric(A, :U)
+            Ac = cholesky(S)
+
+            @test Ac isa Cholesky{T,<:BandedMatrix{T}}
+            @test Matrix(Ac.U) ≈ cholesky(Matrix(S)).U
+
+            b = rand(T, size(A, 1))
+            @test Ac\b ≈ Matrix(S)\b
+        end
+    end
+
+    @testset "Cholesky copy initializes non-isbits storage" begin
+        for wrapper in (Symmetric, Hermitian),
+            bw in ((2, 1), (2, 3), (1, 2), (3, 2)),
+            uplo in (:L, :U), useview in (false, true)
+
+            A = BandedMatrix{BigFloat}(I, (8, 8), bw)
+            P = useview ? view(A, :, :) : A
+            S = wrapper(P, uplo)
+
+            C = LinearAlgebra.cholcopy(S)
+            data = BandedMatrices.bandeddata(parent(C))
+            # Inspect references without reading any undefined entries.
+            @test all(i -> isassigned(data, i), eachindex(data))
+
+            F = cholesky(S)
+            identity = Matrix{BigFloat}(I, 8, 8)
+            @test Matrix(F.L) ≈ identity
+            @test Matrix(F.U) ≈ identity
+        end
+    end
+
     B = BandedMatrix(1=>Float64[1:4;], 0=>Float64[1:5;], -1=>Float64[1:4;])
     @test isposdef(B) == isposdef(Matrix(B)) == false
     B = BandedMatrix(1=>Float64[1:4;], 0=>Float64[2:2:10;], -1=>Float64[1:4;])