Structured opnorm

stdlib/LinearAlgebra/src/bidiag.jl

@@ -202,6 +202,80 @@ function svd(M::Bidiagonal; full::Bool = false)
     svd!(copy(M), full = full)
 end
 
+###################
+#     opnorms     #
+###################
+
+# helpers for upper and lower bidiagonal matrices
+
+function _upopnorm1(A::Bidiagonal{T}) where T
+    n = size(A, 1)
+    Tnorm = typeof(float(real(zero(T))))
+    Tsum = promote_type(Float64, Tnorm)
+    nrm::Tsum = A.dv[1]
+    @inbounds begin
+        for i = 2:n
+            nrmj::Tsum = norm(A.dv[i]) + norm(A.ev[i-1])
+            nrm = max(nrm,nrmj)
+        end
+    end
+    return convert(Tnorm, nrm)
+end
+
+function _loopnorm1(A::Bidiagonal{T}) where T
+    n = size(A, 1)
+    Tnorm = typeof(float(real(zero(T))))
+    Tsum = promote_type(Float64, Tnorm)
+    nrm::Tsum = A.dv[end]
+    @inbounds begin
+        for i = 1:n-1
+            nrmj::Tsum = norm(A.dv[i]) + norm(A.ev[i])
+            nrm = max(nrm,nrmj)
+        end
+    end
+    return convert(Tnorm, nrm)
+end
+
+function _upopnormInf(A::Bidiagonal{T}) where T
+    n = size(A, 1)
+    Tnorm = typeof(float(real(zero(T))))
+    Tsum = promote_type(Float64, Tnorm)
+    nrm::Tsum = A.dv[end]
+    @inbounds begin
+        for i = 1:n-1
+            nrmj::Tsum = norm(A.dv[i]) + norm(A.ev[i])
+            nrm = max(nrm,nrmj)
+        end
+    end
+    return convert(Tnorm, nrm)
+end
+
+function _loopnormInf(A::Bidiagonal{T}) where T
+    n = size(A, 1)
+    Tnorm = typeof(float(real(zero(T))))
+    Tsum = promote_type(Float64, Tnorm)
+    nrm::Tsum = A.dv[1]
+    @inbounds begin
+        for i = 2:n
+            nrmj::Tsum = norm(A.dv[i]) + norm(A.ev[i-1])
+            nrm = max(nrm,nrmj)
+        end
+    end
+    return convert(Tnorm, nrm)
+end
+
+function opnorm(A::Bidiagonal, p::Real=2)
+    if p == 2
+        return opnorm2(A)
+    elseif p == 1
+        A.uplo == 'U' ? _upopnorm1(A) : _loopnorm1(A)
+    elseif p == Inf
+        A.uplo == 'U' ? _upopnormInf(A) : _loopnormInf(A)
+    else
+        throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
+    end
+end
+
 ####################
 # Generic routines #
 ####################

stdlib/LinearAlgebra/src/diagonal.jl

@@ -508,6 +508,20 @@ function svd(D::Diagonal{<:Number})
     return SVD(Up, S[piv], copy(Vp'))
 end
 
+_opnorm1(A::Diagonal) = maximum(norm(x) for x in A.diag)
+_opnormInf(A::Diagonal) = maximum(norm(x) for x in A.diag)
+function opnorm(A::Diagonal, p::Real=2)
+    if p == 2
+        return opnorm2(A)
+    elseif p == 1
+        return _opnorm1(A)
+    elseif p == Inf
+        return _opnormInf(A)
+    else
+        throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
+    end
+end
+
 # disambiguation methods: * of Diagonal and Adj/Trans AbsVec
 *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, D.diag, parent(x)))
 *(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) =

stdlib/LinearAlgebra/src/tridiag.jl

@@ -635,3 +635,90 @@ function SymTridiagonal{T}(M::Tridiagonal) where T
         throw(ArgumentError("Tridiagonal is not symmetric, cannot convert to SymTridiagonal"))
     end
 end
+
+###################
+#     opnorms     #
+###################
+
+# Tridiagonal
+
+function _opnorm1(A::Tridiagonal{T}) where T
+    n = size(A, 1)
+    Tnorm = typeof(float(real(zero(T))))
+    Tsum = promote_type(Float64, Tnorm)
+
+    #first col
+    nrm::Tsum = norm(A.d[1]) + norm(A.dl[1])
+    @inbounds begin
+        for i = 2:n-1
+            nrmj::Tsum = norm(A.d[i]) + norm(A.dl[i]) + norm(A.du[i-1])
+            nrm = max(nrm,nrmj)
+        end
+    end
+
+    # last col
+    @inbounds nrm = max(nrm, norm(A.d[n])+norm(A.du[n-1]))
+    return convert(Tnorm, nrm)
+end
+
+function _opnormInf(A::Tridiagonal{T}) where T
+    n = size(A, 1)
+    Tnorm = typeof(float(real(zero(T))))
+    Tsum = promote_type(Float64, Tnorm)
+
+    #first row
+    nrm::Tsum = norm(A.d[1]) + norm(A.du[1])
+    @inbounds begin
+        for i = 2:n-1
+            nrmj::Tsum = norm(A.d[i]) + norm(A.du[i]) + norm(A.dl[i-1])
+            nrm = max(nrm,nrmj)
+        end
+    end
+
+    # last row
+    @inbounds nrm = max(nrm, norm(A.d[n])+norm(A.dl[n-1]))
+    return convert(Tnorm, nrm)
+end
+
+function opnorm(A::Tridiagonal, p::Real=2)
+    if p == 2
+        return opnorm2(A)
+    elseif p == 1
+        return _opnorm1(A)
+    elseif p == Inf
+        return _opnormInf(A)
+    else
+        throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
+    end
+end
+
+# SymTridiagonal
+
+function _opnormInf1(A::SymTridiagonal{T}) where T
+    n = size(A, 1)
+    Tnorm = typeof(float(real(zero(T))))
+    Tsum = promote_type(Float64, Tnorm)
+
+    #first col/row
+    nrm::Tsum = norm(A.dv[1]) + norm(A.ev[1])
+    @inbounds begin
+        for i = 2:n-1
+            nrmj::Tsum = norm(A.dv[i]) + norm(A.ev[i-1]) + norm(A.ev[i])
+            nrm = max(nrm,nrmj)
+        end
+    end
+
+    # last col/row
+    @inbounds nrm = max(nrm, norm(A.dv[n])+norm(A.ev[n-1]))
+    return convert(Tnorm, nrm)
+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 _opnormInf1(A)
+    else
+        throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
+    end
+end

stdlib/LinearAlgebra/test/bidiag.jl

@@ -401,4 +401,19 @@ end
     @test vcat((Aub\bb)...) ≈ UpperTriangular(A)\b
 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)
+end
+
+
 end # module TestBidiagonal

stdlib/LinearAlgebra/test/diagonal.jl

@@ -490,4 +490,12 @@ 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)
+end
+
 end # module TestDiagonal

stdlib/LinearAlgebra/test/tridiag.jl

@@ -437,4 +437,18 @@ 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)
+end
+
 end # module TestTridiagonal