diff --git a/Project.toml b/Project.toml
index 50b4806..9754928 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "PDMats"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.11.31"
+version = "0.12.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -14,11 +14,19 @@ StaticArrays = "1"
 Test = "<0.0.1, 1"
 julia = "1"
 
+[extensions]
+PDMatsSparseArraysExt = "SparseArrays"
+
 [extras]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BandedMatrices", "StaticArrays", "Random", "Test"]
+test = ["BandedMatrices", "Random", "SparseArrays", "StaticArrays", "Test"]
+
+[weakdeps]
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
diff --git a/README.md b/README.md
index f48baab..13b8ad8 100644
--- a/README.md
+++ b/README.md
@@ -27,16 +27,16 @@ Elemenent types are in princple all Real types, but in practice this is limited
 * `PDMat`: full covariance matrix, defined as
 
 ```julia
-struct PDMat{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T}
+struct PDMat{T<:Real,S<:AbstractMatrix{T},F<:Factorization} <: AbstractPDMat{T}
     mat::S                      # input matrix
-    chol::Cholesky{T,S}         # Cholesky factorization of mat
+    fact::F                     # factorization of mat
 end
 
 # Constructors
 
 PDMat(mat, chol)    # with both the input matrix and a Cholesky factorization
 
-PDMat(mat)          # with the input matrix, of type Matrix or Symmetric
+PDMat(mat)          # with the input matrix
                     # Remarks: the Cholesky factorization will be computed
                     # upon construction.
 
@@ -73,29 +73,6 @@ end
 ScalMat(d, v)        # with dimension d and diagonal value v
 ```
 
-
-* `PDSparseMat`: sparse covariance matrix, defined as
-
-```julia
-struct PDSparseMat{T<:Real,S<:AbstractSparseMatrix} <: AbstractPDMat{T}
-    mat::SparseMatrixCSC           # input matrix
-    chol::CholTypeSparse           # Cholesky factorization of mat
-end
-
-# Constructors
-
-PDSparseMat(mat, chol)    # with both the input matrix and a Cholesky factorization
-
-PDSparseMat(mat)          # with the sparse input matrix, of type SparseMatrixCSC
-                          # Remarks: the Cholesky factorization will be computed
-                          # upon construction.
-
-PDSparseMat(chol)         # with the Cholesky factorization
-                          # Remarks: the sparse matrix 'mat' will be computed upon
-                          # construction.
-```
-
-
 ## Common interface
 
 All subtypes of `AbstractPDMat` share the same API, *i.e.* with the same set of methods to operate on their instances. These methods are introduced below, where `a` is an instance of a subtype of `AbstractPDMat` to represent a positive definite matrix, `x` be a column vector or a matrix with `size(x,1) == size(a, 1) == size(a, 2)`, and `c` be a positive real value.
diff --git a/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl b/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl
new file mode 100644
index 0000000..e3fa624
--- /dev/null
+++ b/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl
@@ -0,0 +1,25 @@
+module PDMatsSparseArraysExt
+
+using PDMats
+using SparseArrays
+
+using PDMats.LinearAlgebra
+
+if isdefined(Base, :get_extension)
+    const HAVE_CHOLMOD = isdefined(SparseArrays, :CHOLMOD)
+else
+    import SuiteSparse
+    const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD)
+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
+
+if HAVE_CHOLMOD
+    include("chol.jl")
+    include("pdsparsemat.jl")
+end
+
+end # module
diff --git a/ext/PDMatsSparseArraysExt/chol.jl b/ext/PDMatsSparseArraysExt/chol.jl
new file mode 100644
index 0000000..195cf4d
--- /dev/null
+++ b/ext/PDMatsSparseArraysExt/chol.jl
@@ -0,0 +1,9 @@
+if isdefined(Base, :get_extension)
+    const CholTypeSparse{T} = SparseArrays.CHOLMOD.Factor{T}
+else
+    const CholTypeSparse{T} = SuiteSparse.CHOLMOD.Factor{T}
+end
+
+# Take into account pivoting!
+PDMats.chol_lower(cf::CholTypeSparse) = cf.PtL
+PDMats.chol_upper(cf::CholTypeSparse) = cf.UP
diff --git a/ext/PDMatsSparseArraysExt/pdsparsemat.jl b/ext/PDMatsSparseArraysExt/pdsparsemat.jl
new file mode 100644
index 0000000..129a33b
--- /dev/null
+++ b/ext/PDMatsSparseArraysExt/pdsparsemat.jl
@@ -0,0 +1,142 @@
+"""
+Sparse positive definite matrix together with a Cholesky factorization object.
+"""
+const PDSparseMat{T<:Real,S<:AbstractSparseMatrix{T},C<:CholTypeSparse} = PDMat{T,S,C}
+
+function PDMats.PDMat(mat::AbstractSparseMatrix, chol::CholTypeSparse)
+    PDMat{eltype(mat),typeof(mat),typeof(chol)}(mat, chol)
+end
+Base.@deprecate PDMat{T,S}(d::Int, m::AbstractSparseMatrix{T}, c::CholTypeSparse) where {T,S} PDSparseMat{T,S,typeof(c)}(m, c)
+
+PDMats.PDMat(mat::SparseMatrixCSC) = PDMat(mat, cholesky(mat))
+PDMats.PDMat(fac::CholTypeSparse) = PDMat(sparse(fac), fac)
+
+PDMats.AbstractPDMat(A::CholTypeSparse) = PDMat(A)
+
+### Conversion
+function Base.convert(::Type{PDMat{T}}, a::PDSparseMat) where {T<:Real}
+    # CholTypeSparse only supports Float64 and ComplexF64 type parameters!
+    # So there is no point in recomputing `cholesky(mat)` and we just reuse
+    # the existing Cholesky factorization
+    mat = convert(AbstractMatrix{T}, a.mat)
+    return PDMat{T,typeof(mat),typeof(a.fact)}(mat, a.fact)
+end
+
+### Arithmetics
+
+Base.:\(a::PDSparseMat{T}, x::AbstractVecOrMat{T}) where {T<:Real} = convert(Array{T},a.fact \ convert(Array{Float64},x)) #to avoid limitations in sparse factorization library CHOLMOD, see e.g., julia issue #14076
+Base.:/(x::AbstractVecOrMat{T}, a::PDSparseMat{T}) where {T<:Real} = convert(Array{T},convert(Array{Float64},x) / a.fact )
+
+### Algebra
+
+Base.inv(a::PDSparseMat) = PDMat(inv(a.mat))
+LinearAlgebra.cholesky(a::PDSparseMat) = a.fact
+Base.sqrt(A::PDSparseMat) = PDMat(sqrt(Hermitian(Matrix(A))))
+
+### whiten and unwhiten
+
+function PDMats.whiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat)
+    PDMats.@check_argdims axes(r) == axes(x)
+    PDMats.@check_argdims a.dim == size(x, 1)
+    # Can't use `ldiv!` due to missing support in SparseArrays
+    return copyto!(r, PDMats.chol_lower(cholesky(a)) \ x)
+end
+
+function PDMats.unwhiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat)
+    PDMats.@check_argdims axes(r) == axes(x)
+    PDMats.@check_argdims a.dim == size(x, 1)
+    # `*` is not defined for `PtL` factor components,
+    # so we can't use `chol_lower(cholesky(a)) * x`
+    C = cholesky(a)
+    PtL = sparse(C.L)[C.p, :]
+    return copyto!(r, PtL * x)
+end
+
+function PDMats.whiten(a::PDSparseMat, x::AbstractVecOrMat)
+    PDMats.@check_argdims a.dim == size(x, 1)
+    return PDMats.chol_lower(cholesky(a)) \ x
+end
+
+function PDMats.unwhiten(a::PDSparseMat, x::AbstractVecOrMat)
+    PDMats.@check_argdims a.dim == size(x, 1)
+    # `*` is not defined for `PtL` factor components,
+    # so we can't use `chol_lower(cholesky(a)) * x`
+    C = cholesky(a)
+    PtL = sparse(C.L)[C.p, :]
+    return PtL * x
+end
+
+### quadratic forms
+
+function PDMats.quad(a::PDSparseMat, x::AbstractVecOrMat)
+    PDMats.@check_argdims a.dim == size(x, 1)
+    # https://github.com/JuliaLang/julia/commit/2425ae760fb5151c5c7dd0554e87c5fc9e24de73
+    if VERSION < v"1.4.0-DEV.92"
+        z = a.mat * x
+        return x isa AbstractVector ? dot(x, z) : map(dot, eachcol(x), eachcol(z))
+    else
+        return x isa AbstractVector ? dot(x, a.mat, x) : map(Base.Fix1(quad, a), eachcol(x))
+    end
+end
+
+function PDMats.quad!(r::AbstractArray, a::PDSparseMat, x::AbstractMatrix)
+    PDMats.@check_argdims eachindex(r) == axes(x, 2)
+    @inbounds for i in axes(x, 2)
+        xi = view(x, :, i)
+        # https://github.com/JuliaLang/julia/commit/2425ae760fb5151c5c7dd0554e87c5fc9e24de73
+        if VERSION < v"1.4.0-DEV.92"
+            # Can't use `lmul!` with buffer due to missing support in SparseArrays
+            r[i] = dot(xi, a.mat * xi)
+        else
+            r[i] = dot(xi, a.mat, xi)
+        end
+    end
+    return r
+end
+
+function PDMats.invquad(a::PDSparseMat, x::AbstractVecOrMat)
+    PDMats.@check_argdims a.dim == size(x, 1)
+    z = cholesky(a) \ x
+    return x isa AbstractVector ? dot(x, z) : map(dot, eachcol(x), eachcol(z))
+end
+
+function PDMats.invquad!(r::AbstractArray, a::PDSparseMat, x::AbstractMatrix)
+    PDMats.@check_argdims eachindex(r) == axes(x, 2)
+    PDMats.@check_argdims a.dim == size(x, 1)
+    # Can't use `ldiv!` with buffer due to missing support in SparseArrays
+    C = cholesky(a)
+    @inbounds for i in axes(x, 2)
+        xi = view(x, :, i)
+        r[i] = dot(xi, C \ xi)
+    end
+    return r
+end
+
+
+### tri products
+
+function PDMats.X_A_Xt(a::PDSparseMat, x::AbstractMatrix{<:Real})
+    PDMats.@check_argdims a.dim == size(x, 2)
+    z = a.mat * transpose(x)
+    return Symmetric(x * z)
+end
+
+
+function PDMats.Xt_A_X(a::PDSparseMat, x::AbstractMatrix{<:Real})
+    PDMats.@check_argdims a.dim == size(x, 1)
+    z = a.mat * x
+    return Symmetric(transpose(x) * z)
+end
+
+
+function PDMats.X_invA_Xt(a::PDSparseMat, x::AbstractMatrix{<:Real})
+    PDMats.@check_argdims a.dim == size(x, 2)
+    z = cholesky(a) \ collect(transpose(x))
+    return Symmetric(x * z)
+end
+
+function PDMats.Xt_invA_X(a::PDSparseMat, x::AbstractMatrix{<:Real})
+    PDMats.@check_argdims a.dim == size(x, 1)
+    z = cholesky(a) \ x
+    return Symmetric(transpose(x) * z)
+end
diff --git a/src/PDMats.jl b/src/PDMats.jl
index c8a2a1c..0366019 100644
--- a/src/PDMats.jl
+++ b/src/PDMats.jl
@@ -1,6 +1,6 @@
 module PDMats
 
