diff --git a/Project.toml b/Project.toml
index d9b4481..f708168 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,12 +8,17 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [compat]
+BandedMatrices = "1"
+Random = "1"
+StaticArrays = "1"
+Test = "1"
 julia = "1"
 
 [extras]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BandedMatrices", "StaticArrays", "Test"]
+test = ["BandedMatrices", "StaticArrays", "Random", "Test"]
diff --git a/src/pdiagmat.jl b/src/pdiagmat.jl
index cb43a70..fad2b7e 100644
--- a/src/pdiagmat.jl
+++ b/src/pdiagmat.jl
@@ -29,6 +29,7 @@ Base.convert(::Type{AbstractPDMat{T}}, a::PDiagMat) where {T<:Real} = convert(PD
 
 Base.size(a::PDiagMat) = (a.dim, a.dim)
 Base.Matrix(a::PDiagMat) = Matrix(Diagonal(a.diag))
+Base.Matrix{T}(a::PDiagMat) where {T} = Matrix{T}(Diagonal(a.diag))
 LinearAlgebra.diag(a::PDiagMat) = copy(a.diag)
 LinearAlgebra.cholesky(a::PDiagMat) = Cholesky(Diagonal(map(sqrt, a.diag)), 'U', 0)
 
@@ -37,11 +38,7 @@ Base.broadcastable(a::PDiagMat) = Base.broadcastable(Diagonal(a.diag))
 
 ### Inheriting from AbstractMatrix
 
-function Base.getindex(a::PDiagMat, i::Integer)
-    ncol, nrow = fldmod1(i, a.dim)
-    ncol == nrow ? a.diag[nrow] : zero(eltype(a))
-end
-Base.getindex(a::PDiagMat{T},i::Integer,j::Integer) where {T} = i == j ? a.diag[i] : zero(T)
+Base.@propagate_inbounds Base.getindex(a::PDiagMat{T}, i::Int, j::Int) where {T} = i == j ? a.diag[i] : zero(T)
 
 ### Arithmetics
 
diff --git a/src/pdmat.jl b/src/pdmat.jl
index abaec69..fc3dca4 100644
--- a/src/pdmat.jl
+++ b/src/pdmat.jl
@@ -43,6 +43,7 @@ Base.convert(::Type{AbstractPDMat{T}}, a::PDMat) where {T<:Real} = convert(PDMat
 
 Base.size(a::PDMat) = (a.dim, a.dim)
 Base.Matrix(a::PDMat) = Matrix(a.mat)
+Base.Matrix{T}(a::PDMat) where {T} = Matrix{T}(a.mat)
 LinearAlgebra.diag(a::PDMat) = diag(a.mat)
 LinearAlgebra.cholesky(a::PDMat) = a.chol
 
@@ -51,8 +52,7 @@ Base.broadcastable(a::PDMat) = Base.broadcastable(a.mat)
 
 ### Inheriting from AbstractMatrix
 
-Base.getindex(a::PDMat, i::Int) = getindex(a.mat, i)
-Base.getindex(a::PDMat, I::Vararg{Int, N}) where {N} = getindex(a.mat, I...)
+Base.@propagate_inbounds Base.getindex(a::PDMat, I::Vararg{Int, 2}) = getindex(a.mat, I...)
 
 ### Arithmetics
 
diff --git a/src/pdsparsemat.jl b/src/pdsparsemat.jl
index 032e0ed..b1b602e 100644
--- a/src/pdsparsemat.jl
+++ b/src/pdsparsemat.jl
@@ -46,13 +46,13 @@ Base.convert(::Type{AbstractPDMat{T}}, a::PDSparseMat) where {T<:Real} = convert
 
 Base.size(a::PDSparseMat) = (a.dim, a.dim)
 Base.Matrix(a::PDSparseMat) = Matrix(a.mat)
+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.getindex(a::PDSparseMat,i::Integer) = getindex(a.mat, i)
-Base.getindex(a::PDSparseMat,I::Vararg{Int, N}) where {N} = getindex(a.mat, I...)
+Base.getindex(a::PDSparseMat, I::Vararg{Int, 2}) = getindex(a.mat, I...)
 
 ### Arithmetics
 
diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl
index c749cf8..99ce6ce 100644
--- a/test/pdmtypes.jl
+++ b/test/pdmtypes.jl
@@ -166,19 +166,28 @@ using Test
 
     @testset "AbstractPDMat constructors (#136)" begin
         x = rand(10, 10)
-        A = x' * x + I
+        A = Array(Hermitian(x' * x + I))
 
         M = @inferred AbstractPDMat(A)
         @test M isa PDMat
         @test Matrix(M) ≈ A
+        Mat32 = Matrix{Float32}(M)
+        @test eltype(Mat32) == Float32
+        @test Mat32 ≈ Float32.(A)
 
         M = @inferred AbstractPDMat(cholesky(A))
         @test M isa PDMat
         @test Matrix(M) ≈ A
+        Mat32 = Matrix{Float32}(M)
+        @test eltype(Mat32) == Float32
+        @test Mat32 ≈ Float32.(A)
 
         M = @inferred AbstractPDMat(Diagonal(A))
         @test M isa PDiagMat
         @test Matrix(M) ≈ Diagonal(A)
+        Mat32 = Matrix{Float32}(M)
+        @test eltype(Mat32) == Float32
+        @test Mat32 ≈ Float32.(Diagonal(A))
 
         M = @inferred AbstractPDMat(Symmetric(Diagonal(A)))
         @test M isa PDiagMat
@@ -191,6 +200,9 @@ using Test
         M = @inferred AbstractPDMat(sparse(A))
         @test M isa PDSparseMat
         @test Matrix(M) ≈ A
+        Mat32 = Matrix{Float32}(M)
+        @test eltype(Mat32) == Float32
+        @test Mat32 ≈ Float32.(A)
 
         if VERSION < v"1.6"
             # inference fails e.g. on Julia 1.0
@@ -205,7 +217,7 @@ using Test
     @testset "properties and fields" begin
         for dim in (1, 5, 10)
             x = rand(dim, dim)
-            M = PDMat(x' * x + I)
+            M = PDMat(Array(Hermitian(x' * x + I)))
             @test fieldnames(typeof(M)) == (:mat, :chol)
             @test propertynames(M) == (fieldnames(typeof(M))..., :dim)
             @test getproperty(M, :dim) === dim
@@ -230,7 +242,7 @@ using Test
 
             if HAVE_CHOLMOD
                 x = sprand(dim, dim, 0.2)
-                M = PDSparseMat(x' * x + I)
+                M = PDSparseMat(sparse(Hermitian(x' * x + I)))
                 @test fieldnames(typeof(M)) == (:mat, :chol)
                 @test propertynames(M) == (fieldnames(typeof(M))..., :dim)
                 @test getproperty(M, :dim) === dim
@@ -243,14 +255,14 @@ using Test
 
     @testset "Incorrect dimensions" begin
         x = rand(10, 10)
-        A = x * x' + I
+        A = Array(Hermitian(x * x' + I))
         C = cholesky(A)
         @test_throws DimensionMismatch PDMat(A[:, 1:(end - 1)], C)
         @test_throws DimensionMismatch PDMat(A[1:(end - 1), 1:(end - 1)], C)
 
         if HAVE_CHOLMOD
             x = sprand(10, 10, 0.2)
-            A = x * x' + I
+            A = sparse(Hermitian(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)
@@ -261,7 +273,7 @@ using Test
         # This falls back to the generic method in Julia based on broadcasting
         dim = 4
         x = rand(dim, dim)
-        A = PDMat(x' * x + I)
+        A = PDMat(Array(Hermitian(x' * x + I)))
         @test Base.broadcastable(A) == A.mat
 
         B = PDiagMat(rand(dim))
diff --git a/test/specialarrays.jl b/test/specialarrays.jl
index be6cbe8..8e67d77 100644
--- a/test/specialarrays.jl
+++ b/test/specialarrays.jl
@@ -4,7 +4,7 @@ using StaticArrays
 @testset "Special matrix types" begin
     @testset "StaticArrays" begin
         # Full matrix
-        S = (x -> x * x' + I)(@SMatrix(randn(4, 7)))
+        S = (x -> SMatrix{4,4}(Hermitian(x * x' + I)))(@SMatrix(randn(4, 7)))
         PDS = PDMat(S)
         @test PDS isa PDMat{Float64, <:SMatrix{4, 4, Float64}}
         @test isbits(PDS)
diff --git a/test/testutils.jl b/test/testutils.jl
index 32bd5a9..e3a6e1e 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -4,7 +4,9 @@
 #       the implementation of a subtype of AbstractPDMat
 #
 
-using PDMats, SuiteSparse, Test
+using PDMats, SuiteSparse, Test, Random
+
+Random.seed!(10)
 
 const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD)
 const PDMatType = HAVE_CHOLMOD ? Union{PDMat, PDSparseMat, PDiagMat} : Union{PDMat, PDiagMat}