combining all of the redundant methods and adding support for block diagonal matrices

Author: mcognetta

stdlib/LinearAlgebra/src/bidiag.jl

@@ -214,16 +214,6 @@ function _opnorm1Inf(B::Bidiagonal, p)
     return max(mapreduce(t -> sum(norm, t), max, zip(view(B.dv, (1:length(B.ev)) .+ !case), B.ev)), normdend)
 end
 
-function opnorm(B::Bidiagonal, p::Real=2)
-    if p == 2
-        return opnorm2(B)
-    elseif p == 1 || p == Inf
-        return _opnorm1Inf(B, p)
-    else
-        throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
-    end
-end
-
 ####################
 # Generic routines #
 ####################

stdlib/LinearAlgebra/src/diagonal.jl

@@ -548,12 +548,11 @@ function svd(D::Diagonal{<:Number})
     return SVD(Up, S[piv], copy(Vp'))
 end
 
-_opnorm1Inf(A::Diagonal) = maximum(norm(x) for x in A.diag)
+_opnorm12Inf(A::Diagonal, p) = maximum(opnorm(x, p) for x in A.diag)
+
 function opnorm(A::Diagonal, p::Real=2)
-    if p == 2
-        return opnorm2(A)
-    elseif p == 1 || p == Inf
-        return _opnorm1Inf(A)
+    if p == 1 || p == Inf || p == 2
+        return _opnorm12Inf(A, p)
     else
         throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
     end

stdlib/LinearAlgebra/src/special.jl

@@ -338,3 +338,16 @@ end
 
 ==(A::Bidiagonal, B::SymTridiagonal) = iszero(B.ev) && iszero(A.ev) && A.dv == B.dv
 ==(B::SymTridiagonal, A::Bidiagonal) = A == B
+
+# op norm for structured matrices other than diagonal
+# the _opnorm1Inf methods are housed in the structured type's respective files
+
+function opnorm(A::Union{Bidiagonal,Tridiagonal,SymTridiagonal}, p::Real=2)
+    if p == 2
+        return opnorm2(A)
+    elseif p == 1 || p == Inf
+        return _opnorm1Inf(A, p)
+    else
+        throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
+    end
+end

stdlib/LinearAlgebra/src/tridiag.jl

@@ -656,19 +656,9 @@ function _opnorm1Inf(A::Tridiagonal, p)
                 normfirst, normend)
 end
 
-function opnorm(A::Tridiagonal, p::Real=2)
-    if p == 2
-        return opnorm2(A)
-    elseif p == 1 || p == Inf
-        _opnorm1Inf(A, p)
-    else
-        throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
-    end
-end
-
 # SymTridiagonal
 
-function _opnorm1Inf(A::SymTridiagonal)
+function _opnorm1Inf(A::SymTridiagonal, p::Real)
     size(A, 1) == 1 && return norm(first(A.dv))
     lowerrange, upperrange = 1:length(A.ev)-1, 2:length(A.ev)
     normfirst, normend = norm(first(A.dv))+norm(first(A.ev)), norm(last(A.ev))+norm(last(A.dv))
@@ -680,13 +670,3 @@ function _opnorm1Inf(A::SymTridiagonal)
                 ),
                 normfirst, normend)
 end
-
-function opnorm(A::SymTridiagonal, p::Real=2)
-    if p == 2
-        return opnorm2(A)
-    elseif p == 1 || p == Inf # these are the same for symmetric matrices
-        return _opnorm1Inf(A)
-    else
-        throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
-    end
-end

stdlib/LinearAlgebra/test/diagonal.jl

@@ -511,6 +511,13 @@ end
     @test opnorm(D, 1) == opnorm(Matrix(D), 1)
     @test opnorm(D, 2) ≈ opnorm(Matrix(D), 2)
     @test opnorm(D, Inf) == opnorm(Matrix(D), Inf)
+
+    # block diagonal matrices
+    D = Diagonal([[1 2; 3 4], [5 6; 7 8]])
+    A = [1 2 0 0; 3 4 0 0; 0 0 5 6; 0 0 7 8] # full matrix of D
+    @test opnorm(D, 1) == opnorm(A, 1)
+    @test opnorm(D, 2) == opnorm(A, 2)
+    @test opnorm(D, Inf) == opnorm(A, Inf)
 end
 
 @testset "eigenvalue sorting" begin