Add FillArrays extension & bump Julia compat to >=1.6

.github/workflows/ci.yml

@@ -12,7 +12,6 @@ jobs:
  fail-fast: false
  matrix:
  version:
- - '1.0'
  - '1.6'
  - '1'
  - nightly

Project.toml

@@ -3,17 +3,26 @@ uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
 version = "0.11.26"
 
 [deps]
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [compat]
-julia = "1"
+FillArrays = "1"
+julia = "1.6"
+
+[extensions]
+PDMatsFillArraysExt = "FillArrays"
 
 [extras]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BandedMatrices", "StaticArrays", "Test"]
+test = ["BandedMatrices", "FillArrays", "StaticArrays", "Test"]
+
+[weakdeps]
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"

ext/PDMatsFillArraysExt.jl

@@ -0,0 +1,11 @@
+module PDMatsFillArraysExt
+
+using PDMats: PDMats, LinearAlgebra
+using FillArrays: FillArrays
+
+function PDMats.AbstractPDMat(a::LinearAlgebra.Diagonal{T,<:FillArrays.AbstractFillVector{T}}) where {T<:Real}
+ dim = size(a, 1)
+ return PDMats.ScalMat(dim, FillArrays.getindex_value(a.diag))
+end
+
+end # module

src/PDMats.jl

@@ -56,4 +56,7 @@ module PDMats
 
  include("deprecates.jl")
 
+ if !isdefined(Base, :get_extension)
+ include("../ext/PDMatsFillArraysExt.jl")
+ end
 end # module

src/scalmat.jl

@@ -78,7 +78,7 @@ LinearAlgebra.sqrt(a::ScalMat) = ScalMat(a.dim, sqrt(a.value))
 function whiten!(r::AbstractVecOrMat, a::ScalMat, x::AbstractVecOrMat)
  @check_argdims axes(r) == axes(x)
  @check_argdims a.dim == size(x, 1)
- _ldiv!(r, sqrt(a.value), x)
+ ldiv!(r, sqrt(a.value), x)
 end
 
 function unwhiten!(r::AbstractVecOrMat, a::ScalMat, x::AbstractVecOrMat)
@@ -174,10 +174,10 @@ end
 
 function X_invA_Xt(a::ScalMat, x::Matrix)
  @check_argdims LinearAlgebra.checksquare(a) == size(x, 2)
- _rdiv!(x * transpose(x), a.value)
+ rdiv!(x * transpose(x), a.value)
 end
 
 function Xt_invA_X(a::ScalMat, x::Matrix)
  @check_argdims LinearAlgebra.checksquare(a) == size(x, 1)
- _rdiv!(transpose(x) * x, a.value)
+ rdiv!(transpose(x) * x, a.value)
 end

src/utils.jl

@@ -103,29 +103,3 @@ function colwise_sumsqinv!(r::AbstractArray, a::AbstractMatrix, c::Real)
  end
  return r
 end
-
-# `rdiv!(::AbstractArray, ::Number)` was introduced in Julia 1.2
-# https://github.com/JuliaLang/julia/pull/31179
-@static if VERSION < v"1.2.0-DEV.385"
- function _rdiv!(X::AbstractArray, s::Number)
- @simd for I in eachindex(X)
- @inbounds X[I] /= s
- end
- X
- end
-else
- _rdiv!(X::AbstractArray, s::Number) = rdiv!(X, s)
-end
-
-# `ldiv!(::AbstractArray, ::Number, ::AbstractArray)` was introduced in Julia 1.4
-# https://github.com/JuliaLang/julia/pull/33806
-@static if VERSION < v"1.4.0-DEV.635"
- _ldiv!(Y::AbstractArray, s::Number, X::AbstractArray) = Y .= s .\ X
-else
- _ldiv!(Y::AbstractArray, s::Number, X::AbstractArray) = ldiv!(Y, s, X)
-end
-
-# https://github.com/JuliaLang/julia/pull/29749
-if VERSION < v"1.1.0-DEV.792"
- eachcol(A::AbstractVecOrMat) = (view(A, :, i) for i in axes(A, 2))
-end

test/ext.jl

@@ -0,0 +1,10 @@
+using FillArrays
+
+@testset "PDMatsFillArraysExt" begin
+ for diag in (Ones(5), Fill(4.1, 8))
+ a = @inferred(AbstractPDMat(Diagonal(diag)))
+ @test a isa ScalMat
+ @test a.dim == length(diag)
+ @test a.value == first(diag)
+ end
+end

test/pdmtypes.jl

@@ -134,13 +134,10 @@ using Test
  @test z ≈ y
  end
 
- # requires https://github.com/JuliaLang/julia/pull/32594
- if VERSION >= v"1.3.0-DEV.562"
- z = x / PDMat(A)
- @test typeof(z) === typeof(y)
- @test size(z) == size(y)
- @test z ≈ y
- end
+ z = x / PDMat(A)
+ @test typeof(z) === typeof(y)
+ @test size(z) == size(y)
+ @test z ≈ y
 
  # right division not defined for CHOLMOD:
  # `rdiv!(::Matrix{Float64}, ::SuiteSparse.CHOLMOD.Factor{Float64})` not defined
