Fix opnorm for banded matrices (#1428)

Author: dkarrasch

src/tridiag.jl

@@ -1166,7 +1166,7 @@ function _opnorm1Inf(A::Tridiagonal, p)
     case = p == Inf
     lowerrange, upperrange = case ? (1:length(A.dl)-1, 2:length(A.dl)) : (2:length(A.dl), 1:length(A.dl)-1)
     normfirst, normend = case ? (norm(first(A.d))+norm(first(A.du)), norm(last(A.dl))+norm(last(A.d))) : (norm(first(A.d))+norm(first(A.dl)), norm(last(A.du))+norm(last(A.d)))
-
+    size(A, 1) == 2 && return max(normfirst, normend)
     return max(
                 mapreduce(t -> sum(norm, t),
                     max,
@@ -1181,7 +1181,7 @@ 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))
-
+    size(A, 1) == 2 && return max(normfirst, normend)
     return max(
                 mapreduce(t -> sum(norm, t),
                     max,

test/bidiag.jl

@@ -1144,29 +1144,17 @@ end
 end
 
 @testset "opnorms" begin
-    B = Bidiagonal([1,-2,3,-4], [1,2,3], 'U')
-
-    @test opnorm(B, 1) == opnorm(Matrix(B), 1)
-    @test opnorm(B, 2) ≈ opnorm(Matrix(B), 2)
-    @test opnorm(B, Inf) == opnorm(Matrix(B), Inf)
-
-    B = Bidiagonal([1,-2,3,-4], [1,2,3], 'L')
-
-    @test opnorm(B, 1) == opnorm(Matrix(B), 1)
-    @test opnorm(B, 2) ≈ opnorm(Matrix(B), 2)
-    @test opnorm(B, Inf) == opnorm(Matrix(B), Inf)
-
-    B = Bidiagonal([2], Int[], 'L')
-
-    @test opnorm(B, 1) == opnorm(Matrix(B), 1)
-    @test opnorm(B, 2) ≈ opnorm(Matrix(B), 2)
-    @test opnorm(B, Inf) == opnorm(Matrix(B), Inf)
-
-    B = Bidiagonal([2], Int[], 'U')
-
-    @test opnorm(B, 1) == opnorm(Matrix(B), 1)
-    @test opnorm(B, 2) ≈ opnorm(Matrix(B), 2)
-    @test opnorm(B, Inf) == opnorm(Matrix(B), Inf)
+    for B in (Bidiagonal([1,-2,3,-4], [1,2,3], 'U'),
+                Bidiagonal([1,-2,3,-4], [1,2,3], 'L'),
+                Bidiagonal([2], Int[], 'L'),
+                Bidiagonal([2], Int[], 'U'),
+                Bidiagonal([1,-2], [-4], 'U'),
+                Bidiagonal([1,-2], [-4], 'L')
+            )
+        @test opnorm(B, 1) == opnorm(Matrix(B), 1)
+        @test opnorm(B, 2) ≈ opnorm(Matrix(B), 2)
+        @test opnorm(B, Inf) == opnorm(Matrix(B), Inf)
+    end
 end
 
 @testset "convert to Bidiagonal" begin

test/tridiag.jl

@@ -1067,27 +1067,17 @@ end
 end
 
 @testset "opnorms" begin
-    T = Tridiagonal([1,2,3], [1,-2,3,-4], [1,2,3])
-
-    @test opnorm(T, 1) == opnorm(Matrix(T), 1)
-    @test_skip opnorm(T, 2) ≈ opnorm(Matrix(T), 2) # currently missing
-    @test opnorm(T, Inf) == opnorm(Matrix(T), Inf)
-
-    S = SymTridiagonal([1,-2,3,-4], [1,2,3])
-
-    @test opnorm(S, 1) == opnorm(Matrix(S), 1)
-    @test_skip opnorm(S, 2) ≈ opnorm(Matrix(S), 2) # currently missing
-    @test opnorm(S, Inf) == opnorm(Matrix(S), Inf)
-
-    T = Tridiagonal(Int[], [-5], Int[])
-    @test opnorm(T, 1) == opnorm(Matrix(T), 1)
-    @test_skip opnorm(T, 2) ≈ opnorm(Matrix(T), 2) # currently missing
-    @test opnorm(T, Inf) == opnorm(Matrix(T), Inf)
-
-    S = SymTridiagonal(T)
-    @test opnorm(S, 1) == opnorm(Matrix(S), 1)
-    @test_skip opnorm(S, 2) ≈ opnorm(Matrix(S), 2) # currently missing
-    @test opnorm(S, Inf) == opnorm(Matrix(S), Inf)
+    for T in (Tridiagonal([1,2,3], [1,-2,3,-4], [1,2,3]),
+                SymTridiagonal([1,-2,3,-4], [1,2,3]),
+                Tridiagonal(Int[], [-5], Int[]),
+                SymTridiagonal([-5], Int[]),
+                Tridiagonal([1], [1,-2], [3]),
+                SymTridiagonal([1,-2], [3])
+             )
+        @test opnorm(T, 1) == opnorm(Matrix(T), 1)
+        @test_skip opnorm(T, 2) ≈ opnorm(Matrix(T), 2) # currently missing
+        @test opnorm(T, Inf) == opnorm(Matrix(T), Inf)
+    end
 end
 
 @testset "block-bidiagonal matrix indexing" begin