Test Diagonal functions against non-dense matrices (#1316)

In the `Diagonal` tests, we often compare results of function
evaluations with that for the corresponding dense matrix. E.g., with `D
= Diagonal(::Vector)` and `DM = Matrix(D)`,
```julia
@test log(Diagonal(abs.(D.diag))) ≈ log(abs.(DM)) atol=n^3*eps(relty)
```
The problem with this is that, in the dense method `log(::Matrix)`, we
check if the matrix is diagonal, and dispatch to the `Diagonal` methods
for performance. This makes this test useless, as we are comparing
identical methods.

In this PR, I've changed the type of `DM` to be `Hermitian` for real
matrices, and `UpperTriangular` for complex ones. Both of these types
have special methods for `log`, at least. The `Hermitian` matrices have
special methods for all the trigonometric functions, and for the complex
matrices, we wrap these in a custom wrapper `NotDiagonal` that avoids
the `isdiag` paths. This wrapper lies about the structure, but it should
be ok for the testing purposes as we ensure that we are comparing
different paths.

Author: jishnub

test/diagonal.jl

@@ -27,6 +27,17 @@ using .Main.ImmutableArrays
 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))
@@ -37,13 +48,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
@@ -83,9 +99,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]
@@ -227,7 +243,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