diff --git a/Project.toml b/Project.toml
index 40abe87..0ad759a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,10 +6,11 @@ version = "0.1.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-julia="1"
+julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/src/DiffTests.jl b/src/DiffTests.jl
index 31d7bd6..be079d4 100644
--- a/src/DiffTests.jl
+++ b/src/DiffTests.jl
@@ -1,6 +1,7 @@
 module DiffTests
 
-using LinearAlgebra: det, norm, dot, tr
+using LinearAlgebra: det, norm, dot, tr, Diagonal, LowerTriangular, UpperTriangular
+using SparseArrays: sparse
 using Statistics: mean
 
 #=
@@ -249,6 +250,43 @@ end
 const INPLACE_ARRAY_TO_ARRAY_FUNCS = (chebyquad!, brown_almost_linear!, trigonometric!,
                                       mutation_test_1!, mutation_test_2!)
 
+############################
+# f(x::VecOrMat)::VecOrMat #
+############################
+
+diag_matrix(::Type{T}, n::Integer) where T<:Real =
+    Diagonal(LinRange(convert(T, 0.01), convert(T, 100.0), n))
+diag_matrix(x::VecOrMat) = diag_matrix(Float64, size(x, 1))
+
+lehmer_matrix(::Type{T}, n::Integer) where T<:Real =
+    [convert(T, min(i, j)/max(i, j)) for i in 1:n, j in 1:n]
+lehmer_matrix(x::VecOrMat) = lehmer_matrix(Float64, size(x, 1))
+
+test_matrix = lehmer_matrix
+
+# left multiplication by a constant matrix
+diag_lmul(x::VecOrMat) = diag_matrix(x) * x
+dense_lmul(x::VecOrMat) = test_matrix(x) * x
+utriag_lmul(x::VecOrMat) =  UpperTriangular(test_matrix(x)) * x
+ltriag_lmul(x::VecOrMat) =  LowerTriangular(test_matrix(x)) * x
+sparse_lmul(x::VecOrMat) = sparse(test_matrix(x)) * x
+sp_utriag_lmul(x::VecOrMat) = UpperTriangular(sparse(test_matrix(x))) * x
+sp_ltriag_lmul(x::VecOrMat) = LowerTriangular(sparse(test_matrix(x))) * x
+
+# left division by a constant matrix
+diag_ldiv(x::VecOrMat) = diag_matrix(x) \ x
+dense_ldiv(x::VecOrMat) = test_matrix(x) \ x
+utriag_ldiv(x::VecOrMat) =  UpperTriangular(test_matrix(x)) \ x
+ltriag_ldiv(x::VecOrMat) =  LowerTriangular(test_matrix(x)) \ x
+sparse_ldiv(x::VecOrMat) = sparse(test_matrix(x)) \ x
+sp_utriag_ldiv(x::VecOrMat) = UpperTriangular(sparse(test_matrix(x))) \ x
+sp_ltriag_ldiv(x::VecOrMat) = LowerTriangular(sparse(test_matrix(x))) \ x
+
+const VECTOR_TO_VECTOR_FUNCS = (diag_lmul, dense_lmul, utriag_lmul, ltriag_lmul,
+                                sparse_lmul, sp_utriag_lmul, sp_ltriag_lmul,
+                                diag_ldiv, utriag_ldiv, ltriag_ldiv,
+                                sparse_ldiv, sp_utriag_ldiv, sp_ltriag_ldiv,)
+
 ######################
 # f(x::Array)::Array #
 ######################
@@ -274,6 +312,10 @@ const ARRAY_TO_ARRAY_FUNCS = (-, chebyquad, brown_almost_linear, trigonometric,
 # f(::Matrix)::Matrix #
 #######################
 
-const MATRIX_TO_MATRIX_FUNCS = (inv,)
+const MATRIX_TO_MATRIX_FUNCS = (inv,
+                                diag_lmul, dense_lmul, utriag_lmul, ltriag_lmul,
+                                sparse_lmul, sp_utriag_lmul, sp_ltriag_lmul,
+                                diag_ldiv, utriag_ldiv, ltriag_ldiv,
+                                sparse_ldiv, sp_utriag_ldiv, sp_ltriag_ldiv,)
 
 end # module
diff --git a/test/runtests.jl b/test/runtests.jl
index fe5836d..d217acf 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -34,12 +34,16 @@ for f in DiffTests.ARRAY_TO_ARRAY_FUNCS
     @test isa(f(y), Array)
 end
 
+for f in DiffTests.VECTOR_TO_VECTOR_FUNCS
+    @test isa(f(y), Vector)
+end
+
 for f in DiffTests.MATRIX_TO_MATRIX_FUNCS
-    @test isa(f(A), Array)
+    @test isa(f(A), Matrix)
 end
 
 for f in DiffTests.BINARY_MATRIX_TO_MATRIX_FUNCS
-    @test isa(f(A, B), Array)
+    @test isa(f(A, B), Matrix)
 end
 
 # f! returns Nothing