Coercing constructors & convert (#212)

This adds a "constructor cascade" that allows one to coerce the inputs
to the chosen type. This is frequently useful in generic programming
where you want an object to be of the same type as something else,
for example when appending to a list of objects.

Co-authored-by: David Widmann <[email protected]>

Author: timholy

Project.toml

@@ -1,6 +1,6 @@
 name = "PDMats"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.11.31"
+version = "0.11.32"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/chol.jl

@@ -73,7 +73,7 @@ function invquad(A::Cholesky, x::AbstractVector)
     @check_argdims size(A, 1) == size(x, 1)
     return sum(abs2, chol_lower(A) \ x)
 end
-function invquad(A::Cholesky, X::AbstractMatrix) 
+function invquad(A::Cholesky, X::AbstractMatrix)
     @check_argdims size(A, 1) == size(X, 1)
     Z = chol_lower(A) \ X
     return vec(sum(abs2, Z; dims=1))

src/pdmat.jl

@@ -1,23 +1,54 @@
 """
 Full positive definite matrix together with a Cholesky factorization object.
 """
-struct PDMat{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T}
+struct PDMat{T<:Real,S<:AbstractMatrix{T}} <: AbstractPDMat{T}
     mat::S
     chol::Cholesky{T,S}
 
-    PDMat{T,S}(m::AbstractMatrix{T},c::Cholesky{T,S}) where {T,S} = new{T,S}(m,c)
-end
-
-function PDMat(mat::AbstractMatrix,chol::Cholesky{T,S}) where {T,S}
-    d = LinearAlgebra.checksquare(mat)
-    if size(chol, 1) != d
-        throw(DimensionMismatch("Dimensions of mat and chol are inconsistent."))
+    function PDMat{T,S}(m::AbstractMatrix, c::Cholesky) where {T<:Real,S<:AbstractMatrix{T}}
+        d = LinearAlgebra.checksquare(m)
+        if size(c, 1) != d
+            throw(DimensionMismatch("Dimensions of mat and chol are inconsistent."))
+        end
+        # in principle we might want to check that `c` is a Cholesky factorization of `m`,
+        # but that's slow
+        return new{T,S}(m,c)
     end
-    PDMat{T,S}(convert(S, mat), chol)
 end
+function PDMat{T}(m::AbstractMatrix, c::Cholesky) where T<:Real
+    c = convert(Cholesky{T}, c)
+    return PDMat{T,typeof(c.factors)}(m, c)
+end
+PDMat(mat::AbstractMatrix,chol::Cholesky{T,S}) where {T<:Real,S<:AbstractMatrix{T}} = PDMat{T,S}(mat, chol)
 
+# Construction from another PDMat
+PDMat{T,S}(pdm::PDMat{T,S}) where {T<:Real,S<:AbstractMatrix{T}} = pdm  # since PDMat doesn't support `setindex!` it's not mutable (xref https://docs.julialang.org/en/v1/manual/conversion-and-promotion/#Mutable-collections)
+PDMat{T,S}(pdm::PDMat) where {T<:Real,S<:AbstractMatrix{T}} = PDMat{T,S}(pdm.mat, pdm.chol)
+PDMat{T}(pdm::PDMat{T}) where T<:Real = pdm
+PDMat{T}(pdm::PDMat) where T<:Real = PDMat{T}(pdm.mat, pdm.chol)
+PDMat(pdm::PDMat) = pdm
+
+# Construction from an AbstractMatrix
+function PDMat{T,S}(mat::AbstractMatrix) where {T<:Real,S<:AbstractMatrix{T}}
+    mat = convert(S, mat)
+    return PDMat{T,S}(mat, cholesky(mat))
+end
+function PDMat{T}(mat::AbstractMatrix) where T<:Real
+    mat = convert(AbstractMatrix{T}, mat)
+    return PDMat{T}(mat, cholesky(mat))
+end
 PDMat(mat::AbstractMatrix) = PDMat(mat, cholesky(mat))
-PDMat(fac::Cholesky) = PDMat(AbstractMatrix(fac), fac)
+
+# Construction from a Cholesky factorization
+function PDMat{T,S}(c::Cholesky) where {T<:Real,S<:AbstractMatrix{T}}
+    c = convert(Cholesky{T,S}, c)
+    return PDMat{T,S}(AbstractMatrix(c), c)
+end
+function PDMat{T}(c::Cholesky) where T<:Real
+    c = convert(Cholesky{T}, c)
+    return PDMat{T}(AbstractMatrix(c), c)
+end
+PDMat(c::Cholesky) = PDMat(AbstractMatrix(c), c)
 
 function Base.getproperty(a::PDMat, s::Symbol)
     if s === :dim
@@ -30,13 +61,11 @@ Base.propertynames(::PDMat) = (:mat, :chol, :dim)
 AbstractPDMat(A::Cholesky) = PDMat(A)
 
 ### Conversion
-Base.convert(::Type{PDMat{T}}, a::PDMat{T}) where {T<:Real} = a
-function Base.convert(::Type{PDMat{T}}, a::PDMat) where {T<:Real}
-    chol = convert(Cholesky{T}, a.chol)
-    S = typeof(chol.factors)
-    mat = convert(S, a.mat)
-    return PDMat{T,S}(mat, chol)
-end
+# This next method isn't needed because PDMat{T}(a) returns `a` directly
+# Base.convert(::Type{PDMat{T}}, a::PDMat{T}) where {T<:Real} = a
+Base.convert(::Type{PDMat{T}}, a::PDMat) where {T<:Real} = PDMat{T}(a)
+Base.convert(::Type{PDMat{T,S}}, a::PDMat) where {T<:Real,S<:AbstractMatrix{T}} = PDMat{T,S}(a)
+
 Base.convert(::Type{AbstractPDMat{T}}, a::PDMat) where {T<:Real} = convert(PDMat{T}, a)
 
 ### Basics
@@ -55,7 +84,7 @@ Base.IndexStyle(::Type{PDMat{T,S}}) where {T,S} = Base.IndexStyle(S)
 # Linear Indexing
 Base.@propagate_inbounds Base.getindex(a::PDMat, i::Int) = getindex(a.mat, i)
 # Cartesian Indexing
-Base.@propagate_inbounds Base.getindex(a::PDMat, I::Vararg{Int, 2}) = getindex(a.mat, I...)
+Base.@propagate_inbounds Base.getindex(a::PDMat, i::Int, j::Int) = getindex(a.mat, i, j)
 
 ### Arithmetics
 

test/pdmtypes.jl

@@ -18,13 +18,36 @@ using Test
             M = convert(Array{T,2}, [4. -2. -1.; -2. 5. -1.; -1. -1. 6.])
             V = convert(Array{T,1}, [1.5, 2.5, 2.0])
             X = convert(T,2.0)
+            f64M = Float64.(M)
 
             @testset "PDMat from Matrix" begin
+                pdf64M = PDMat(f64M)
                 test_pdmat(PDMat(M), M,                        cmat_eq=true, verbose=1)
+                test_pdmat(PDMat{Float64}(M), f64M,            cmat_eq=true, verbose=1)
+                test_pdmat(PDMat{Float64,Matrix{Float64}}(M), f64M, cmat_eq=true, verbose=1)
+                @test_throws TypeError PDMat{Float32,Matrix{Float64}}(M)
+            end
+            @testset "PDMat from PDMat" begin
+                pdM = PDMat(M)
+                pdf64M = PDMat(f64M)
+                @test pdM === PDMat(pdM)
+                @test pdf64M === PDMat{Float64}(pdf64M) === PDMat{Float64,Matrix{Float64}}(pdf64M)
+                test_pdmat(PDMat(pdM), M,                      cmat_eq=true, verbose=1)
+                test_pdmat(PDMat{Float64}(pdf64M), f64M,       cmat_eq=true, verbose=1)
+                test_pdmat(PDMat{Float64,Matrix{Float64}}(pdf64M), f64M, cmat_eq=true, verbose=1)
+                if Base.VERSION >= v"1.12.0-DEV.1654"   # julia #56562
+                    @test isa(convert(typeof(pdf64M), pdM), typeof(pdf64M))
+                end
+                @test_throws TypeError PDMat{Float32,Matrix{Float64}}(pdM)
             end
             @testset "PDMat from Cholesky" begin
                 cholL = Cholesky(Matrix(transpose(cholesky(M).factors)), 'L', 0)
+                cholLf64 = Cholesky(Matrix(transpose(cholesky(f64M).factors)), 'L', 0)
                 test_pdmat(PDMat(cholL), M,                    cmat_eq=true, verbose=1)
+                test_pdmat(PDMat{Float64}(cholLf64), f64M,     cmat_eq=true, verbose=1)
+                if Base.VERSION >= v"1.12.0-DEV.1654"   # julia #56562
+                    test_pdmat(PDMat{Float64,Matrix{Float64}}(cholLf64), f64M, cmat_eq=true, verbose=1)
+                end
             end
             @testset "PDiagMat" begin
                 test_pdmat(PDiagMat(V), Matrix(Diagonal(V)),   cmat_eq=true, verbose=1)
@@ -142,7 +165,7 @@ using Test
         # right division not defined for CHOLMOD:
         # `rdiv!(::Matrix{Float64}, ::SuiteSparse.CHOLMOD.Factor{Float64})` not defined
         if !HAVE_CHOLMOD
-            z = x / PDSparseMat(sparse(first(A), 1, 1)) 
+            z = x / PDSparseMat(sparse(first(A), 1, 1))
             @test typeof(z) === typeof(y)
             @test size(z) == size(y)
             @test z ≈ y

test/specialarrays.jl

@@ -15,6 +15,14 @@ using StaticArrays
         @test PDMat(S, C) === PDS
         @test @allocated(PDMat(S)) == @allocated(PDMat(C)) == @allocated(PDMat(S, C))
 
+        if Base.VERSION >= v"1.12.0-DEV.1654"    # julia #56562
+            A = PDMat(Matrix{Float64}(I, 2, 2))
+            B = PDMat(SMatrix{2,2,Float64}(I))
+            @test !isa(A.mat, typeof(B.mat))
+            S = convert(typeof(B), A)
+            @test  isa(S.mat, typeof(B.mat))
+        end
+
         # Diagonal matrix
         D = PDiagMat(@SVector(rand(4)))
         @test D isa PDiagMat{Float64, <:SVector{4, Float64}}