diff --git a/test/symmetric.jl b/test/symmetric.jl
index 739d983e..50c72e04 100644
--- a/test/symmetric.jl
+++ b/test/symmetric.jl
@@ -876,21 +876,67 @@ end
     end
 end
 
-@testset "Multiplications symmetric/hermitian for $T and $S" for T in
-        (Float16, Float32, Float64, BigFloat), S in (ComplexF16, ComplexF32, ComplexF64)
-    let A = transpose(Symmetric(rand(S, 3, 3))), Bv = Vector(rand(T, 3)), Bm = Matrix(rand(T, 3,3))
+@testset "Multiplications symmetric/hermitian for T=$T and S=$S for size n=$n" for T in
+        (Float16, Float32, Float64, BigFloat, Quaternion{Float64}),
+        S in (T <: Quaternion ? (Quaternion{Float64},) : (ComplexF16, ComplexF32, ComplexF64, Quaternion{Float64})),
+        n in (2, 3, 4)
+    let A = transpose(Symmetric(rand(S, n, n))), Bv = Vector(rand(T, n)), Bm = Matrix(rand(T, n,n))
         @test A * Bv ≈ Matrix(A) * Bv
         @test A * Bm ≈ Matrix(A) * Bm
+        @test A * transpose(Bm) ≈ Matrix(A) * transpose(Bm)
+        @test A * adjoint(Bm) ≈ Matrix(A) * adjoint(Bm)
         @test Bm * A ≈ Bm * Matrix(A)
+        @test transpose(Bm) * A ≈ transpose(Bm) * Matrix(A)
+        @test adjoint(Bm) * A ≈ adjoint(Bm) * Matrix(A)
+        C = similar(Bm, promote_type(T, S))
+        @test mul!(C, A, Bm) ≈ A * Bm
+        @test mul!(adjoint(C), A, adjoint(Bm)) ≈ A * adjoint(Bm)
+        @test mul!(transpose(C), A, transpose(Bm)) ≈ A * transpose(Bm)
+        rand!(C)
+        @test mul!(copy(C), A, Bm, 2, 3) ≈ A * Bm * 2 + C * 3
+        @test mul!(copy(C), Bm, A, 2, 3) ≈ Bm * A * 2 + C * 3
+        @test mul!(adjoint(copy(C)), A, adjoint(Bm), 2, 3) ≈ A * adjoint(Bm) * 2 + adjoint(C) * 3
+        @test mul!(adjoint(copy(C)), adjoint(Bm), A, 2, 3) ≈ adjoint(Bm) * A * 2 + adjoint(C) * 3
+        @test mul!(transpose(copy(C)), A, transpose(Bm), 2, 3) ≈ A * transpose(Bm) * 2 + transpose(C) * 3
+        @test mul!(transpose(copy(C)), transpose(Bm), A, 2, 3) ≈ transpose(Bm) * A * 2 + transpose(C) * 3
+        if eltype(C) <: Complex
+            alpha, beta  = 4+2im, 3+im
+            @test mul!(adjoint(copy(C)), A, adjoint(Bm), alpha, beta) ≈ A * adjoint(Bm) * alpha + adjoint(C) * beta
+            @test mul!(adjoint(copy(C)), adjoint(Bm), A, alpha, beta) ≈ adjoint(Bm) * A * alpha + adjoint(C) * beta
+            @test mul!(transpose(copy(C)), A, transpose(Bm), alpha, beta) ≈ A * transpose(Bm) * alpha + transpose(C) * beta
+            @test mul!(transpose(copy(C)), transpose(Bm), A, alpha, beta) ≈ transpose(Bm) * A * alpha + transpose(C) * beta
+        end
     end
-    let A = adjoint(Hermitian(rand(S, 3,3))), Bv = Vector(rand(T, 3)), Bm = Matrix(rand(T, 3,3))
+    let A = adjoint(Hermitian(rand(S, n,n))), Bv = Vector(rand(T, n)), Bm = Matrix(rand(T, n,n))
         @test A * Bv ≈ Matrix(A) * Bv
         @test A * Bm ≈ Matrix(A) * Bm
+        @test A * transpose(Bm) ≈ Matrix(A) * transpose(Bm)
+        @test A * adjoint(Bm) ≈ Matrix(A) * adjoint(Bm)
         @test Bm * A ≈ Bm * Matrix(A)
