diff --git a/src/PDMats.jl b/src/PDMats.jl
index 0aa0fd1..2fb1405 100644
--- a/src/PDMats.jl
+++ b/src/PDMats.jl
@@ -37,7 +37,7 @@ module PDMats
     """
     The base type for positive definite matrices.
     """
-    abstract type AbstractPDMat{T<:Real} end
+    abstract type AbstractPDMat{T<:Real} <: AbstractMatrix{T} end
 
     const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD)
 
diff --git a/src/pdiagmat.jl b/src/pdiagmat.jl
index f291b5a..651f35d 100644
--- a/src/pdiagmat.jl
+++ b/src/pdiagmat.jl
@@ -28,6 +28,12 @@ Base.Matrix(a::PDiagMat) = Matrix(Diagonal(a.diag))
 LinearAlgebra.diag(a::PDiagMat) = copy(a.diag)
 LinearAlgebra.cholesky(a::PDiagMat) = cholesky(Diagonal(a.diag))
 
+### Inheriting from AbstractMatrix
+
+Base.size(a::PDiagMat) = (a.dim,a.dim)
+Base.getindex(a::PDiagMat{T},i::Integer) where {T} = i%a.dim == (i÷a.dim+1) ? a.diag[i%a.dim] : zero(T)
+Base.getindex(a::PDiagMat{T},i::Integer,j::Integer) where {T} = i == j ? a.diag[i] : zero(T)
+
 ### Arithmetics
 
 function pdadd!(r::Matrix, a::Matrix, b::PDiagMat, c)
diff --git a/src/pdmat.jl b/src/pdmat.jl
index 8205c20..175615e 100644
--- a/src/pdmat.jl
+++ b/src/pdmat.jl
@@ -32,6 +32,11 @@ Base.Matrix(a::PDMat) = copy(a.mat)
 LinearAlgebra.diag(a::PDMat) = diag(a.mat)
 LinearAlgebra.cholesky(a::PDMat) = a.chol
 
+### Inheriting from AbstractMatrix
+
+Base.size(a::PDMat) = size(a.mat)
+Base.getindex(a::PDMat,i::Int) = getindex(a.mat,i)
+Base.getindex(a::PDMat,i::Int,j::Int) = getindex(a.mat,i,j)
 
 ### Arithmetics
 
@@ -41,7 +46,8 @@ function pdadd!(r::Matrix, a::Matrix, b::PDMat, c)
 end
 
 *(a::PDMat{S}, c::T) where {S<:Real, T<:Real} = PDMat(a.mat * c)
-*(a::PDMat, x::AbstractVecOrMat) = a.mat * x
+*(a::PDMat, x::AbstractVector{T}) where {T} = a.mat * x
+*(a::PDMat, x::AbstractMatrix{T}) where {T} = a.mat * x
 \(a::PDMat, x::AbstractVecOrMat) = a.chol \ x
 
 
diff --git a/src/pdsparsemat.jl b/src/pdsparsemat.jl
index cf5be98..17d955b 100644
--- a/src/pdsparsemat.jl
+++ b/src/pdsparsemat.jl
@@ -32,6 +32,11 @@ Base.Matrix(a::PDSparseMat) = Matrix(a.mat)
 LinearAlgebra.diag(a::PDSparseMat) = diag(a.mat)
 LinearAlgebra.cholesky(a::PDSparseMat) = a.chol
 
+### Inheriting from AbstractMatrix
+
+Base.size(a::PDSparseMat) = (a.dim,a.dim)
+Base.getindex(a::PDSparseMat,i::Integer) = getindex(a.mat,i)
+Base.getindex(a::PDSparseMat,i::Integer,j::Integer) = getindex(a.mat,i,j)
 
 ### Arithmetics
 
diff --git a/src/scalmat.jl b/src/scalmat.jl
index a8c49ff..bf3b9b0 100644
--- a/src/scalmat.jl
+++ b/src/scalmat.jl
@@ -20,6 +20,12 @@ Base.Matrix(a::ScalMat) = Matrix(Diagonal(fill(a.value, a.dim)))
 LinearAlgebra.diag(a::ScalMat) = fill(a.value, a.dim)
 LinearAlgebra.cholesky(a::ScalMat) = cholesky(Diagonal(fill(a.value, a.dim)))
 
+### Inheriting from AbstractMatrix
+
+Base.size(a::ScalMat) = (a.dim,a.dim)
+Base.getindex(a::ScalMat{T},i::Integer) where {T} = i%a.dim == (i÷a.dim+1) ? a.value : zero(T)
+Base.getindex(a::ScalMat{T},i::Integer,j::Integer) where {T} = i == j ? a.value : zero(T)
+
 ### Arithmetics
 
 function pdadd!(r::Matrix, a::Matrix, b::ScalMat, c)