Skip to content

Commit c497aaa

Browse files
mcognettamcognetta
committed
changing the opnorms for structured matrices to use mapreduce
1 parent dc2313e commit c497aaa

File tree

6 files changed

+57
-123
lines changed

6 files changed

+57
-123
lines changed

‎stdlib/LinearAlgebra/src/bidiag.jl‎

Lines changed: 2 additions & 60 deletions
Original file line numberDiff line numberDiff line change
@@ -203,69 +203,11 @@ function svd(M::Bidiagonal; full::Bool = false)
203203
end
204204

205205
###################
206-
# opnorms #
206+
# opnorm #
207207
###################
208208

209-
# helpers for upper and lower bidiagonal matrices
210-
211-
# function _upopnorm1(A::Bidiagonal{T}) where T
212-
# n = size(A, 1)
213-
# Tnorm = typeof(float(real(zero(T))))
214-
# Tsum = promote_type(Float64, Tnorm)
215-
# nrm::Tsum = A.dv[1]
216-
# @inbounds begin
217-
# for i = 2:n
218-
# nrmj::Tsum = norm(A.dv[i]) + norm(A.ev[i-1])
219-
# nrm = max(nrm,nrmj)
220-
# end
221-
# end
222-
# return convert(Tnorm, nrm)
223-
# end
224-
225-
# function _loopnorm1(A::Bidiagonal{T}) where T
226-
# n = size(A, 1)
227-
# Tnorm = typeof(float(real(zero(T))))
228-
# Tsum = promote_type(Float64, Tnorm)
229-
# nrm::Tsum = A.dv[end]
230-
# @inbounds begin
231-
# for i = 1:n-1
232-
# nrmj::Tsum = norm(A.dv[i]) + norm(A.ev[i])
233-
# nrm = max(nrm,nrmj)
234-
# end
235-
# end
236-
# return convert(Tnorm, nrm)
237-
# end
238-
239-
# function _upopnormInf(A::Bidiagonal{T}) where T
240-
# n = size(A, 1)
241-
# Tnorm = typeof(float(real(zero(T))))
242-
# Tsum = promote_type(Float64, Tnorm)
243-
# nrm::Tsum = A.dv[end]
244-
# @inbounds begin
245-
# for i = 1:n-1
246-
# nrmj::Tsum = norm(A.dv[i]) + norm(A.ev[i])
247-
# nrm = max(nrm,nrmj)
248-
# end
249-
# end
250-
# return convert(Tnorm, nrm)
251-
# end
252-
253-
# function _loopnormInf(A::Bidiagonal{T}) where T
254-
# n = size(A, 1)
255-
# Tnorm = typeof(float(real(zero(T))))
256-
# Tsum = promote_type(Float64, Tnorm)
257-
# nrm::Tsum = A.dv[1]
258-
# @inbounds begin
259-
# for i = 2:n
260-
# nrmj::Tsum = norm(A.dv[i]) + norm(A.ev[i-1])
261-
# nrm = max(nrm,nrmj)
262-
# end
263-
# end
264-
# return convert(Tnorm, nrm)
265-
# end
266-
267209
function _opnorm1Inf(B::Bidiagonal, p)
268-
size(B, 1) == 1 || return max(B)
210+
size(B, 1) == 1 && return norm(first(B.dv))
269211
case = xor(p == 1, istriu(B))
270212
normd1, normdend = norm(first(B.dv)), norm(last(B.dv))
271213
normd1, normdend = case ? (zero(normd1), normdend) : (normd1, zero(normdend))

‎stdlib/LinearAlgebra/src/diagonal.jl‎

