From d8e10c7edead31916ab2d970b338438033022119 Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Mon, 19 Aug 2024 11:07:57 +0100 Subject: [PATCH 1/2] Fix layout of `cholesky(Symmetric(view(BandedMatrix)))` (#451) * Fix 450 * Change tests for LTS * Change tests for LTS * Should also test that U is actually correct * Fix check of LTS --- Project.toml | 2 +- src/symbanded/BandedCholesky.jl | 2 ++ test/test_symbanded.jl | 17 +++++++++++++++++ 3 files changed, 20 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 388aa786..a1238e8b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "BandedMatrices" uuid = "aae01518-5342-5314-be14-df237901396f" -version = "1.7.2" +version = "1.7.3" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" diff --git a/src/symbanded/BandedCholesky.jl b/src/symbanded/BandedCholesky.jl index 5f25e811..0e8549a1 100644 --- a/src/symbanded/BandedCholesky.jl +++ b/src/symbanded/BandedCholesky.jl @@ -98,5 +98,7 @@ ArrayLayouts._cholesky(::Union{SymmetricLayout{<:AbstractBandedLayout},Hermitian cholcopy(A::RealHermSymComplexHerm{<:Any,<:BandedMatrix}) = copyto!(similar(A, LinearAlgebra.choltype(A)), A) +cholcopy(A::RealHermSymComplexHerm{<:Any,<:SubArray{<:Number, 2, <:BandedMatrix}}) = + copyto!(similar(A, LinearAlgebra.choltype(A)), 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 c1ac0fb1..2cf802a3 100644 --- a/test/test_symbanded.jl +++ b/test/test_symbanded.jl @@ -402,6 +402,23 @@ end @test isposdef(B) == isposdef(Matrix(B)) == false B = BandedMatrix(1=>Float64[1:4;], 0=>Float64[2:2:10;], -1=>Float64[1:4;]) @test isposdef(B) == isposdef(Matrix(B)) == true + + @testset "Issue #450: Cholesky of a BandedMatrix's view" begin + L = brand(10, 10, 2, 0) + B = BandedMatrix(Symmetric(L * L')); + Bv = Symmetric(view(B, :, :)) + chol = cholesky(Bv) + if VERSION >= v"1.7" + @test MemoryLayout(chol.U) isa TriangularLayout{'U', 'N', typeof(MemoryLayout(B))} + @test MemoryLayout(chol.L) isa TriangularLayout{'L', 'N', typeof(MemoryLayout(B'))} + @test isbanded(chol.L) + end + @test isbanded(chol.U) + cc = LinearAlgebra.cholcopy(Bv) + @test cc isa Symmetric{Float64, typeof(B)} + @test chol.U ≈ cholesky(Matrix(B)).U + @test cc == B + end end end # module From c6c39129db2662bbb5e8fa24263ff5204fdf9777 Mon Sep 17 00:00:00 2001 From: Maxim Vassiliev <76599693+max-vassili3v@users.noreply.github.com> Date: Wed, 21 Aug 2024 16:37:19 +0100 Subject: [PATCH 2/2] Colon-integer indexing bug (#452) * fix bug * fix unit tests --- src/banded/BandedMatrix.jl | 4 ++-- test/test_indexing.jl | 32 ++++++++++++++++++-------------- 2 files changed, 20 insertions(+), 16 deletions(-) diff --git a/src/banded/BandedMatrix.jl b/src/banded/BandedMatrix.jl index 4d1909fb..2137bfbe 100644 --- a/src/banded/BandedMatrix.jl +++ b/src/banded/BandedMatrix.jl @@ -420,7 +420,7 @@ end @propagate_inbounds function getindex(A::BandedMatrix, ::Colon, j::Int) @boundscheck checkbounds(A, axes(A,1), j) r = similar(A, axes(A,1)) - r[firstindex(r):colstart(A,j)-1] .= zero(eltype(r)) + r[firstindex(r):min(size(A, 1), colstart(A,j)-1)] .= zero(eltype(r)) # broadcasted assignment is currently faster than setindex # see https://github.com/JuliaLang/julia/issues/40962#issuecomment-1921340377 # may need revisiting in the future @@ -439,7 +439,7 @@ end @propagate_inbounds function getindex(A::BandedMatrix, k::Int, ::Colon) @boundscheck checkbounds(A, k, axes(A,2)) r = similar(A, axes(A,2)) - r[firstindex(r):rowstart(A,k)-1] .= zero(eltype(r)) + r[firstindex(r):min(size(A, 2), rowstart(A,k)-1)] .= zero(eltype(r)) r[rowrange(A,k)] = @view A.data[data_rowrange(A,k)] r[rowstop(A,k)+1:end] .= zero(eltype(r)) return r diff --git a/test/test_indexing.jl b/test/test_indexing.jl index f514453d..e0a94a88 100644 --- a/test/test_indexing.jl +++ b/test/test_indexing.jl @@ -208,7 +208,7 @@ import BandedMatrices: rowstart, rowstop, colstart, colstop, end @testset "vector - BandRange/Colon - integer" begin - a = BandedMatrix(Ones{Int}(5, 7), (2, 1)) + a = BandedMatrix(Ones{Int}(5, 8), (2, 1)) # 5x7 BandedMatrices.BandedMatrix{Float64}: # 1.0 1.0 0 0 0 0 0 0 # 1.0 1.0 1.0 0 0 0 0 0 @@ -224,11 +224,11 @@ import BandedMatrices: rowstart, rowstop, colstart, colstop, a[BandRange, 5] = [15, 16] a[BandRange, 6] = [17] - @test a == [ 1 4 0 0 0 0 0; - 2 5 8 0 0 0 0; - 3 6 9 12 0 0 0; - 0 7 10 13 15 0 0; - 0 0 11 14 16 17 0] + @test a == [ 1 4 0 0 0 0 0 0; + 2 5 8 0 0 0 0 0; + 3 6 9 12 0 0 0 0; + 0 7 10 13 15 0 0 0; + 0 0 11 14 16 17 0 0] @test a[BandRange, 1] == @view(a[BandRange, 1]) == [1, 2, 3] @test a[BandRange, 2] == @view(a[BandRange, 2]) == [4, 5, 6, 7] @@ -237,6 +237,7 @@ import BandedMatrices: rowstart, rowstop, colstart, colstop, @test a[BandRange, 5] == @view(a[BandRange, 5]) == [15, 16] @test a[BandRange, 6] == @view(a[BandRange, 6]) == [17] @test a[BandRange, 7] == @view(a[BandRange, 7]) == Int[] + @test a[BandRange, 8] == @view(a[BandRange, 8]) == Int[] @test a[:, 1] == view(a, :, 1) == [1,2,3,0,0] @test a[:, 2] == view(a, :, 2) == [4,5,6,7,0] @@ -245,11 +246,12 @@ import BandedMatrices: rowstart, rowstop, colstart, colstop, @test a[:, 5] == view(a, :, 5) == [0,0,0,15,16] @test a[:, 6] == view(a, :, 6) == [0,0,0,0,17] @test a[:, 7] == view(a, :, 7) == [0,0,0,0,0] + @test a[:, 8] == view(a, :, 8) == [0,0,0,0,0] @test_throws BoundsError a[:, 0] = [1, 2, 3] @test_throws DimensionMismatch a[:, 1] = [1, 2, 3] @test_throws BoundsError a[BandRange, 0] = [1, 2, 3] - @test_throws BoundsError a[BandRange, 8] = [1, 2, 3] + @test_throws BoundsError a[BandRange, 9] = [1, 2, 3] @test_throws DimensionMismatch a[BandRange, 1] = [1, 2] end @@ -299,7 +301,7 @@ import BandedMatrices: rowstart, rowstop, colstart, colstop, end @testset "vector - integer - BandRange/Colon" begin - a = BandedMatrix(Ones{Int}(7, 5), (1, 2)) + a = BandedMatrix(Ones{Int}(8, 5), (1, 2)) # 5x7 BandedMatrices.BandedMatrix{Float64}: # 1.0 1.0 0 0 0 0 0 0 # 1.0 1.0 1.0 0 0 0 0 0 @@ -315,11 +317,11 @@ import BandedMatrices: rowstart, rowstop, colstart, colstop, a[5, BandRange] = [15, 16] a[6, BandRange] = [17] - @test a == [ 1 4 0 0 0 0 0; - 2 5 8 0 0 0 0; - 3 6 9 12 0 0 0; - 0 7 10 13 15 0 0; - 0 0 11 14 16 17 0]' + @test a == [ 1 4 0 0 0 0 0 0; + 2 5 8 0 0 0 0 0; + 3 6 9 12 0 0 0 0; + 0 7 10 13 15 0 0 0; + 0 0 11 14 16 17 0 0]' @test a[1, BandRange] == @view(a[1, BandRange]) == [1, 2, 3] @test a[2, BandRange] == @view(a[2, BandRange]) == [4, 5, 6, 7] @@ -328,6 +330,7 @@ import BandedMatrices: rowstart, rowstop, colstart, colstop, @test a[5, BandRange] == @view(a[5, BandRange]) == [15, 16] @test a[6, BandRange] == @view(a[6, BandRange]) == [17] @test a[7, BandRange] == @view(a[7, BandRange]) == Int[] + @test a[8, BandRange] == @view(a[7, BandRange]) == Int[] @test a[1, :] == @view(a[1, :]) == [1,2,3,0,0] @test a[2, :] == @view(a[2, :]) == [4,5,6,7,0] @@ -336,11 +339,12 @@ import BandedMatrices: rowstart, rowstop, colstart, colstop, @test a[5, :] == @view(a[5, :]) == [0,0,0,15,16] @test a[6, :] == @view(a[6, :]) == [0,0,0,0,17] @test a[7, :] == @view(a[7, :]) == [0,0,0,0,0] + @test a[8, :] == @view(a[7, :]) == [0,0,0,0,0] @test_throws BoundsError a[0, :] = [1, 2, 3] @test_throws DimensionMismatch a[1, :] = [1, 2, 3] @test_throws BoundsError a[0, BandRange] = [1, 2, 3] - @test_throws BoundsError a[8, BandRange] = [1, 2, 3] + @test_throws BoundsError a[9, BandRange] = [1, 2, 3] @test_throws DimensionMismatch a[1, BandRange] = [1, 2] end