@@ -153,15 +150,11 @@ using Test
  end
 
  @testset "PDMat from Cholesky decomposition of diagonal matrix (#137)" begin
- # U'*U where U isa UpperTriangular etc.
- # requires https://github.com/JuliaLang/julia/pull/33334
- if VERSION >= v"1.4.0-DEV.286"
- x = rand(10, 10)
- A = Diagonal(x' * x)
- M = PDMat(cholesky(A))
- @test M isa PDMat{Float64, typeof(A)}
- @test Matrix(M) ≈ A
- end
+ x = rand(10, 10)
+ A = Diagonal(x' * x)
+ M = PDMat(cholesky(A))
+ @test M isa PDMat{Float64, typeof(A)}
+ @test Matrix(M) ≈ A
  end
 
  @testset "AbstractPDMat constructors (#136)" begin
@@ -192,12 +185,7 @@ using Test
  @test M isa PDSparseMat
  @test Matrix(M) ≈ A
 
- if VERSION < v"1.6"
- # inference fails e.g. on Julia 1.0
- M = AbstractPDMat(cholesky(sparse(A)))
- else
- M = @inferred AbstractPDMat(cholesky(sparse(A)))
- end
+ M = @inferred AbstractPDMat(cholesky(sparse(A)))
  @test M isa PDSparseMat
  @test Matrix(M) ≈ A
  end

test/runtests.jl

@@ -1,5 +1,5 @@
 include("testutils.jl")
-tests = ["pdmtypes", "abstracttypes", "addition", "generics", "kron", "chol", "specialarrays", "sqrt"]
+tests = ["pdmtypes", "abstracttypes", "addition", "generics", "kron", "chol", "specialarrays", "sqrt", "ext"]
 println("Running tests ...")
 
 for t in tests

test/specialarrays.jl

@@ -47,16 +47,16 @@ using StaticArrays
  @test A \ Y ≈ Matrix(A) \ Matrix(Y)
 
  @test whiten(A, x) isa SVector{4, Float64}
- @test whiten(A, x) ≈ cholesky(Matrix(A)).L \ Vector(x)
+ @test whiten(A, x) ≈ cholesky(Symmetric(Matrix(A))).L \ Vector(x)
 
  @test whiten(A, Y) isa SMatrix{4, 10, Float64}
- @test whiten(A, Y) ≈ cholesky(Matrix(A)).L \ Matrix(Y)
+ @test whiten(A, Y) ≈ cholesky(Symmetric(Matrix(A))).L \ Matrix(Y)
 
  @test unwhiten(A, x) isa SVector{4, Float64}
- @test unwhiten(A, x) ≈ cholesky(Matrix(A)).L * Vector(x)
+ @test unwhiten(A, x) ≈ cholesky(Symmetric(Matrix(A))).L * Vector(x)
 
  @test unwhiten(A, Y) isa SMatrix{4, 10, Float64}
- @test unwhiten(A, Y) ≈ cholesky(Matrix(A)).L * Matrix(Y)
+ @test unwhiten(A, Y) ≈ cholesky(Symmetric(Matrix(A))).L * Matrix(Y)
 
  @test quad(A, x) isa Float64
  @test quad(A, x) ≈ Vector(x)' * Matrix(A) * Vector(x)
@@ -95,10 +95,7 @@ using StaticArrays
  Y = rand(5, 2)
  @test P * x ≈ x
  @test P * Y ≈ Y
- # Right division with Cholesky requires https://github.com/JuliaLang/julia/pull/32594
- if VERSION >= v"1.3.0-DEV.562"
- @test X / P ≈ X
- end
+ @test X / P ≈ X
  @test P \ x ≈ x
  @test P \ Y ≈ Y
  @test X_A_Xt(P, X) ≈ X * X'

test/testutils.jl

@@ -237,10 +237,9 @@ function pdtest_div(C, Imat::Matrix, X::Matrix, verbose::Int)
  @assert d == size(C, 1) == size(C, 2)
  @assert size(Imat) == size(C)
  @test C \ X ≈ Imat * X
- # Right division with Choleskyrequires https://github.com/JuliaLang/julia/pull/32594
  # CHOLMOD throws error since no method is found for
  # `rdiv!(::Matrix{Float64}, ::SuiteSparse.CHOLMOD.Factor{Float64})`
- check_rdiv = !(C isa PDMat && VERSION < v"1.3.0-DEV.562") && !(C isa PDSparseMat && HAVE_CHOLMOD)
+ check_rdiv = !(C isa PDSparseMat && HAVE_CHOLMOD)
  check_rdiv && @test Matrix(X') / C ≈ (C \ X)'
 
  for i = 1:n