Lines changed: 3 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -548,15 +548,12 @@ function svd(D::Diagonal{<:Number})
548548
return SVD(Up, S[piv], copy(Vp'))
549549
end
550550

551-
_opnorm1(A::Diagonal) = maximum(norm(x) for x in A.diag)
552-
_opnormInf(A::Diagonal) = maximum(norm(x) for x in A.diag)
551+
_opnorm1Inf(A::Diagonal) = maximum(norm(x) for x in A.diag)
553552
function opnorm(A::Diagonal, p::Real=2)
554553
if p == 2
555554
return opnorm2(A)
556-
elseif p == 1
557-
return _opnorm1(A)
558-
elseif p == Inf
559-
return _opnormInf(A)
555+
elseif p == 1 || p == Inf
556+
return _opnorm1Inf(A)
560557
else
561558
throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
562559
end

‎stdlib/LinearAlgebra/src/tridiag.jl‎

Lines changed: 25 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -642,82 +642,50 @@ end
642642

643643
# Tridiagonal
644644

645-
function _opnorm1(A::Tridiagonal{T}) where T
646-
n = size(A, 1)
647-
Tnorm = typeof(float(real(zero(T))))
648-
Tsum = promote_type(Float64, Tnorm)
649-
650-
#first col
651-
nrm::Tsum = norm(A.d[1]) + norm(A.dl[1])
652-
@inbounds begin
653-
for i = 2:n-1
654-
nrmj::Tsum = norm(A.d[i]) + norm(A.dl[i]) + norm(A.du[i-1])
655-
nrm = max(nrm,nrmj)
656-
end
657-
end
658-
659-
# last col
660-
@inbounds nrm = max(nrm, norm(A.d[n])+norm(A.du[n-1]))
661-
return convert(Tnorm, nrm)
662-
end
663-
664-
function _opnormInf(A::Tridiagonal{T}) where T
665-
n = size(A, 1)
666-
Tnorm = typeof(float(real(zero(T))))
667-
Tsum = promote_type(Float64, Tnorm)
668-
669-
#first row
670-
nrm::Tsum = norm(A.d[1]) + norm(A.du[1])
671-
@inbounds begin
672-
for i = 2:n-1
673-
nrmj::Tsum = norm(A.d[i]) + norm(A.du[i]) + norm(A.dl[i-1])
674-
nrm = max(nrm,nrmj)
675-
end
676-
end
677-
678-
# last row
679-
@inbounds nrm = max(nrm, norm(A.d[n])+norm(A.dl[n-1]))
680-
return convert(Tnorm, nrm)
645+
function _opnorm1Inf(A::Tridiagonal, p)
646+
size(A, 1) == 1 && return norm(first(A.d))
647+
case = p == Inf
648+
lowerrange, upperrange = case ? (1:length(A.dl)-1, 2:length(A.dl)) : (2:length(A.dl), 1:length(A.dl)-1)
649+
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)))
650+
651+
return max(
652+
mapreduce(t -> sum(norm, t),
653+
max,
654+
zip(view(A.d, (2:length(A.d)-1)), view(A.dl, lowerrange), view(A.du, upperrange))
655+
),
656+
normfirst, normend)
681657
end
682658

683659
function opnorm(A::Tridiagonal, p::Real=2)
684660
if p == 2
685661
return opnorm2(A)
686-
elseif p == 1
687-
return _opnorm1(A)
688-
elseif p == Inf
689-
return _opnormInf(A)
662+
elseif p == 1 || p == Inf
663+
_opnorm1Inf(A, p)
690664
else
691665
throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
692666
end
693667
end
694668

695669
# SymTridiagonal
696670

697-
function _opnormInf1(A::SymTridiagonal{T}) where T
698-
n = size(A, 1)
699-
Tnorm = typeof(float(real(zero(T))))
700-
Tsum = promote_type(Float64, Tnorm)
701-
702-
#first col/row
703-
nrm::Tsum = norm(A.dv[1]) + norm(A.ev[1])
704-
@inbounds begin
705-
for i = 2:n-1
706-
nrmj::Tsum = norm(A.dv[i]) + norm(A.ev[i-1]) + norm(A.ev[i])
707-
nrm = max(nrm,nrmj)
708-
end
709-
end
671+
function _opnorm1Inf(A::SymTridiagonal)
672+
size(A, 1) == 1 && return norm(first(A.dv))
673+
lowerrange, upperrange = 1:length(A.ev)-1, 2:length(A.ev)
674+
normfirst, normend = norm(first(A.dv))+norm(first(A.ev)), norm(last(A.ev))+norm(last(A.dv))
710675