+        @test transpose(Bm) * A ≈ transpose(Bm) * Matrix(A)
+        @test adjoint(Bm) * A ≈ adjoint(Bm) * Matrix(A)
+        C = similar(Bm, promote_type(T, S))
+        @test mul!(C, A, Bm) ≈ A * Bm
+        @test mul!(adjoint(C), A, adjoint(Bm)) ≈ A * adjoint(Bm)
+        @test mul!(transpose(C), A, transpose(Bm)) ≈ A * transpose(Bm)
+        rand!(C)
+        @test mul!(copy(C), A, Bm, 2, 3) ≈ A * Bm * 2 + C * 3
+        @test mul!(copy(C), Bm, A, 2, 3) ≈ Bm * A * 2 + C * 3
+        @test mul!(adjoint(copy(C)), A, adjoint(Bm), 2, 3) ≈ A * adjoint(Bm) * 2 + adjoint(C) * 3
+        @test mul!(adjoint(copy(C)), adjoint(Bm), A, 2, 3) ≈ adjoint(Bm) * A * 2 + adjoint(C) * 3
+        @test mul!(transpose(copy(C)), A, transpose(Bm), 2, 3) ≈ A * transpose(Bm) * 2 + transpose(C) * 3
+        @test mul!(transpose(copy(C)), transpose(Bm), A, 2, 3) ≈ transpose(Bm) * A * 2 + transpose(C) * 3
+        if eltype(C) <: Complex
+            alpha, beta  = 4+2im, 3+im
+            @test mul!(adjoint(copy(C)), A, adjoint(Bm), alpha, beta) ≈ A * adjoint(Bm) * alpha + adjoint(C) * beta
+            @test mul!(adjoint(copy(C)), adjoint(Bm), A, alpha, beta) ≈ adjoint(Bm) * A * alpha + adjoint(C) * beta
+            @test mul!(transpose(copy(C)), A, transpose(Bm), alpha, beta) ≈ A * transpose(Bm) * alpha + transpose(C) * beta
+            @test mul!(transpose(copy(C)), transpose(Bm), A, alpha, beta) ≈ transpose(Bm) * A * alpha + transpose(C) * beta
+        end
     end
-    let Ahrs = transpose(Hermitian(Symmetric(rand(T, 3, 3)))),
-        Acs = transpose(Symmetric(rand(S, 3, 3))),
-        Ahcs = transpose(Hermitian(Symmetric(rand(S, 3, 3))))
+    let Ahrs = transpose(Hermitian(Symmetric(rand(T, n, n)))),
+        Acs = transpose(Symmetric(rand(S, n, n))),
+        Ahcs = transpose(Hermitian(Symmetric(rand(S, n, n))))
 
         @test Ahrs * Ahrs ≈ Ahrs * Matrix(Ahrs)
         @test Ahrs * Acs ≈ Ahrs * Matrix(Acs)
@@ -899,9 +945,9 @@ end
         @test Ahrs * Ahcs ≈ Matrix(Ahrs) * Ahcs
         @test Ahcs * Ahrs ≈ Ahcs * Matrix(Ahrs)
     end
-    let Ahrs = adjoint(Hermitian(Symmetric(rand(T, 3, 3)))),
-        Acs = adjoint(Symmetric(rand(S, 3, 3))),
-        Ahcs = adjoint(Hermitian(Symmetric(rand(S, 3, 3))))
+    let Ahrs = adjoint(Hermitian(Symmetric(rand(T, n, n)))),
+        Acs = adjoint(Symmetric(rand(S, n, n))),
+        Ahcs = adjoint(Hermitian(Symmetric(rand(S, n, n))))
 
         @test Ahrs * Ahrs ≈ Ahrs * Matrix(Ahrs)
         @test Ahcs * Ahcs ≈ Matrix(Ahcs) * Matrix(Ahcs)
diff --git a/test/testhelpers/Quaternions.jl b/test/testhelpers/Quaternions.jl
index b1a41426..49539574 100644
--- a/test/testhelpers/Quaternions.jl
+++ b/test/testhelpers/Quaternions.jl
@@ -16,6 +16,7 @@ struct Quaternion{T<:Real} <: Number
 end
 Quaternion{T}(s::Real) where {T<:Real} = Quaternion{T}(T(s), zero(T), zero(T), zero(T))
 Quaternion(s::Real, v1::Real, v2::Real, v3::Real) = Quaternion(promote(s, v1, v2, v3)...)
+Quaternion{T}(q::Quaternion) where {T<:Real} = Quaternion{T}(T(q.s), T(q.v1), T(q.v2), T(q.v3))
 Base.convert(::Type{Quaternion{T}}, s::Real) where {T <: Real} =
     Quaternion{T}(convert(T, s), zero(T), zero(T), zero(T))
 Base.promote_rule(::Type{Quaternion{T}}, ::Type{S}) where {T <: Real, S <: Real} =