-    using LinearAlgebra, SparseArrays, SuiteSparse
+    using LinearAlgebra
 
     import Base: +, *, \, /, ==, convert, inv, Matrix, kron
 
@@ -8,7 +8,6 @@ module PDMats
         # Types
         AbstractPDMat,
         PDMat,
-        PDSparseMat,
         PDiagMat,
         ScalMat,
 
@@ -35,19 +34,12 @@ module PDMats
     """
     abstract type AbstractPDMat{T<:Real} <: AbstractMatrix{T} end
 
-    const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD)
-
     # source files
 
     include("utils.jl")
     include("chol.jl")
 
     include("pdmat.jl")
-
-    if HAVE_CHOLMOD
-        include("pdsparsemat.jl")
-    end
-
     include("pdiagmat.jl")
     include("scalmat.jl")
 
@@ -56,4 +48,8 @@ module PDMats
 
     include("deprecates.jl")
 
+    # Support for sparse arrays
+    if !isdefined(Base, :get_extension)
+        include("../ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl")
+    end
 end # module
diff --git a/src/addition.jl b/src/addition.jl
index ab4cbc0..05f4b84 100644
--- a/src/addition.jl
+++ b/src/addition.jl
@@ -4,38 +4,18 @@
 +(a::PDMat, b::AbstractPDMat) = PDMat(a.mat + Matrix(b))
 +(a::PDiagMat, b::AbstractPDMat) = PDMat(_adddiag!(Matrix(b), a.diag, true))
 +(a::ScalMat, b::AbstractPDMat) = PDMat(_adddiag!(Matrix(b), a.value))
-if HAVE_CHOLMOD
-    +(a::PDSparseMat, b::AbstractPDMat) = PDMat(a.mat + Matrix(b))
-end
 
 +(a::PDMat, b::PDMat) = PDMat(a.mat + b.mat)
 +(a::PDMat, b::PDiagMat) = PDMat(_adddiag(a.mat, b.diag))
 +(a::PDMat, b::ScalMat) = PDMat(_adddiag(a.mat, b.value))
-if HAVE_CHOLMOD
-    +(a::PDMat, b::PDSparseMat) = PDMat(a.mat + b.mat)
-end
 
 +(a::PDiagMat, b::PDMat) = PDMat(_adddiag(b.mat, a.diag))
 +(a::PDiagMat, b::PDiagMat) = PDiagMat(a.diag + b.diag)
 +(a::PDiagMat, b::ScalMat) = PDiagMat(a.diag .+ b.value)
-if HAVE_CHOLMOD
-    +(a::PDiagMat, b::PDSparseMat) = PDSparseMat(_adddiag(b.mat, a.diag))
-end
 
 +(a::ScalMat, b::PDMat) = PDMat(_adddiag(b.mat, a.value))
 +(a::ScalMat, b::PDiagMat) = PDiagMat(a.value .+ b.diag)
 +(a::ScalMat, b::ScalMat) = ScalMat(a.dim, a.value + b.value)
-if HAVE_CHOLMOD
-    +(a::ScalMat, b::PDSparseMat) = PDSparseMat(_adddiag(b.mat, a.value))
-end
-
-if HAVE_CHOLMOD
-    +(a::PDSparseMat, b::PDMat) = PDMat(Matrix(a) + b.mat)
-    +(a::PDSparseMat, b::PDiagMat) = PDSparseMat(_adddiag(a.mat, b.diag))
-    +(a::PDSparseMat, b::ScalMat) = PDSparseMat(_adddiag(a.mat, b.value))
-    +(a::PDSparseMat, b::PDSparseMat) = PDSparseMat(a.mat + b.mat)
-end
-
 
 # between pdmat and uniformscaling (multiple of identity)
 
@@ -45,34 +25,16 @@ end
 pdadd(a::PDMat, b::AbstractPDMat, c::Real) = PDMat(a.mat + Matrix(b * c))
 pdadd(a::PDiagMat, b::AbstractPDMat, c::Real) = PDMat(_adddiag!(Matrix(b * c), a.diag, one(c)))
 pdadd(a::ScalMat, b::AbstractPDMat, c::Real) = PDMat(_adddiag!(Matrix(b * c), a.value))
-if HAVE_CHOLMOD
-    pdadd(a::PDSparseMat, b::AbstractPDMat, c::Real) = PDMat(a.mat + Matrix(b * c))
-end
 
 pdadd(a::PDMat, b::PDMat, c::Real) = PDMat(a.mat + b.mat * c)
 pdadd(a::PDMat, b::PDiagMat, c::Real) = PDMat(_adddiag(a.mat, b.diag, c))
 pdadd(a::PDMat, b::ScalMat, c::Real) = PDMat(_adddiag(a.mat, b.value * c))
-if HAVE_CHOLMOD
-    pdadd(a::PDMat, b::PDSparseMat, c::Real) = PDMat(a.mat + b.mat * c)
-end
 
 pdadd(a::PDiagMat, b::PDMat, c::Real) = PDMat(_adddiag!(b.mat * c, a.diag, one(c)))
 pdadd(a::PDiagMat, b::PDiagMat, c::Real) = PDiagMat(a.diag + b.diag * c)
 pdadd(a::PDiagMat, b::ScalMat, c::Real) = PDiagMat(a.diag .+ b.value * c)
-if HAVE_CHOLMOD
-    pdadd(a::PDiagMat, b::PDSparseMat, c::Real) = PDSparseMat(_adddiag!(b.mat * c, a.diag, one(c)))
-end
 
 pdadd(a::ScalMat, b::PDMat, c::Real) = PDMat(_adddiag!(b.mat * c, a.value))
 pdadd(a::ScalMat, b::PDiagMat, c::Real) = PDiagMat(a.value .+ b.diag * c)
 pdadd(a::ScalMat, b::ScalMat, c::Real) = ScalMat(a.dim, a.value + b.value * c)
-if HAVE_CHOLMOD
-    pdadd(a::ScalMat, b::PDSparseMat, c::Real) = PDSparseMat(_adddiag!(b.mat * c, a.value))
-end
 
-if HAVE_CHOLMOD
-    pdadd(a::PDSparseMat, b::PDMat, c::Real) = PDMat(a.mat + b.mat * c)
-    pdadd(a::PDSparseMat, b::PDiagMat, c::Real) = PDSparseMat(_adddiag(a.mat, b.diag, c))
-    pdadd(a::PDSparseMat, b::ScalMat, c::Real) = PDSparseMat(_adddiag(a.mat, b.value * c))
-    pdadd(a::PDSparseMat, b::PDSparseMat, c::Real) = PDSparseMat(a.mat + b.mat * c)
-end
diff --git a/src/chol.jl b/src/chol.jl
index ec3ec0c..cb08cf0 100644
--- a/src/chol.jl
+++ b/src/chol.jl
@@ -17,14 +17,6 @@ chol_lower(a::Matrix) = cholesky(Symmetric(a, :L)).L
 # type-restricted to be Real, so they are equivalent.
 chol_upper(a::Matrix) = cholesky(Symmetric(a, :U)).U
 
-if HAVE_CHOLMOD
-    CholTypeSparse{T} = SuiteSparse.CHOLMOD.Factor{T}
-
-    # Take into account pivoting!
-    chol_lower(cf::CholTypeSparse) = cf.PtL
-    chol_upper(cf::CholTypeSparse) = cf.UP
-end
-
 # Interface for `Cholesky`
 
 dim(A::Cholesky) = LinearAlgebra.checksquare(A)
diff --git a/src/deprecates.jl b/src/deprecates.jl
index 385206a..603ce41 100644
--- a/src/deprecates.jl
+++ b/src/deprecates.jl
@@ -17,7 +17,7 @@ using Base: @deprecate
 
 @deprecate dim(a::AbstractMatrix) LinearAlgebra.checksquare(a)
 
-@deprecate PDMat{T,S}(d::Int, m::AbstractMatrix{T}, c::Cholesky{T,S}) where {T,S} PDMat{T,S}(m, c)
+@deprecate PDMat{T,S}(d::Int, m::AbstractMatrix{T}, c::Cholesky{T,S}) where {T,S} PDMat{T,S,Cholesky{T,S}}(m, c)
 
 @deprecate PDiagMat(dim::Int, diag::AbstractVector{<:Real}) PDiagMat(diag)
 @deprecate PDiagMat{T,V}(dim, diag) where {T<:Real, V<:AbstractVector{T}} PDiagMat{T,V}(diag)
diff --git a/src/pdiagmat.jl b/src/pdiagmat.jl
index 70ab2c5..e542f0e 100644
--- a/src/pdiagmat.jl
+++ b/src/pdiagmat.jl
@@ -68,7 +68,7 @@ function /(x::AbstractVecOrMat, a::PDiagMat)
     @check_argdims a.dim == size(x, 2)
     if VERSION < v"1.9-"
         # return matrix for 1-element vectors `x`, consistent with LinearAlgebra < 1.9
-        return reshape(x, Val(2)) ./ permutedims(a.diag) # = (a' \ x')'
+        return reshape(x, Val(2)) ./ permutedims(a.diag) # = (a' \ x')
     else
         return x ./ (x isa AbstractVector ? a.diag : a.diag')
     end
diff --git a/src/pdmat.jl b/src/pdmat.jl
index 7f4fcf0..61e7681 100644
--- a/src/pdmat.jl
+++ b/src/pdmat.jl
@@ -1,19 +1,21 @@
 """