711-
# last col/row
712-
@inbounds nrm = max(nrm, norm(A.dv[n])+norm(A.ev[n-1]))
713-
return convert(Tnorm, nrm)
676+
return max(
677+
mapreduce(t -> sum(norm, t),
678+
max,
679+
zip(view(A.dv, (2:length(A.dv)-1)), view(A.ev, lowerrange), view(A.ev, upperrange))
680+
),
681+
normfirst, normend)
714682
end
715683

716684
function opnorm(A::SymTridiagonal, p::Real=2)
717685
if p == 2
718686
return opnorm2(A)
719687
elseif p == 1 || p == Inf # these are the same for symmetric matrices
720-
return _opnormInf1(A)
688+
return _opnorm1Inf(A)
721689
else
722690
throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
723691
end

‎stdlib/LinearAlgebra/test/bidiag.jl‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -413,6 +413,18 @@ end
413413
@test opnorm(B, 1) == opnorm(Matrix(B), 1)
414414
@test opnorm(B, 2) ≈ opnorm(Matrix(B), 2)
415415
@test opnorm(B, Inf) == opnorm(Matrix(B), Inf)
416+
417+
B = Bidiagonal([2], Int[], 'L')
418+
419+
@test opnorm(B, 1) == opnorm(Matrix(B), 1)
420+
@test opnorm(B, 2) ≈ opnorm(Matrix(B), 2)
421+
@test opnorm(B, Inf) == opnorm(Matrix(B), Inf)
422+
423+
B = Bidiagonal([2], Int[], 'U')
424+
425+
@test opnorm(B, 1) == opnorm(Matrix(B), 1)
426+
@test opnorm(B, 2) ≈ opnorm(Matrix(B), 2)
427+
@test opnorm(B, Inf) == opnorm(Matrix(B), Inf)
416428
end
417429

418430

‎stdlib/LinearAlgebra/test/diagonal.jl‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -504,6 +504,11 @@ end
504504
@test opnorm(D, 1) == opnorm(Matrix(D), 1)
505505
@test opnorm(D, 2) ≈ opnorm(Matrix(D), 2)
506506
@test opnorm(D, Inf) == opnorm(Matrix(D), Inf)
507+
508+
D = Diagonal([-11])
509+
@test opnorm(D, 1) == opnorm(Matrix(D), 1)
510+
@test opnorm(D, 2) ≈ opnorm(Matrix(D), 2)
511+
@test opnorm(D, Inf) == opnorm(Matrix(D), Inf)
507512
end
508513

509514
end # module TestDiagonal

‎stdlib/LinearAlgebra/test/tridiag.jl‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -449,6 +449,16 @@ end
449449
@test opnorm(S, 1) == opnorm(Matrix(S), 1)
450450
@test_skip opnorm(S, 2) ≈ opnorm(Matrix(S), 2) # currently missing
451451
@test opnorm(S, Inf) == opnorm(Matrix(S), Inf)
452+
453+
T = Tridiagonal(Int[], [-5], Int[])
454+
@test opnorm(T, 1) == opnorm(Matrix(T), 1)
455+
@test_skip opnorm(T, 2) ≈ opnorm(Matrix(T), 2) # currently missing
456+
@test opnorm(T, Inf) == opnorm(Matrix(T), Inf)
457+
458+
S = SymTridiagonal(T)
459+
@test opnorm(S, 1) == opnorm(Matrix(S), 1)
460+
@test_skip opnorm(S, 2) ≈ opnorm(Matrix(S), 2) # currently missing
461+
@test opnorm(S, Inf) == opnorm(Matrix(S), Inf)
452462
end
453463

454464
end # module TestTridiagonal

0 commit comments

Comments
 (0)