diff --git a/src/PDMats.jl b/src/PDMats.jl
index 2c450e5..0aa0fd1 100644
--- a/src/PDMats.jl
+++ b/src/PDMats.jl
@@ -2,7 +2,7 @@ module PDMats
 
     using LinearAlgebra, SparseArrays, SuiteSparse
 
-    import Base: +, *, \, /, ==, convert, inv, Matrix
+    import Base: +, *, \, /, ==, convert, inv, Matrix, kron
 
     export
         # Types
diff --git a/src/pdmat.jl b/src/pdmat.jl
index 28a6ec2..158a5e6 100644
--- a/src/pdmat.jl
+++ b/src/pdmat.jl
@@ -49,7 +49,7 @@ Base.inv(a::PDMat) = PDMat(inv(a.chol))
 LinearAlgebra.logdet(a::PDMat) = logdet(a.chol)
 LinearAlgebra.eigmax(a::PDMat) = eigmax(a.mat)
 LinearAlgebra.eigmin(a::PDMat) = eigmin(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))
 
 ### whiten and unwhiten
 
diff --git a/test/kron.jl b/test/kron.jl
new file mode 100644
index 0000000..401f48d
--- /dev/null
+++ b/test/kron.jl
@@ -0,0 +1,22 @@
+using PDMats
+using Test
+
+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
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index cda3a0d..8366a7d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,4 @@
-tests = ["pdmtypes", "addition", "generics"]
+tests = ["pdmtypes", "addition", "generics", "kron"]
 println("Running tests ...")
 
 for t in tests