Specialize norm-related functions for Diagonal (#1457)

For banded matrices such as `Diagonal`, we may skip the non-stored
elements in computing the norm, as the contribution of these would be
zero.

This improves performance.
On master
```julia
julia> D1 = Diagonal(rand(1000));

julia> D2 = Diagonal(rand(1000));

julia> @Btime norm($D1);
  2.668 ms (0 allocations: 0 bytes)

julia> @Btime isapprox($D1, $D2);
  8.007 ms (3 allocations: 7.88 KiB)
```
vs this PR
```julia
julia> @Btime norm($D1);
  247.196 ns (0 allocations: 0 bytes)

julia> @Btime isapprox($D1, $D2);
  1.696 μs (0 allocations: 0 bytes)
```

---------

Co-authored-by: Steven G. Johnson <[email protected]>

Author: jishnub

src/diagonal.jl

@@ -1111,6 +1111,23 @@ end
 /(u::AdjointAbsVec, D::Diagonal) = (D' \ u')'
 /(u::TransposeAbsVec, D::Diagonal) = transpose(transpose(D) \ transpose(u))
 
+# norm
+function generic_normMinusInf(D::Diagonal)
+    norm_diag = norm(D.diag, -Inf)
+    return size(D,1) > 1 ? min(norm_diag, zero(norm_diag)) : norm_diag
+end
+generic_normInf(D::Diagonal) = norm(D.diag, Inf)
+generic_norm1(D::Diagonal) = norm(D.diag, 1)
+generic_norm2(D::Diagonal) = norm(D.diag)
+function generic_normp(D::Diagonal, p)
+    v = norm(D.diag, p)
+    if size(D,1) > 1 && p < 0
+        v = norm(zero(v), p)
+    end
+    return v
+end
+norm_x_minus_y(D1::Diagonal, D2::Diagonal) = norm_x_minus_y(D1.diag, D2.diag)
+
 _opnorm1(A::Diagonal) = maximum(norm(x) for x in A.diag)
 _opnormInf(A::Diagonal) = maximum(norm(x) for x in A.diag)
 _opnorm12Inf(A::Diagonal, p) = maximum(opnorm(x, p) for x in A.diag)

test/diagonal.jl

@@ -1576,4 +1576,23 @@ end
     @test D == D2
 end
 
+@testset "norm" begin
+    # test x ≈ y but also ensure that the types are identical
+    function test_isapprox_and_type(x::T, y::T) where {T}
+        @test x ≈ y
+    end
+    @testset "size(D,1) = $(size(D,1))" for D in ( Diagonal(1:3), Diagonal(1:1), Diagonal(1:0) )
+        A = Array(D)
+        @testset for p in -2:2
+            p == 0 && continue
+            s = sum(float.(abs.(D)).^p)^(1/p)
+            test_isapprox_and_type(norm(D, p), isempty(D) ? zero(s) : s)
+            test_isapprox_and_type(norm(D, p), norm(A, p))
+        end
+        test_isapprox_and_type(norm(D, Inf), norm(A, Inf))
+        test_isapprox_and_type(norm(D, -Inf), norm(A, -Inf))
+        test_isapprox_and_type(norm(D, 0), norm(A, 0))
+    end
+end
+
 end # module TestDiagonal