diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl index 5b7264558f9ae..2fb1415fff590 100644 --- a/stdlib/LinearAlgebra/src/bidiag.jl +++ b/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 # #################### diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 243df4d82eec2..ec1cd89c05781 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -1038,6 +1038,18 @@ end *(A::AbstractTriangular, D::Diagonal) = @invoke *(A::AbstractMatrix, D::Diagonal) *(D::Diagonal, A::AbstractTriangular) = @invoke *(D::Diagonal, A::AbstractMatrix) +_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) + +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) diff --git a/stdlib/LinearAlgebra/src/special.jl b/stdlib/LinearAlgebra/src/special.jl index c61586a810140..636fd8d1a9f43 100644 --- a/stdlib/LinearAlgebra/src/special.jl +++ b/stdlib/LinearAlgebra/src/special.jl @@ -593,3 +593,15 @@ end _istril(A::LowerTriangular{<:Any, <:BandedMatrix}, k) = istril(parent(A), k) _istriu(A::UpperTriangular{<:Any, <:BandedMatrix}, k) = istriu(parent(A), k) _istriu(A::UpperHessenberg{<:Any, <:BandedMatrix}, k) = istriu(parent(A), k) +# 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 diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl index 0d73e6dd46fdb..a24cc50b42b9f 100644 --- a/stdlib/LinearAlgebra/src/tridiag.jl +++ b/stdlib/LinearAlgebra/src/tridiag.jl @@ -1097,3 +1097,38 @@ function show(io::IO, S::SymTridiagonal) show(io, S.ev) print(io, ")") 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 diff --git a/stdlib/LinearAlgebra/test/bidiag.jl b/stdlib/LinearAlgebra/test/bidiag.jl index df30748e042b5..d87af970e8f2e 100644 --- a/stdlib/LinearAlgebra/test/bidiag.jl +++ b/stdlib/LinearAlgebra/test/bidiag.jl @@ -1138,4 +1138,31 @@ 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 diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index 16f3d2287f317..c29728b518b32 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -1452,4 +1452,24 @@ end @test kron(D2, D) == kron(Array{eltype(D2)}(D2), D) 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 + end # module TestDiagonal diff --git a/stdlib/LinearAlgebra/test/tridiag.jl b/stdlib/LinearAlgebra/test/tridiag.jl index dc14ddb1d1b27..4b592a87a6b19 100644 --- a/stdlib/LinearAlgebra/test/tridiag.jl +++ b/stdlib/LinearAlgebra/test/tridiag.jl @@ -1075,4 +1075,28 @@ end @test all(iszero, diag(A,1)) 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 + end # module TestTridiagonal