diff --git a/test/diagonal.jl b/test/diagonal.jl
index 2b7cfd4d..b69fe655 100644
--- a/test/diagonal.jl
+++ b/test/diagonal.jl
@@ -24,6 +24,17 @@ using .Main.SizedArrays
 const n=12 # Size of matrix problem to test
 Random.seed!(1)
 
+# wrapper to avoid dispatching to diagonal methods
+struct NotDiagonal{T,A<:AbstractMatrix{T}} <: AbstractMatrix{T}
+    a :: A
+end
+Base.size(N::NotDiagonal) = size(N.a)
+Base.getindex(N::NotDiagonal, i::Int, j::Int) = N.a[i, j]
+LinearAlgebra.isdiag(N::NotDiagonal) = false # this contradicts `getindex`
+LinearAlgebra.ishermitian(N::NotDiagonal) = ishermitian(N.a)
+LinearAlgebra.istriu(N::NotDiagonal) = istriu(N.a)
+LinearAlgebra.istril(N::NotDiagonal) = istril(N.a)
+
 @testset for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
     dd=convert(Vector{elty}, randn(n))
     vv=convert(Vector{elty}, randn(n))
@@ -34,13 +45,18 @@ Random.seed!(1)
         UU+=im*convert(Matrix{elty}, randn(n,n))
     end
     D = Diagonal(dd)
-    DM = Matrix(Diagonal(dd))
+    M = Matrix(D)
+    # we can't directly compare with a Matrix, since the dense methods often dispatch
+    # to Diagonal ones. We therefore compare with other structured matrix types
+    # which have their own implementations.
+    # We wrap the complex matrices in NotDiagonal to avoid falling back to Diagonal methods
+    DM = elty <: Real ? Hermitian(M) : NotDiagonal(UpperTriangular(M))
 
     @testset "constructor" begin
         for x in (dd, GenericArray(dd))
-            @test Diagonal(x)::Diagonal{elty,typeof(x)} == DM
+            @test Diagonal(x)::Diagonal{elty,typeof(x)} == M
             @test Diagonal(x).diag === x
-            @test Diagonal{elty}(x)::Diagonal{elty,typeof(x)} == DM
+            @test Diagonal{elty}(x)::Diagonal{elty,typeof(x)} == M
             @test Diagonal{elty}(x).diag === x
             @test Diagonal{elty}(D) === D
         end
@@ -80,9 +96,9 @@ Random.seed!(1)
         @test typeof(convert(Diagonal{ComplexF32},D)) <: Diagonal{ComplexF32}
         @test typeof(convert(AbstractMatrix{ComplexF32},D)) <: Diagonal{ComplexF32}
 
-        @test Array(real(D)) == real(DM)
-        @test Array(abs.(D)) == abs.(DM)
-        @test Array(imag(D)) == imag(DM)
+        @test Array(real(D)) == real(M)
+        @test Array(abs.(D)) == abs.(M)
+        @test Array(imag(D)) == imag(M)
 
         @test parent(D) == dd
         @test D[1,1] == dd[1]
@@ -224,7 +240,7 @@ Random.seed!(1)
             if elty <: Real
                 @test Array(abs.(D)^a) ≈ abs.(DM)^a
             else
-                @test Array(D^a) ≈ DM^a
+                @test Array(D^a) ≈ DM^a rtol=max(eps(relty), 1e-15) # TODO: improve precision
             end
             @test Diagonal(1:100)^2 == Diagonal((1:100).^2)
             p = 3