From f913946efe7a042a1b682b9291bef225b5fb7641 Mon Sep 17 00:00:00 2001 From: Giuseppe Ragusa Date: Thu, 15 Oct 2026 03:25:05 +0200 Subject: [PATCH 1/8] Fix lower banded Cholesky with upper bands cholesky(Symmetric(A, :L)) and solves with its factor passed bandeddata(A) to LAPACK's pbtrf!/pbtrs! with uplo 'L'. LAPACK's lower band storage starts at the diagonal, which a BandedMatrix with upper bandwidth u stores in row u + 1, so for u > 0 LAPACK read a superdiagonal as the diagonal. Both calls now pass the rows from the diagonal down. Fixes #141 --- src/symbanded/BandedCholesky.jl | 11 +++++++++-- test/test_symbanded.jl | 15 +++++++++++++++ 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/src/symbanded/BandedCholesky.jl b/src/symbanded/BandedCholesky.jl index 0e8549a1..24ec052a 100644 --- a/src/symbanded/BandedCholesky.jl +++ b/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 diff --git a/test/test_symbanded.jl b/test/test_symbanded.jl index 7a579e20..aecb47c2 100644 --- a/test/test_symbanded.jl +++ b/test/test_symbanded.jl @@ -398,6 +398,21 @@ 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 + 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;]) From daf80804796b251d184db6fb9e822d556cc7187f Mon Sep 17 00:00:00 2001 From: Giuseppe Ragusa Date: Thu, 15 Oct 2026 08:40:34 +0200 Subject: [PATCH 2/8] Zero-initialize the banded Cholesky copy cholcopy copied a Symmetric or Hermitian banded matrix into a `similar` one, which leaves the parent's band storage outside the stored triangle uninitialized. For BigFloat those entries are #undef, and on Julia 1.10, where Matrix(::LowerTriangular) copies the whole parent, Matrix(cholesky(Symmetric(A, :L)).L) threw an UndefRefError when A also stored upper bands (likewise for :U with lower bands). cholcopy now zeroes the parent before copying. Assisted-by: Claude Opus 5.5 --- src/symbanded/BandedCholesky.jl | 13 ++++++++++--- test/test_symbanded.jl | 15 +++++++++++++++ 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/src/symbanded/BandedCholesky.jl b/src/symbanded/BandedCholesky.jl index 24ec052a..3fe91f44 100644 --- a/src/symbanded/BandedCholesky.jl +++ b/src/symbanded/BandedCholesky.jl @@ -103,9 +103,16 @@ 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) +# `copyto!` into a `Symmetric`/`Hermitian` need not write every entry of the +# parent's band storage, so the parent is zeroed first: otherwise the factor holds +# uninitialized entries (`#undef` for non-isbits element types such as `BigFloat`). +function _banded_cholcopy(A) + B = similar(A, LinearAlgebra.choltype(A)) + fill!(parent(B), zero(eltype(B))) + copyto!(B, A) +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) diff --git a/test/test_symbanded.jl b/test/test_symbanded.jl index aecb47c2..2ceeaecb 100644 --- a/test/test_symbanded.jl +++ b/test/test_symbanded.jl @@ -413,6 +413,21 @@ end 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 + 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;]) From 30c0326cfdafcd6a661a3528196b9050b5b86252 Mon Sep 17 00:00:00 2001 From: Giuseppe Ragusa Date: Thu, 15 Oct 2026 22:46:54 +0200 Subject: [PATCH 3/8] Improve copying speed for bittypes matrices. --- src/symbanded/BandedCholesky.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/symbanded/BandedCholesky.jl b/src/symbanded/BandedCholesky.jl index 3fe91f44..83c09bbe 100644 --- a/src/symbanded/BandedCholesky.jl +++ b/src/symbanded/BandedCholesky.jl @@ -108,8 +108,12 @@ ArrayLayouts._cholesky(::Union{SymmetricLayout{<:AbstractBandedLayout},Hermitian # uninitialized entries (`#undef` for non-isbits element types such as `BigFloat`). function _banded_cholcopy(A) B = similar(A, LinearAlgebra.choltype(A)) - fill!(parent(B), zero(eltype(B))) + P = parent(B) + u = bandwidth(B, 2) + + fill!(view(P, 1:u, :), zero(eltype(B))) copyto!(B, A) + B end cholcopy(A::RealHermSymComplexHerm{<:Any,<:BandedMatrix}) = _banded_cholcopy(A) cholcopy(A::RealHermSymComplexHerm{<:Any,<:SubArray{<:Number, 2, <:BandedMatrix}}) = From b2e856f1a0b9ed220063dec932bfc0078e5ad17b Mon Sep 17 00:00:00 2001 From: Giuseppe Ragusa Date: Thu, 15 Oct 2026 23:17:09 +0200 Subject: [PATCH 4/8] Fix Cholesky copy initialization for non-isbits types Restore the full-parent fill for non-isbits types while skipping it for isbits types. Add regression coverage for both wrappers, triangle orientations, asymmetric bandwidths, and views. Full package tests pass on Julia 1.10.12 and 1.13.0. --- src/symbanded/BandedCholesky.jl | 12 ++++++------ test/test_symbanded.jl | 21 +++++++++++++++++++++ 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/src/symbanded/BandedCholesky.jl b/src/symbanded/BandedCholesky.jl index 83c09bbe..303c28dc 100644 --- a/src/symbanded/BandedCholesky.jl +++ b/src/symbanded/BandedCholesky.jl @@ -103,15 +103,15 @@ 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) -# `copyto!` into a `Symmetric`/`Hermitian` need not write every entry of the -# parent's band storage, so the parent is zeroed first: otherwise the factor holds -# uninitialized entries (`#undef` for non-isbits element types such as `BigFloat`). function _banded_cholcopy(A) B = similar(A, LinearAlgebra.choltype(A)) - P = parent(B) - u = bandwidth(B, 2) - fill!(view(P, 1:u, :), zero(eltype(B))) + # 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 diff --git a/test/test_symbanded.jl b/test/test_symbanded.jl index 2ceeaecb..d1050a77 100644 --- a/test/test_symbanded.jl +++ b/test/test_symbanded.jl @@ -428,6 +428,27 @@ end 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;]) From 20160e8b49859ba09e4b49b4665a5ce04f91be96 Mon Sep 17 00:00:00 2001 From: Giuseppe Ragusa Date: Thu, 15 Oct 2026 23:44:45 +0200 Subject: [PATCH 5/8] Apply suggestion from @dlfivefifty Co-authored-by: Sheehan Olver --- src/symbanded/BandedCholesky.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/symbanded/BandedCholesky.jl b/src/symbanded/BandedCholesky.jl index 303c28dc..947b2fdd 100644 --- a/src/symbanded/BandedCholesky.jl +++ b/src/symbanded/BandedCholesky.jl @@ -104,6 +104,9 @@ ArrayLayouts._cholesky(::Union{SymmetricLayout{<:AbstractBandedLayout},Hermitian cholesky!(cholcopy(A); check = check) function _banded_cholcopy(A) + zero!(view(parent(B), 1:bandwidth(B, 2), :)) + copyto!(similar(A, LinearAlgebra.choltype(A)), A) +end B = similar(A, LinearAlgebra.choltype(A)) # Copying the wrapper may leave some parent storage untouched. Initialize it From 3f430c55f52985adfa93bc7b3038b6816d6e3fa0 Mon Sep 17 00:00:00 2001 From: Giuseppe Ragusa Date: Thu, 15 Oct 2026 23:45:47 +0200 Subject: [PATCH 6/8] Refactor _banded_cholcopy to simplify logic Removed unnecessary initialization and copying steps in the _banded_cholcopy function. --- src/symbanded/BandedCholesky.jl | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/symbanded/BandedCholesky.jl b/src/symbanded/BandedCholesky.jl index 947b2fdd..bb9243f6 100644 --- a/src/symbanded/BandedCholesky.jl +++ b/src/symbanded/BandedCholesky.jl @@ -106,17 +106,6 @@ ArrayLayouts._cholesky(::Union{SymmetricLayout{<:AbstractBandedLayout},Hermitian function _banded_cholcopy(A) zero!(view(parent(B), 1:bandwidth(B, 2), :)) copyto!(similar(A, LinearAlgebra.choltype(A)), A) -end - 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}}) = From 403c5281cb5698347c5280d8e74642ccc69f6210 Mon Sep 17 00:00:00 2001 From: Giuseppe Ragusa Date: Fri, 16 Oct 2026 01:22:52 +0200 Subject: [PATCH 7/8] Fix zeroing in _banded_cholcopy function --- src/symbanded/BandedCholesky.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/symbanded/BandedCholesky.jl b/src/symbanded/BandedCholesky.jl index bb9243f6..7693f155 100644 --- a/src/symbanded/BandedCholesky.jl +++ b/src/symbanded/BandedCholesky.jl @@ -104,7 +104,7 @@ ArrayLayouts._cholesky(::Union{SymmetricLayout{<:AbstractBandedLayout},Hermitian cholesky!(cholcopy(A); check = check) function _banded_cholcopy(A) - zero!(view(parent(B), 1:bandwidth(B, 2), :)) + zero!(view(parent(A), 1:bandwidth(A, 2), :)) copyto!(similar(A, LinearAlgebra.choltype(A)), A) end cholcopy(A::RealHermSymComplexHerm{<:Any,<:BandedMatrix}) = _banded_cholcopy(A) From e7e2645f3ae78db45096a8e6e456bc749a24d833 Mon Sep 17 00:00:00 2001 From: Giuseppe Ragusa Date: Fri, 16 Oct 2026 01:34:49 +0200 Subject: [PATCH 8/8] Enhance _banded_cholcopy to handle non-isbits types Initialize parent storage for non-isbits types in _banded_cholcopy function. --- src/symbanded/BandedCholesky.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/symbanded/BandedCholesky.jl b/src/symbanded/BandedCholesky.jl index 7693f155..0d26f876 100644 --- a/src/symbanded/BandedCholesky.jl +++ b/src/symbanded/BandedCholesky.jl @@ -104,8 +104,14 @@ ArrayLayouts._cholesky(::Union{SymmetricLayout{<:AbstractBandedLayout},Hermitian cholesky!(cholcopy(A); check = check) function _banded_cholcopy(A) - zero!(view(parent(A), 1:bandwidth(A, 2), :)) - copyto!(similar(A, LinearAlgebra.choltype(A)), 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}}) =