-Full positive definite matrix together with a Cholesky factorization object.
+Full positive definite matrix together with a factorization object.
 """
-struct PDMat{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T}
+struct PDMat{T<:Real,S<:AbstractMatrix{T},F<:Factorization} <: AbstractPDMat{T}
     mat::S
-    chol::Cholesky{T,S}
+    fact::F
 
-    PDMat{T,S}(m::AbstractMatrix{T},c::Cholesky{T,S}) where {T,S} = new{T,S}(m,c)
+    function PDMat{T,S,F}(mat::S, fact::F) where {T,S<:AbstractMatrix{T},F<:Factorization}
+        d = LinearAlgebra.checksquare(mat)
+        if size(fact) != (d, d)
+            throw(DimensionMismatch("Dimensions of the matrix and the factorization are inconsistent."))
+        end
+        new{T,S,F}(mat, fact)
+    end
 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."))
-    end
-    PDMat{T,S}(convert(S, mat), chol)
+function PDMat(mat::AbstractMatrix{T}, chol::Cholesky) where {T<:Real}
+    PDMat{T,typeof(mat),typeof(chol)}(mat, chol)
 end
 
 PDMat(mat::AbstractMatrix) = PDMat(mat, cholesky(mat))
@@ -25,17 +27,19 @@ function Base.getproperty(a::PDMat, s::Symbol)
     end
     return getfield(a, s)
 end
-Base.propertynames(::PDMat) = (:mat, :chol, :dim)
+Base.propertynames(::PDMat) = (:mat, :fact, :dim)
 
 AbstractPDMat(A::Cholesky) = PDMat(A)
 
+### Abstract type alias
+const PDMatCholesky{T<:Real,S<:AbstractMatrix{T}} = PDMat{T,S,<:Cholesky}
+
 ### 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)
+    chol = convert(Factorization{float(T)}, a.fact)
+    mat = convert(AbstractMatrix{T}, a.mat)
+    return PDMat{T,typeof(mat),typeof(chol)}(mat, chol)
 end
 Base.convert(::Type{AbstractPDMat{T}}, a::PDMat) where {T<:Real} = convert(PDMat{T}, a)
 
@@ -44,7 +48,7 @@ Base.convert(::Type{AbstractPDMat{T}}, a::PDMat) where {T<:Real} = convert(PDMat
 Base.size(a::PDMat) = (a.dim, a.dim)
 Base.Matrix{T}(a::PDMat) where {T} = Matrix{T}(a.mat)
 LinearAlgebra.diag(a::PDMat) = diag(a.mat)
-LinearAlgebra.cholesky(a::PDMat) = a.chol
+LinearAlgebra.cholesky(a::PDMatCholesky) = a.fact
 
 ### Work with the underlying matrix in broadcasting
 Base.broadcastable(a::PDMat) = Base.broadcastable(a.mat)
@@ -67,58 +71,57 @@ end
 *(a::PDMat, c::Real) = PDMat(a.mat * c)
 *(a::PDMat, x::AbstractVector) = a.mat * x
 *(a::PDMat, x::AbstractMatrix) = a.mat * x
-\(a::PDMat, x::AbstractVecOrMat) = a.chol \ x
+\(a::PDMat, x::AbstractVecOrMat) = a.fact \ x
 function /(x::AbstractVecOrMat, a::PDMat)
     # /(::AbstractVector, ::Cholesky) is not defined
     if VERSION < v"1.9-"
         # return matrix for 1-element vectors `x`, consistent with LinearAlgebra
-        return reshape(x, Val(2)) / a.chol
+        return reshape(x, Val(2)) / a.fact
     else
-        if x isa AbstractVector
-            return vec(reshape(x, Val(2)) / a.chol)
-        else
-            return x / a.chol
-        end
+        return x isa AbstractVector ? vec(reshape(x, Val(2)) / a.fact) : x / a.fact
     end
 end
 
 ### Algebra
 
-Base.inv(a::PDMat) = PDMat(inv(a.chol))
-LinearAlgebra.det(a::PDMat) = det(a.chol)
-LinearAlgebra.logdet(a::PDMat) = logdet(a.chol)
+Base.inv(a::PDMat) = PDMat(inv(a.fact))
+LinearAlgebra.det(a::PDMat) = det(a.fact)
+LinearAlgebra.logdet(a::PDMat) = logdet(a.fact)
 LinearAlgebra.eigmax(a::PDMat) = eigmax(Symmetric(a.mat))
 LinearAlgebra.eigmin(a::PDMat) = eigmin(Symmetric(a.mat))
-Base.kron(A::PDMat, B::PDMat) = PDMat(kron(A.mat, B.mat), Cholesky(kron(A.chol.U, B.chol.U), 'U', A.chol.info))
-LinearAlgebra.sqrt(A::PDMat) = PDMat(sqrt(Hermitian(A.mat)))
+Base.sqrt(A::PDMat) = PDMat(sqrt(Hermitian(A.mat)))
+
+function Base.kron(A::PDMatCholesky, B::PDMatCholesky)
+    return PDMat(kron(A.mat, B.mat), Cholesky(kron(A.fact.U, B.fact.U), 'U', A.fact.info))
+end
 
 ### (un)whitening
 
-function whiten!(r::AbstractVecOrMat, a::PDMat, x::AbstractVecOrMat)
+function whiten!(r::AbstractVecOrMat, a::PDMatCholesky, x::AbstractVecOrMat)
     @check_argdims axes(r) == axes(x)
     @check_argdims a.dim == size(x, 1)
     v = _rcopy!(r, x)
     return ldiv!(chol_lower(cholesky(a)), v)
 end
-function unwhiten!(r::AbstractVecOrMat, a::PDMat, x::AbstractVecOrMat)
+function unwhiten!(r::AbstractVecOrMat, a::PDMatCholesky, x::AbstractVecOrMat)
     @check_argdims axes(r) == axes(x)
     @check_argdims a.dim == size(x, 1)
     v = _rcopy!(r, x)
     return lmul!(chol_lower(cholesky(a)), v)
 end
 
-function whiten(a::PDMat, x::AbstractVecOrMat)
+function whiten(a::PDMatCholesky, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     return chol_lower(cholesky(a)) \ x
 end
-function unwhiten(a::PDMat, x::AbstractVecOrMat)
+function unwhiten(a::PDMatCholesky, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     return chol_lower(cholesky(a)) * x
 end
 
 ## quad/invquad
 
-function quad(a::PDMat, x::AbstractVecOrMat)
+function quad(a::PDMatCholesky, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     aU_x = chol_upper(cholesky(a)) * x
     if x isa AbstractVector
@@ -128,7 +131,7 @@ function quad(a::PDMat, x::AbstractVecOrMat)
     end
 end
 
-function quad!(r::AbstractArray, a::PDMat, x::AbstractMatrix)
+function quad!(r::AbstractArray, a::PDMatCholesky, x::AbstractMatrix)
     @check_argdims eachindex(r) == axes(x, 2)
     @check_argdims a.dim == size(x, 1)
     aU = chol_upper(cholesky(a))
@@ -141,7 +144,7 @@ function quad!(r::AbstractArray, a::PDMat, x::AbstractMatrix)
     return r
 end
 
-function invquad(a::PDMat, x::AbstractVecOrMat)
+function invquad(a::PDMatCholesky, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     inv_aL_x = chol_lower(cholesky(a)) \ x
     if x isa AbstractVector
@@ -151,7 +154,7 @@ function invquad(a::PDMat, x::AbstractVecOrMat)
     end
 end
 
-function invquad!(r::AbstractArray, a::PDMat, x::AbstractMatrix)
+function invquad!(r::AbstractArray, a::PDMatCholesky, x::AbstractMatrix)
     @check_argdims eachindex(r) == axes(x, 2)
     @check_argdims a.dim == size(x, 1)
     aL = chol_lower(cholesky(a))
@@ -166,41 +169,41 @@ end
 
 ### tri products
 
-function X_A_Xt(a::PDMat, x::AbstractMatrix{<:Real})
+function X_A_Xt(a::PDMatCholesky, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 2)
-    z = x * chol_lower(a.chol)
+    z = x * chol_lower(cholesky(a))
     return Symmetric(z * transpose(z))
 end
 
-function Xt_A_X(a::PDMat, x::AbstractMatrix{<:Real})
+function Xt_A_X(a::PDMatCholesky, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 1)
-    z = chol_upper(a.chol) * x
+    z = chol_upper(cholesky(a)) * x
     return Symmetric(transpose(z) * z)
 end
 
-function X_invA_Xt(a::PDMat, x::AbstractMatrix{<:Real})
+function X_invA_Xt(a::PDMatCholesky, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 2)
-    z = x / chol_upper(a.chol)
+    z = x / chol_upper(cholesky(a))
     return Symmetric(z * transpose(z))
 end
 
-function Xt_invA_X(a::PDMat, x::AbstractMatrix{<:Real})
+function Xt_invA_X(a::PDMatCholesky, x::AbstractMatrix{<:Real})
     @check_argdims a.dim == size(x, 1)
-    z = chol_lower(a.chol) \ x
+    z = chol_lower(cholesky(a)) \ x
     return Symmetric(transpose(z) * z)
 end
 
 ### Specializations for `Array` arguments with reduced allocations
 
-function quad(a::PDMat{<:Real,<:Vector}, x::Matrix)
+function quad(a::PDMatCholesky{T,Matrix{T}}, x::Matrix) where {T<:Real}
     @check_argdims a.dim == size(x, 1)
-    T = typeof(zero(eltype(a)) * abs2(zero(eltype(x))))
-    return quad!(Vector{T}(undef, size(x, 2)), a, x)
+    S = typeof(zero(T) * abs2(zero(eltype(x))))
+    return quad!(Vector{S}(undef, size(x, 2)), a, x)
 end
 
-function invquad(a::PDMat{<:Real,<:Vector}, x::Matrix)
+function invquad(a::PDMatCholesky{T,Matrix{T}}, x::Matrix) where {T<:Real}
     @check_argdims a.dim == size(x, 1)
-    T = typeof(abs2(zero(eltype(x))) / zero(eltype(a)))
-    return invquad!(Vector{T}(undef, size(x, 2)), a, x)
+    S = typeof(abs2(zero(eltype(x))) / zero(T))
+    return invquad!(Vector{S}(undef, size(x, 2)), a, x)
 end
 
diff --git a/src/pdsparsemat.jl b/src/pdsparsemat.jl
deleted file mode 100644
index eb5f640..0000000
--- a/src/pdsparsemat.jl
+++ /dev/null
@@ -1,186 +0,0 @@
-"""
-Sparse positive definite matrix together with a Cholesky factorization object.
-"""
-struct PDSparseMat{T<:Real,S<:AbstractSparseMatrix} <: AbstractPDMat{T}
-    mat::S
-    chol::CholTypeSparse
-
-    PDSparseMat{T,S}(m::AbstractSparseMatrix{T},c::CholTypeSparse) where {T,S} =
-        new{T,S}(m,c) #add {T} to CholTypeSparse argument once #14076 is implemented
-end
-@deprecate PDSparseMat{T,S}(d::Int, m::AbstractSparseMatrix{T}, c::CholTypeSparse) where {T,S} PDSparseMat{T,S}(m, c)
-
-function PDSparseMat(mat::AbstractSparseMatrix,chol::CholTypeSparse)
-    d = LinearAlgebra.checksquare(mat)
-    size(chol, 1) == d ||
-        throw(DimensionMismatch("Dimensions of mat and chol are inconsistent."))
-    PDSparseMat{eltype(mat),typeof(mat)}(mat, chol)
-end
-
-PDSparseMat(mat::SparseMatrixCSC) = PDSparseMat(mat, cholesky(mat))
-PDSparseMat(fac::CholTypeSparse) = PDSparseMat(sparse(fac), fac)
-
-function Base.getproperty(a::PDSparseMat, s::Symbol)
-    if s === :dim
-        return size(getfield(a, :mat), 1)
-    end
-    return getfield(a, s)
-end
-Base.propertynames(::PDSparseMat) = (:mat, :chol, :dim)
-
-AbstractPDMat(A::SparseMatrixCSC) = PDSparseMat(A)
-AbstractPDMat(A::CholTypeSparse) = PDSparseMat(A)
-
-### Conversion
-Base.convert(::Type{PDSparseMat{T}}, a::PDSparseMat{T}) where {T<:Real} = a
-function Base.convert(::Type{PDSparseMat{T}}, a::PDSparseMat) where {T<:Real}
-    # CholTypeSparse only supports Float64 and ComplexF64 type parameters!
-    # So there is no point in recomputing `cholesky(mat)` and we just reuse
-    # the existing Cholesky factorization
-    mat = convert(AbstractMatrix{T}, a.mat)
-    return PDSparseMat{T,typeof(mat)}(mat, a.chol)
-end
-Base.convert(::Type{AbstractPDMat{T}}, a::PDSparseMat) where {T<:Real} = convert(PDSparseMat{T}, a)
-
-### Basics
-
-Base.size(a::PDSparseMat) = (a.dim, a.dim)
-Base.Matrix{T}(a::PDSparseMat) where {T} = Matrix{T}(a.mat)
-LinearAlgebra.diag(a::PDSparseMat) = diag(a.mat)
-LinearAlgebra.cholesky(a::PDSparseMat) = a.chol
-
-### Inheriting from AbstractMatrix
-
-Base.IndexStyle(::Type{PDSparseMat{T,S}}) where {T,S} = IndexStyle(S)
-# Linear Indexing
-Base.@propagate_inbounds Base.getindex(a::PDSparseMat, i::Int) = getindex(a.mat, i)
-# Cartesian Indexing
-Base.@propagate_inbounds Base.getindex(a::PDSparseMat, I::Vararg{Int, 2}) = getindex(a.mat, I...)
-
-### Arithmetics
-
-# add `a * c` to a dense matrix `m` of the same size inplace.
-function pdadd!(r::Matrix, a::Matrix, b::PDSparseMat, c)
-    @check_argdims size(r) == size(a) == size(b)
-    _addscal!(r, a, b.mat, c)
-end
-
-*(a::PDSparseMat, c::Real) = PDSparseMat(a.mat * c)
-*(a::PDSparseMat, x::AbstractMatrix) = a.mat * x  # defining these seperately to avoid
-*(a::PDSparseMat, x::AbstractVector) = a.mat * x  # ambiguity errors
-\(a::PDSparseMat{T}, x::AbstractVecOrMat{T}) where {T<:Real} = convert(Array{T},a.chol \ convert(Array{Float64},x)) #to avoid limitations in sparse factorization library CHOLMOD, see e.g., julia issue #14076
-/(x::AbstractVecOrMat{T}, a::PDSparseMat{T}) where {T<:Real} = convert(Array{T},convert(Array{Float64},x) / a.chol )
-
-### Algebra
-
-Base.inv(a::PDSparseMat{T}) where {T<:Real} = PDMat(inv(a.mat))
-LinearAlgebra.det(a::PDSparseMat) = det(a.chol)
-LinearAlgebra.logdet(a::PDSparseMat) = logdet(a.chol)
-LinearAlgebra.sqrt(A::PDSparseMat) = PDMat(sqrt(Hermitian(Matrix(A))))
-
-### whiten and unwhiten
-
-function whiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat)
-    @check_argdims axes(r) == axes(x)
-    @check_argdims a.dim == size(x, 1)
-    # Can't use `ldiv!` due to missing support in SparseArrays
-    return copyto!(r, chol_lower(a.chol) \ x)
-end
-
-function unwhiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat)
-    @check_argdims axes(r) == axes(x)
-    @check_argdims a.dim == size(x, 1)
-    # `*` is not defined for `PtL` factor components,
-    # so we can't use `chol_lower(a.chol) * x`
-    C = a.chol
-    PtL = sparse(C.L)[C.p, :]
-    return copyto!(r, PtL * x)
-end
-
-function whiten(a::PDSparseMat, x::AbstractVecOrMat)
-    @check_argdims a.dim == size(x, 1)
-    return chol_lower(cholesky(a)) \ x
-end
-
-function unwhiten(a::PDSparseMat, x::AbstractVecOrMat)
-    @check_argdims a.dim == size(x, 1)
-    # `*` is not defined for `PtL` factor components,
-    # so we can't use `chol_lower(a.chol) * x`
-    C = a.chol
-    PtL = sparse(C.L)[C.p, :]
-    return PtL * x
-end
-
-### quadratic forms
-
-function quad(a::PDSparseMat, x::AbstractVecOrMat)
-    @check_argdims a.dim == size(x, 1)
-    # https://github.com/JuliaLang/julia/commit/2425ae760fb5151c5c7dd0554e87c5fc9e24de73
-    if VERSION < v"1.4.0-DEV.92"
-        z = a.mat * x
-        return x isa AbstractVector ? dot(x, z) : map(dot, eachcol(x), eachcol(z))
-    else
-        return x isa AbstractVector ? dot(x, a.mat, x) : map(Base.Fix1(quad, a), eachcol(x))
-    end
-end
-
-function quad!(r::AbstractArray, a::PDSparseMat, x::AbstractMatrix)
-    @check_argdims eachindex(r) == axes(x, 2)
-    @inbounds for i in axes(x, 2)
-        xi = view(x, :, i)
-        # https://github.com/JuliaLang/julia/commit/2425ae760fb5151c5c7dd0554e87c5fc9e24de73
-        if VERSION < v"1.4.0-DEV.92"
-            # Can't use `lmul!` with buffer due to missing support in SparseArrays
-            r[i] = dot(xi, a.mat * xi)
-        else
-            r[i] = dot(xi, a.mat, xi)
-        end
-    end
-    return r
-end
-
-function invquad(a::PDSparseMat, x::AbstractVecOrMat)
-    @check_argdims a.dim == size(x, 1)
-    z = a.chol \ x
-    return x isa AbstractVector ? dot(x, z) : map(dot, eachcol(x), eachcol(z))
-end
-
-function invquad!(r::AbstractArray, a::PDSparseMat, x::AbstractMatrix)
-    @check_argdims eachindex(r) == axes(x, 2)
-    @check_argdims a.dim == size(x, 1)
-    # Can't use `ldiv!` with buffer due to missing support in SparseArrays
-    @inbounds for i in axes(x, 2)
-        xi = view(x, :, i)
-        r[i] = dot(xi, a.chol \ xi)
-    end
-    return r
-end
-
-
-### tri products
-
-function X_A_Xt(a::PDSparseMat, x::AbstractMatrix{<:Real})
-    @check_argdims a.dim == size(x, 2)
-    z = a.mat * transpose(x)
-    return Symmetric(x * z)
-end
-
-
-function Xt_A_X(a::PDSparseMat, x::AbstractMatrix{<:Real})
-    @check_argdims a.dim == size(x, 1)
-    z = a.mat * x
-    return Symmetric(transpose(x) * z)
-end
-
-
-function X_invA_Xt(a::PDSparseMat, x::AbstractMatrix{<:Real})
-    @check_argdims a.dim == size(x, 2)
-    z = a.chol \ collect(transpose(x))
-    return Symmetric(x * z)
-end
-
-function Xt_invA_X(a::PDSparseMat, x::AbstractMatrix{<:Real})
-    @check_argdims a.dim == size(x, 1)
-    z = a.chol \ x
-    return Symmetric(transpose(x) * z)
-end
diff --git a/src/utils.jl b/src/utils.jl
index 398bbcc..fecc3c8 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -10,27 +10,27 @@ end
 _rcopy!(r, x) = (r === x || copyto!(r, x); r)
 
 
-function _addscal!(r::Matrix, a::Matrix, b::Union{Matrix, SparseMatrixCSC}, c::Real)
-    if c == one(c)
-        for i in eachindex(a)
+function _addscal!(r::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, c::Real)
+    if isone(c)
+        for i in eachindex(r, a, b)
             @inbounds r[i] = a[i] + b[i]
         end
     else
-        for i in eachindex(a)
-            @inbounds r[i] = a[i] + b[i] * c
+        for i in eachindex(r, a, b)
+            @inbounds r[i] = muladd(b[i], c, a[i])
         end
     end
     return r
 end
 
-function _adddiag!(a::Union{Matrix, SparseMatrixCSC}, v::Real)
+function _adddiag!(a::AbstractMatrix, v::Real)
     for i in diagind(a)
         @inbounds a[i] += v
     end
     return a
 end
 
-function _adddiag!(a::Union{Matrix, SparseMatrixCSC}, v::AbstractVector, c::Real)
+function _adddiag!(a::AbstractMatrix, v::AbstractVector, c::Real)
     @check_argdims eachindex(v) == axes(a, 1) == axes(a, 2)
     if c == one(c)
         for i in eachindex(v)
@@ -44,9 +44,9 @@ function _adddiag!(a::Union{Matrix, SparseMatrixCSC}, v::AbstractVector, c::Real
     return a
 end
 
-_adddiag(a::Union{Matrix, SparseMatrixCSC}, v::Real) = _adddiag!(copy(a), v)
-_adddiag(a::Union{Matrix, SparseMatrixCSC}, v::AbstractVector, c::Real) = _adddiag!(copy(a), v, c)
-_adddiag(a::Union{Matrix, SparseMatrixCSC}, v::AbstractVector{T}) where {T<:Real} = _adddiag!(copy(a), v, one(T))
+_adddiag(a::AbstractMatrix, v::Real) = _adddiag!(copy(a), v)
+_adddiag(a::AbstractMatrix, v::AbstractVector, c::Real) = _adddiag!(copy(a), v, c)
+_adddiag(a::AbstractMatrix, v::AbstractVector{T}) where {T<:Real} = _adddiag!(copy(a), v, one(T))
 
 
 function wsumsq(w::AbstractVector, a::AbstractVector)
@@ -125,7 +125,3 @@ 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
diff --git a/test/addition.jl b/test/addition.jl
index 32b48e1..6902a5e 100644
--- a/test/addition.jl
+++ b/test/addition.jl
@@ -25,7 +25,7 @@ Base.getindex(a::ScalMat3D, i::Int, j::Int) = i == j ? a.value : zero(a.value)
         pm2 = PDiagMat(V)
         pm3 = PDiagMat(sparse(V))
         pm4 = ScalMat(3, X)
-        pm5 = PDSparseMat(sparse(M))
+        pm5 = PDMat(sparse(M))
         pm6 = ScalMat3D(X)
 
         pmats = Any[pm1, pm2, pm3, pm4, pm5, pm6]
diff --git a/test/chol.jl b/test/chol.jl
index a327af1..c33b610 100644
--- a/test/chol.jl
+++ b/test/chol.jl
@@ -47,7 +47,7 @@ using PDMats: chol_lower, chol_upper
             @test sum(abs2, PDMats.chol_upper(ch_dense)' \ x) ≈ b
 
             # sparse version
-            if PDMats.HAVE_CHOLMOD
+            if HAVE_CHOLMOD
                 ch_sparse = cholesky(Symmetric(sparse(A), uplo))
                 @test sum(abs2, PDMats.chol_lower(ch_sparse) \ x) ≈ b
                 @test sum(abs2, PDMats.chol_upper(ch_sparse)' \ x) ≈ b
diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl
index a8f7530..2d44b80 100644
--- a/test/pdmtypes.jl
+++ b/test/pdmtypes.jl
@@ -1,4 +1,4 @@
-using LinearAlgebra, PDMats, SparseArrays, SuiteSparse
+using LinearAlgebra, PDMats, SparseArrays
 using Test
 
 @testset "pd matrix types" begin
@@ -10,8 +10,10 @@ using Test
             @test @test_deprecated(PDiagMat(d, d)) == PDiagMat(d)
             x = one(T)
             @test @test_deprecated(ScalMat(2, x, x)) == ScalMat(2, x)
-            s = SparseMatrixCSC{T}(I, 2, 2)
-            @test PDSparseMat(s, cholesky(s)).mat == PDSparseMat(s).mat == PDSparseMat(cholesky(s)).mat
+            if HAVE_CHOLMOD
+                s = SparseMatrixCSC{T}(I, 2, 2)
+                @test PDMat(s, cholesky(s)).mat == PDMat(s).mat == PDMat(cholesky(s)).mat
+            end
         end
 
         @testset "test the functionality" begin
@@ -32,8 +34,10 @@ using Test
             @testset "ScalMat" begin
                 test_pdmat(ScalMat(3,X), X*Matrix{T}(I, 3, 3), cmat_eq=true, verbose=1)
             end
-            @testset "PDSparseMat" begin
-                test_pdmat(PDSparseMat(sparse(M)), M,          cmat_eq=true, verbose=1, t_eig=false)
+            if HAVE_CHOLMOD
+                @testset "PDMat from sparse matrix" begin
+                    test_pdmat(PDMat(sparse(M)), M, cmat_eq=true, verbose=1, t_eig=false)
+                end
             end
         end
 
@@ -45,7 +49,7 @@ using Test
             @test @test_deprecated(PDiagMat(2, d)) == @test_deprecated(PDiagMat{T,Vector{T}}(2, d)) == PDiagMat(d)
             if HAVE_CHOLMOD
                 s = SparseMatrixCSC{T}(I, 2, 2)
-                @test @test_deprecated(PDSparseMat{T, typeof(s)}(2, s, cholesky(s))) == PDSparseMat(s)
+                @test @test_deprecated(PDMat{T, typeof(s)}(2, s, cholesky(s))) == PDMat(s)
             end
         end
     end
@@ -65,7 +69,7 @@ using Test
                 @test B == A
                 @test (B === A) === (S === T)
                 @test (B.mat === A.mat) === (S === T)
-                @test (B.chol === A.chol) === (S === T)
+                @test (B.fact === A.fact) === (S === T)
             end
 
             A = PDiagMat(ones(T, 2))
@@ -87,16 +91,16 @@ using Test
             end
 
             if HAVE_CHOLMOD
-                A = PDSparseMat(SparseMatrixCSC{T}(I, 2, 2))
-                for R in (AbstractArray{S}, AbstractMatrix{S}, AbstractPDMat{S}, PDSparseMat{S})
+                A = PDMat(SparseMatrixCSC{T}(I, 2, 2))
+                for R in (AbstractArray{S}, AbstractMatrix{S}, AbstractPDMat{S}, PDMat{S})
                     B = @inferred(convert(R, A))
-                    @test B isa PDSparseMat{S}
+                    @test B isa PDMat{S}
                     @test B == A
                     @test (B === A) === (S === T)
                     @test (B.mat === A.mat) === (S === T)
                     # CholMOD only supports Float64 and ComplexF64 type parameters!
                     # Hence the Cholesky factorization is reused
-                    @test B.chol === A.chol
+                    @test B.fact === A.fact
                 end
             end
         end
@@ -143,9 +147,9 @@ using Test
         end
 
         # right division not defined for CHOLMOD:
-        # `rdiv!(::Matrix{Float64}, ::SuiteSparse.CHOLMOD.Factor{Float64})` not defined
+        # `rdiv!(::Matrix{Float64}, ::SparseArrays.CHOLMOD.Factor{Float64})` not defined
         if !HAVE_CHOLMOD
-            z = x / PDSparseMat(sparse(first(A), 1, 1)) 
+            z = x / PDMat(sparse(first(A), 1, 1)) 
             @test typeof(z) === typeof(y)
             @test size(z) == size(y)
             @test z ≈ y
@@ -170,6 +174,7 @@ using Test
 
         M = @inferred AbstractPDMat(A)
         @test M isa PDMat
+        @test cholesky(M) isa Cholesky
         @test Matrix(M) ≈ A
         Mat32 = @inferred Matrix{Float32}(M)
         @test eltype(Mat32) == Float32
@@ -177,6 +182,7 @@ using Test
 
         M = @inferred AbstractPDMat(cholesky(A))
         @test M isa PDMat
+        @test cholesky(M) isa Cholesky
         @test Matrix(M) ≈ A
         Mat32 = @inferred Matrix{Float32}(M)
         @test Mat32 isa Matrix{Float32}
@@ -198,7 +204,8 @@ using Test
         @test Matrix(M) ≈ Diagonal(A)
 
         M = @inferred AbstractPDMat(sparse(A))
-        @test M isa PDSparseMat
+        @test M isa PDMat
+        @test cholesky(M) isa CHOLMOD.Factor
         @test Matrix(M) ≈ A
         Mat32 = @inferred Matrix{Float32}(M)
         @test Mat32 isa Matrix{Float32}
@@ -210,7 +217,8 @@ using Test
         else
             M = @inferred AbstractPDMat(cholesky(sparse(A)))
         end
-        @test M isa PDSparseMat
+        @test M isa PDMat
+        @test cholesky(M) isa CHOLMOD.Factor
         @test Matrix(M) ≈ A
     end
 
@@ -218,7 +226,7 @@ using Test
         for dim in (1, 5, 10)
             x = rand(dim, dim)
             M = PDMat(Array(Symmetric(x' * x + I)))
-            @test fieldnames(typeof(M)) == (:mat, :chol)
+            @test fieldnames(typeof(M)) == (:mat, :fact)
             @test propertynames(M) == (fieldnames(typeof(M))..., :dim)
             @test getproperty(M, :dim) === dim
             for p in fieldnames(typeof(M))
@@ -242,8 +250,8 @@ using Test
 
             if HAVE_CHOLMOD
                 x = sprand(dim, dim, 0.2)
-                M = PDSparseMat(sparse(Symmetric(x' * x + I)))
-                @test fieldnames(typeof(M)) == (:mat, :chol)
+                M = PDMat(sparse(Symmetric(x' * x + I)))
+                @test fieldnames(typeof(M)) == (:mat, :fact)
                 @test propertynames(M) == (fieldnames(typeof(M))..., :dim)
                 @test getproperty(M, :dim) === dim
                 for p in fieldnames(typeof(M))
@@ -264,8 +272,8 @@ using Test
             x = sprand(10, 10, 0.2)
             A = sparse(Symmetric(x * x' + I))
             C = cholesky(A)
-            @test_throws DimensionMismatch PDSparseMat(A[:, 1:(end - 1)], C)
-            @test_throws DimensionMismatch PDSparseMat(A[1:(end - 1), 1:(end - 1)], C)
+            @test_throws DimensionMismatch PDMat(A[:, 1:(end - 1)], C)
+            @test_throws DimensionMismatch PDMat(A[1:(end - 1), 1:(end - 1)], C)
         end
     end
 
@@ -296,4 +304,15 @@ using Test
         @test_broken C - B isa Diagonal{Float64, Vector{Float64}}
         @test C - B ≈ Matrix(C) - Matrix(B)
     end
+
+    # Ref https://github.com/JuliaStats/PDMats.jl/pull/207
+    # Note: `cholesky(::SymTridiagonal)` requires https://github.com/JuliaLang/julia/pull/44076
+    if VERSION >= v"1.8.0-DEV.1526"
+        @testset "PDMat from SymTridiagonal" begin
+            S = SymTridiagonal(fill(4, 4), fill(1, 3))
+            M = @inferred(PDMat(S))
+            @test M isa PDMat{Int,<:SymTridiagonal,<:Cholesky}
+            @test M == S
+        end
+    end
 end
diff --git a/test/specialarrays.jl b/test/specialarrays.jl
index 9292c2e..d2e1c2f 100644
--- a/test/specialarrays.jl
+++ b/test/specialarrays.jl
@@ -106,7 +106,7 @@ using StaticArrays
         # Full matrix
         A = Symmetric(BandedMatrix(Eye(5), (1, 1)))
         P = PDMat(A)
-        @test P isa PDMat{Float64, <:BandedMatrix{Float64}}
+        @test P isa PDMat{Float64, <:Symmetric{Float64, <:BandedMatrix{Float64}}}
 
         x = rand(5)
         X = rand(2, 5)
diff --git a/test/testutils.jl b/test/testutils.jl
index e3a6e1e..e1abbb2 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -4,12 +4,28 @@
 #       the implementation of a subtype of AbstractPDMat
 #
 
-using PDMats, SuiteSparse, Test, Random
+using PDMats, LinearAlgebra, SparseArrays, Test, Random
 
 Random.seed!(10)
 
-const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD)
-const PDMatType = HAVE_CHOLMOD ? Union{PDMat, PDSparseMat, PDiagMat} : Union{PDMat, PDiagMat}
+if isdefined(Base, :get_extension)
+    const HAVE_CHOLMOD = isdefined(SparseArrays, :CHOLMOD)
+else
+    import SuiteSparse
+    const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD)
+end
+const PDMatCholesky{T<:Real, S<:AbstractMatrix, C<:Cholesky} = PDMat{T, S, C}
+if HAVE_CHOLMOD
+    if isdefined(Base, :get_extension)
+        const CHOLMOD = SparseArrays.CHOLMOD
+    else
+        const CHOLMOD = SuiteSparse.CHOLMOD
+    end
+    const PDSparseMat{T<:Real, S<:AbstractSparseMatrix, C<:CHOLMOD.Factor} = PDMat{T, S, C}
+    const PDMatType = Union{PDMatCholesky, PDSparseMat, PDiagMat, ScalMat}
+else
+    const PDMatType = Union{PDMatCholesky, PDiagMat, ScalMat}
+end
 
 ## driver function
 function test_pdmat(C, Cmat::Matrix;
@@ -127,7 +143,7 @@ function pdtest_diag(C, Cmat::Matrix, cmat_eq::Bool, verbose::Int)
     end
 end
 
-function pdtest_cholesky(C::Union{PDMat, PDiagMat, ScalMat}, Cmat::Matrix, cmat_eq::Bool, verbose::Int)
+function pdtest_cholesky(C::PDMatType, Cmat::Matrix, cmat_eq::Bool, verbose::Int)
     _pdt(verbose, "cholesky")
     if cmat_eq
         @test cholesky(C).U == cholesky(Cmat).U
@@ -146,8 +162,8 @@ end
 if HAVE_CHOLMOD
     function pdtest_cholesky(C::PDSparseMat, Cmat::Matrix, cmat_eq::Bool, verbose::Int)
         _pdt(verbose, "cholesky")
-        # We special case PDSparseMat because we can't perform equality checks on
-        # `SuiteSparse.CHOLMOD.Factor`s and `SuiteSparse.CHOLMOD.FactorComponent`s
+        # We handle this case specially because we can't perform equality checks on
+        # `SparseArrays.CHOLMOD.Factor`s and `SparseArrays.CHOLMOD.FactorComponent`s
         @test diag(cholesky(C)) ≈ diag(cholesky(Cmat).U)
         # NOTE: `==` also doesn't work because `diag(cholesky(C))` will return `Vector{Float64}`
         # even if the inputs are `Float32`s.
@@ -182,8 +198,8 @@ function pdtest_det(C, Cmat::Matrix, verbose::Int)
     @test det(C) ≈ det(Cmat)
 
     # generic fallback in LinearAlgebra performs LU decomposition
-    if C isa Union{PDMat,PDiagMat,ScalMat}
-	@test iszero(@allocated det(C))
+    if C isa Union{PDMatCholesky,PDiagMat,ScalMat}
+        @test iszero(@allocated det(C))
     end
 end
 
@@ -192,8 +208,8 @@ function pdtest_logdet(C, Cmat::Matrix, verbose::Int)
     @test logdet(C) ≈ logdet(Cmat)
 
     # generic fallback in LinearAlgebra performs LU decomposition
-    if C isa Union{PDMat,PDiagMat,ScalMat}
-	@test iszero(@allocated logdet(C))
+    if C isa Union{PDMatCholesky,PDiagMat,ScalMat}
+        @test iszero(@allocated logdet(C))
     end
 end
 
@@ -244,10 +260,10 @@ 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
+    # Right division with Cholesky requires 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)
+    # `rdiv!(::Matrix{Float64}, ::SparseArrays.CHOLMOD.Factor{Float64})`
+    check_rdiv = !(C isa PDMatCholesky && VERSION < v"1.3.0-DEV.562") && !(HAVE_CHOLMOD && C isa PDSparseMat)
     check_rdiv && @test Matrix(X') / C ≈ (C \ X)'
 
     for i = 1:n
@@ -371,7 +387,9 @@ end
 _randPDMat(T, n) = (X = randn(T, n, n); PDMat(X * X' + LinearAlgebra.I))
 _randPDiagMat(T, n) = PDiagMat(rand(T, n))
 _randScalMat(T, n) = ScalMat(n, rand(T))
-_randPDSparseMat(T, n) = (X = T.(sprand(n, 1, 0.5)); PDSparseMat(X * X' + LinearAlgebra.I))
+if HAVE_CHOLMOD
+    _randPDSparseMat(T, n) = (X = sprand(T, n, 1, 0.5); PDMat(X * X' + LinearAlgebra.I))
+end
 
 function _pd_compare(A::AbstractPDMat, B::AbstractPDMat)
     @test size(A) == size(B)