Merge branch 'myopnorm' into fixed-opnorm

Author: dgleich

stdlib/LinearAlgebra/src/bidiag.jl

@@ -248,6 +248,18 @@ function svd(M::Bidiagonal; kw...)
     svd!(copy(M), kw...)
 end
 
+###################
+#     opnorm      #
+###################
+
+function _opnorm1Inf(B::Bidiagonal, p)
+    size(B, 1) == 1 && return norm(first(B.dv))
+    case = xor(p == 1, istriu(B))
+    normd1, normdend = norm(first(B.dv)), norm(last(B.dv))
+    normd1, normdend = case ? (zero(normd1), normdend) : (normd1, zero(normdend))
+    return max(mapreduce(t -> sum(norm, t), max, zip(view(B.dv, (1:length(B.ev)) .+ !case), B.ev)), normdend)
+end
+
 ####################
 # Generic routines #
 ####################

stdlib/LinearAlgebra/src/diagonal.jl

@@ -1036,6 +1036,16 @@ end
 *(A::AbstractTriangular, D::Diagonal) = @invoke *(A::AbstractMatrix, D::Diagonal)
 *(D::Diagonal, A::AbstractTriangular) = @invoke *(D::Diagonal, A::AbstractMatrix)
 
+_opnorm12Inf(A::Diagonal, p) = maximum(opnorm(x, p) for x in A.diag)
+
+function opnorm(A::Diagonal, p::Real=2)
+    if p == 1 || p == Inf || p == 2
+        return _opnorm12Inf(A, p)
+    else
+        throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
+    end
+end 
+
 dot(x::AbstractVector, D::Diagonal, y::AbstractVector) = _mapreduce_prod(dot, x, D, y)
 
 dot(A::Diagonal, B::Diagonal) = dot(A.diag, B.diag)

stdlib/LinearAlgebra/src/special.jl

@@ -466,6 +466,19 @@ end
 ==(A::Bidiagonal, B::SymTridiagonal) = iszero(_evview(B)) && iszero(A.ev) && A.dv == _diagview(B)
 ==(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
+
 # TODO: remove these deprecations (used by SparseArrays in the past)
 const _DenseConcatGroup = Union{}
 const _SpecialArrays = Union{}

stdlib/LinearAlgebra/src/tridiag.jl

@@ -886,6 +886,41 @@ function SymTridiagonal{T}(M::Tridiagonal) where T
     end
 end
 
+###################
+#     opnorms     #
+###################
+
+# Tridiagonal
+
+function _opnorm1Inf(A::Tridiagonal, p)
+    size(A, 1) == 1 && return norm(first(A.d))
+    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)))
+
+    return max(
+                mapreduce(t -> sum(norm, t),
+                    max,
+                    zip(view(A.d, (2:length(A.d)-1)), view(A.dl, lowerrange), view(A.du, upperrange))
+                ),
+                normfirst, normend)
+end
+
+# 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))
+
+    return max(
+                mapreduce(t -> sum(norm, t),
+                    max,
+                    zip(view(A.dv, (2:length(A.dv)-1)), view(A.ev, lowerrange), view(A.ev, upperrange))
+                ),
+                normfirst, normend)
+end
+
 Base._sum(A::Tridiagonal, ::Colon) = sum(A.d) + sum(A.dl) + sum(A.du)
 function Base._sum(A::SymTridiagonal, ::Colon)
     se = sum(_evview(A))

stdlib/LinearAlgebra/test/bidiag.jl

@@ -1138,4 +1138,30 @@ end
     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)
+end
+
 end # module TestBidiagonal

stdlib/LinearAlgebra/test/diagonal.jl

@@ -989,6 +989,26 @@ end
     end
 end
 
+@testset "opnorms" begin
+    D = Diagonal([1,-2,3,-4])
+
+    @test opnorm(D, 1) == opnorm(Matrix(D), 1)
+    @test opnorm(D, 2) ≈ opnorm(Matrix(D), 2)
+    @test opnorm(D, Inf) == opnorm(Matrix(D), Inf)
+
+    D = Diagonal([-11])
+    @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 "(Sym)Tridiagonal division by Diagonal" begin
     for K in (5, 1), elty in (Float64, ComplexF32), overlength in (1, 0)
         S = SymTridiagonal(randn(elty, K), randn(elty, K-overlength))

stdlib/LinearAlgebra/test/tridiag.jl

@@ -602,6 +602,30 @@ end
     @test A90\b90 ≈ inv(A90)*b90
 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)
+end
+
 @testset "singular values of SymTridiag" begin
     @test svdvals(SymTridiagonal([-4,2,3], [0,0])) ≈ [4,3,2]
     @test svdvals(SymTridiagonal(collect(0.:10.), zeros(10))) ≈ reverse(0:10)