diff --git a/src/symbanded/BandedCholesky.jl b/src/symbanded/BandedCholesky.jl index 0e8549a1..0d26f876 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 @@ -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) diff --git a/test/test_symbanded.jl b/test/test_symbanded.jl index 7a579e20..d1050a77 100644 --- a/test/test_symbanded.jl +++ b/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;])