diff --git a/Project.toml b/Project.toml
index 5e107d6..592b883 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "PDMats"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.9.9"
+version = "0.9.10"
 
 [deps]
 Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
diff --git a/src/generics.jl b/src/generics.jl
index 8bd482b..5dacb00 100644
--- a/src/generics.jl
+++ b/src/generics.jl
@@ -25,6 +25,7 @@ pdadd(a::Matrix{T}, b::AbstractPDMat{S}) where {T<:Real, S<:Real} = pdadd!(simil
 *(a::AbstractPDMat, c::T) where {T<:Real} = a * c
 *(c::T, a::AbstractPDMat) where {T<:Real} = a * c
 /(a::AbstractPDMat, c::T) where {T<:Real} = a * inv(c)
+Base.kron(A::AbstractPDMat, B::AbstractPDMat) = PDMat(kron(Matrix(A), Matrix(B)))
 
 
 ## whiten and unwhiten
@@ -58,7 +59,7 @@ PDMat{Float64,Array{Float64,2}}(2, [4.0 10.0; 10.0 30.0], Cholesky{Float64,Array
 
 julia> W = whiten(a, X)
 2×4 Array{Float64,2}:
-  0.5       0.5       0.5       0.5    
+  0.5       0.5       0.5       0.5
  -0.67082  -0.223607  0.223607  0.67082
 
 julia> W * W'
diff --git a/src/pdiagmat.jl b/src/pdiagmat.jl
index a7f4325..651f35d 100644
--- a/src/pdiagmat.jl
+++ b/src/pdiagmat.jl
@@ -49,7 +49,7 @@ end
 *(a::PDiagMat, c::T) where {T<:Real} = PDiagMat(a.diag * c)
 *(a::PDiagMat, x::StridedVecOrMat) = a.diag .* x
 \(a::PDiagMat, x::StridedVecOrMat) = a.inv_diag .* x
-
+Base.kron(A::PDiagMat, B::PDiagMat) = PDiagMat( vcat([A.diag[i] * B.diag for i in 1:dim(A)]...) )
 
 ### Algebra
 
diff --git a/src/scalmat.jl b/src/scalmat.jl
index cf2221b..bf3b9b0 100644
--- a/src/scalmat.jl
+++ b/src/scalmat.jl
@@ -42,7 +42,7 @@ end
 /(a::ScalMat{T}, c::T) where {T<:Real} = ScalMat(a.dim, a.value / c)
 *(a::ScalMat, x::StridedVecOrMat) = a.value * x
 \(a::ScalMat, x::StridedVecOrMat) = a.inv_value * x
-
+Base.kron(A::ScalMat, B::ScalMat) = ScalMat( dim(A) * dim(B), A.value * B.value )
 
 ### Algebra
 
diff --git a/src/testutils.jl b/src/testutils.jl
index f6d1d6a..745a151 100644
--- a/src/testutils.jl
+++ b/src/testutils.jl
@@ -6,6 +6,8 @@
 
 using Test: @test
 
+const PDMatType = HAVE_CHOLMOD ? Union{PDMat, PDSparseMat, PDiagMat} : Union{PDMat, PDiagMat}
+
 ## driver function
 function test_pdmat(C::AbstractPDMat, Cmat::Matrix;
                     verbose::Int=2,             # the level to display intermediate steps
@@ -30,7 +32,7 @@ function test_pdmat(C::AbstractPDMat, Cmat::Matrix;
     pdtest_cmat(C, Cmat, cmat_eq, verbose)
 
     t_diag && pdtest_diag(C, Cmat, cmat_eq, verbose)
-    isa(C, Union{PDMat, PDSparseMat, PDiagMat}) && t_cholesky && pdtest_cholesky(C, Cmat, cmat_eq, verbose)
+    isa(C, PDMatType) && t_cholesky && pdtest_cholesky(C, Cmat, cmat_eq, verbose)
     t_scale && pdtest_scale(C, Cmat, verbose)
     t_add && pdtest_add(C, Cmat, verbose)
     t_logdet && pdtest_logdet(C, Cmat, verbose)
@@ -107,13 +109,15 @@ function pdtest_cholesky(C::Union{PDMat, PDiagMat}, Cmat::Matrix, cmat_eq::Bool,
     end
 end
 
-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
-    @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.
+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
+        @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.
+    end
 end
 
 function pdtest_scale(C::AbstractPDMat, Cmat::Matrix, verbose::Int)
diff --git a/test/kron.jl b/test/kron.jl
index 401f48d..2655078 100644
--- a/test/kron.jl
+++ b/test/kron.jl
@@ -1,22 +1,34 @@
 using PDMats
 using Test
 
+_randPDMat(T, n) = (X = randn(T, n, n); PDMat(X * X'))
+_randPDiagMat(T, n) = PDiagMat(rand(T, n))
+_randScalMat(T, n) = ScalMat(n, rand(T))
+function _pd_compare(A::AbstractPDMat, B::AbstractPDMat)
+    @test dim(A) == dim(B)
+    @test Matrix(A) ≈ Matrix(B)
+    @test cholesky(A).L ≈ cholesky(B).L
+    @test cholesky(A).U ≈ cholesky(B).U
+    nothing
+end
+function _pd_kron_compare(A::AbstractPDMat, B::AbstractPDMat)
+    PDAkB1 = kron(A, B)
+    PDAkB2 = PDMat( kron( Matrix(A), Matrix(B) ) )
+    _pd_compare(PDAkB1, PDAkB2)
+    nothing
+end
+
 n = 4
 m = 7
 
-for T in [Float64, Float32]
-    X = randn(T, n, n)
-    Y = randn(T, m, m)
-    A = X * X'
-    B = Y * Y'
-    AkB = kron(A, B)
-    PDA = PDMat(A)
-    PDB = PDMat(B)
-    PDAkB1 = PDMat(AkB)
-    PDAkB2 = kron(PDA, PDB)
-    @test PDAkB1.dim == PDAkB2.dim
-    @test PDAkB1.mat ≈ PDAkB2.mat
-    @test PDAkB1.chol.L ≈ PDAkB2.chol.L
-    @test PDAkB1.chol.U ≈ PDAkB2.chol.U
-    @test Matrix(PDAkB2.chol) ≈ PDAkB2.mat
+for T in [Float32, Float64]
+    _pd_kron_compare( _randPDMat(T, n),    _randPDMat(T, m) )
+    _pd_kron_compare( _randPDiagMat(T, n), _randPDiagMat(T, m) )
+    _pd_kron_compare( _randScalMat(T, n),  _randScalMat(T, m) )
+    _pd_kron_compare( _randPDMat(T, n),    _randPDiagMat(T, m) )
+    _pd_kron_compare( _randPDiagMat(T, m), _randPDMat(T, n) )
+    _pd_kron_compare( _randPDMat(T, n),    _randScalMat(T, m) )
+    _pd_kron_compare( _randScalMat(T, m),  _randPDMat(T, n) )
+    _pd_kron_compare( _randPDiagMat(T, n), _randScalMat(T, m) )
+    _pd_kron_compare( _randScalMat(T, m),  _randPDiagMat(T, n) )
 end