Skip to content

Commit d2209fe

Browse files
OlivierHntOlivierHnt
committed
Add a LinearAlgebra extension
1 parent c54a880 commit d2209fe

File tree

3 files changed

+131
-129
lines changed

3 files changed

+131
-129
lines changed

‎Project.toml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
2121
IntervalArithmeticDiffRulesExt = "DiffRules"
2222
IntervalArithmeticForwardDiffExt = "ForwardDiff"
2323
IntervalArithmeticIntervalSetsExt = "IntervalSets"
24+
IntervalArithmeticLinearAlgebraExt = "LinearAlgebra"
2425
IntervalArithmeticRecipesBaseExt = "RecipesBase"
2526
IntervalArithmeticSparseArraysExt = "SparseArrays"
2627

‎ext/IntervalArithmeticLinearAlgebra.jl

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
module IntervalArithmeticLinearAlgebraExt
2+
3+
using IntervalArithmetic
4+
import LinearAlgebra
5+
6+
# contructor for `UniformScaling`
7+
8+
IntervalArithmetic.interval(::Type{T}, J::LinearAlgebra.UniformScaling, d::IntervalArithmetic.Decoration = com; format::Symbol = :infsup) where {T} =
9+
LinearAlgebra.UniformScaling(interval(T, J.λ, d; format = format))
10+
IntervalArithmetic.interval(J::LinearAlgebra.UniformScaling, d::IntervalArithmetic.Decoration = com; format::Symbol = :infsup) =
11+
LinearAlgebra.UniformScaling(interval(J.λ, d; format = format))
12+
13+
# by-pass generic `opnorm` from LinearAlgebra to prevent NG flag
14+
15+
function LinearAlgebra.opnorm1(A::AbstractMatrix{T}) where {T<:RealOrComplexI}
16+
LinearAlgebra.require_one_based_indexing(A)
17+
m, n = size(A)
18+
Tnorm = typeof(float(real(zero(T))))
19+
Tsum = promote_type(Float64, Tnorm)
20+
nrm = zero(Tsum)
21+
@inbounds begin
22+
for j = 1:n
23+
nrmj = zero(Tsum)
24+
for i = 1:m
25+
nrmj += LinearAlgebra.norm(A[i,j])
26+
end
27+
nrm = max(nrm, nrmj)
28+
end
29+
end
30+
return convert(Tnorm, nrm)
31+
end
32+
33+
function LinearAlgebra.opnormInf(A::AbstractMatrix{T}) where {T<:RealOrComplexI}
34+
LinearAlgebra.require_one_based_indexing(A)
35+
m, n = size(A)
36+
Tnorm = typeof(float(real(zero(T))))
37+
Tsum = promote_type(Float64, Tnorm)
38+
nrm = zero(Tsum)
39+
@inbounds begin
40+
for i = 1:m
41+
nrmi = zero(Tsum)
42+
for j = 1:n
43+
nrmi += LinearAlgebra.norm(A[i,j])
44+
end
45+
nrm = max(nrm, nrmi)
46+
end
47+
end
48+
return convert(Tnorm, nrm)
49+
end
50+
51+
# matrix eigenvalues
52+
53+
function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby)
54+
# note: this function does not overwrite `A`
55+
v = _eigvals(A, permute, scale, sortby)
56+
isreal(v) && return v
57+
_fold_conjugate!(v)
58+
isreal(v) && return real(v)
59+
return v
60+
end
61+
62+
LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
63+
# note: this function does not overwrite `A`
64+
_eigvals(A, permute, scale, sortby)
65+
66+
function _eigvals(A, permute, scale, sortby)
67+
# Gershgorin circle theorem
68+
B = _similarity_transform(A, permute, scale, sortby)
69+
v = LinearAlgebra.diag(B)
70+
T = eltype(v)
71+
for j ∈ axes(B, 1)
72+
r = zero(T)
73+
for i ∈ axes(B, 2)
74+
if i ≠ j
75+
r += abs(B[i,j])
76+
end
77+
end
78+
v[j] = interval(v[j], r; format = :midpoint)
79+
end
80+
return v
81+
end
82+
83+
function _similarity_transform(A, permute, scale, sortby)
84+
mA = mid.(A)
85+
mλ, mV = LinearAlgebra.eigen(mA; permute = permute, scale = scale, sortby = sortby)
86+
mλ .+= LinearAlgebra.diag(mV \ (mA * mV - mV * LinearAlgebra.Diagonal(mλ)))
87+
Λ = LinearAlgebra.Diagonal(interval(mλ))
88+
V = interval(mV)
89+
V .= Λ .+ inv(V) * (A * V - V * Λ)
90+
return V
91+
end
92+
93+
function _fold_conjugate!(v)
94+
for i ∈ eachindex(v)
95+
vᵢ = v[i]
96+
idxs = findall(j -> (j ≠ i) & !isdisjoint_interval(conj(vᵢ), v[j]), eachindex(v))
97+
if isempty(idxs)
98+
v[i] = real(vᵢ)
99+
else
100+
w = view(v, idxs)
101+
z = conj(intersect_interval(conj(vᵢ), reduce(intersect_interval, w)))
102+
z = complex(IntervalArithmetic.setdecoration(real(z), min(decoration(real(vᵢ)), minimum(decoration ∘ real, w))), IntervalArithmetic.setdecoration(imag(z), min(decoration(imag(vᵢ)), minimum(decoration ∘ imag, w))))
103+
v[i] = z
104+
end
105+
end
106+
return v
107+
end
108+
109+
# matrix determinant
110+
111+
LinearAlgebra.det(A::AbstractMatrix{<:Interval}) = real(reduce(*, LinearAlgebra.eigvals(A)))
112+
LinearAlgebra.det(A::AbstractMatrix{<:Complex{<:Interval}}) = reduce(*, LinearAlgebra.eigvals(A))
113+
114+
end

‎src/matmul.jl

Lines changed: 16 additions & 129 deletions
Original file line numberDiff line numberDiff line change
@@ -1,50 +1,3 @@
1-
interval(::Type{T}, J::LinearAlgebra.UniformScaling, d::Decoration = com; format::Symbol = :infsup) where {T} =
2-
LinearAlgebra.UniformScaling(interval(T, J.λ, d; format = format))
3-
interval(J::LinearAlgebra.UniformScaling, d::Decoration = com; format::Symbol = :infsup) =
4-
LinearAlgebra.UniformScaling(interval(J.λ, d; format = format))
5-
6-
7-
8-
# by-pass generic `opnorm` from LinearAlgebra to prevent NG flag
9-
10-
function LinearAlgebra.opnorm1(A::AbstractMatrix{T}) where {T<:RealOrComplexI}
11-
LinearAlgebra.require_one_based_indexing(A)
12-
m, n = size(A)
13-
Tnorm = typeof(float(real(zero(T))))
14-
Tsum = promote_type(Float64, Tnorm)
15-
nrm = zero(Tsum)
16-
@inbounds begin
17-
for j = 1:n
18-
nrmj = zero(Tsum)
19-
for i = 1:m
20-
nrmj += LinearAlgebra.norm(A[i,j])
21-
end
22-
nrm = max(nrm, nrmj)
23-
end
24-
end
25-
return convert(Tnorm, nrm)
26-
end
27-
28-
function LinearAlgebra.opnormInf(A::AbstractMatrix{T}) where {T<:RealOrComplexI}
29-
LinearAlgebra.require_one_based_indexing(A)
30-
m, n = size(A)
31-
Tnorm = typeof(float(real(zero(T))))
32-
Tsum = promote_type(Float64, Tnorm)
33-
nrm = zero(Tsum)
34-
@inbounds begin
35-
for i = 1:m
36-
nrmi = zero(Tsum)
37-
for j = 1:n
38-
nrmi += LinearAlgebra.norm(A[i,j])
39-
end
40-
nrm = max(nrm, nrmi)
41-
end
42-
end
43-
return convert(Tnorm, nrm)
44-
end
45-
46-
47-
481
# matrix inversion
492
# note: use the contraction mapping theorem, only works when the entries of A have small radii
503

@@ -63,88 +16,7 @@ function Base.inv(A::Matrix{<:RealOrComplexI})
6316
return A⁻¹
6417
end
6518

66-
67-
68-
# matrix eigenvalues
69-
70-
function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby)
71-
# note: this function does not overwrite `A`
72-
v = _eigvals(A, permute, scale, sortby)
73-
isreal(v) && return v
74-
_fold_conjugate!(v)
75-
isreal(v) && return real(v)
76-
return v
77-
end
78-
79-
LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
80-
# note: this function does not overwrite `A`
81-
_eigvals(A, permute, scale, sortby)
82-
83-
function _eigvals(A, permute, scale, sortby)
84-
# Gershgorin circle theorem
85-
B = _similarity_transform(A, permute, scale, sortby)
86-
v = LinearAlgebra.diag(B)
87-
T = eltype(v)
88-
for j ∈ axes(B, 1)
89-
r = zero(T)
90-
for i ∈ axes(B, 2)
91-
if i ≠ j
92-
r += abs(B[i,j])
93-
end
94-
end
95-
v[j] = interval(v[j], r; format = :midpoint)
96-
end
97-
return v
98-
end
99-
100-
function _similarity_transform(A, permute, scale, sortby)
101-
mA = mid.(A)
102-
mλ, mV = LinearAlgebra.eigen(mA; permute = permute, scale = scale, sortby = sortby)
103-
mλ .+= LinearAlgebra.diag(mV \ (mA * mV - mV * LinearAlgebra.Diagonal(mλ)))
104-
Λ = LinearAlgebra.Diagonal(interval(mλ))
105-
V = interval(mV)
106-
V .= Λ .+ inv(V) * (A * V - V * Λ)
107-
return V
108-
end
109-
110-
function _fold_conjugate!(v)
111-
for i ∈ eachindex(v)
112-
vᵢ = v[i]
113-
idxs = findall(j -> (j ≠ i) & !isdisjoint_interval(conj(vᵢ), v[j]), eachindex(v))
114-
if isempty(idxs)
115-
v[i] = real(vᵢ)
116-
else
117-
w = view(v, idxs)
118-
z = conj(intersect_interval(conj(vᵢ), reduce(intersect_interval, w)))
119-
z = complex(setdecoration(real(z), min(decoration(real(vᵢ)), minimum(decoration ∘ real, w))), setdecoration(imag(z), min(decoration(imag(vᵢ)), minimum(decoration ∘ imag, w))))
120-
v[i] = z
121-
end
122-
end
123-
return v
124-
end
125-
126-
127-
128-
# matrix determinant
129-
130-
LinearAlgebra.det(A::AbstractMatrix{<:Interval}) = real(reduce(*, LinearAlgebra.eigvals(A)))
131-
LinearAlgebra.det(A::AbstractMatrix{<:Complex{<:Interval}}) = reduce(*, LinearAlgebra.eigvals(A))
132-
133-
134-
135-
# matrix multiplication
136-
137-
"""
138-
MatMulMode
139-
140-
Matrix multiplication type.
141-
142-
Available mode types:
143-
- `:slow` (default): generic algorithm.
144-
- `:fast` : Rump's algorithm.
145-
"""
146-
struct MatMulMode{T} end
147-
19+
#
14820
# by-pass `similar` methods defined in array.jl
14921
# note: written in this form to avoid by-passing the default behaviour for `Union{}`
15022
Base.similar(a::Array{Interval{T},1}) where {T<:NumTypes} = zeros(Interval{T}, size(a, 1))
@@ -169,6 +41,21 @@ Base.similar(::Array, S::Type{Interval{T}}, dims::Dims) where {T<:NumTy
16941
Base.similar(::Array, S::Type{Complex{Interval{T}}}, dims::Dims) where {T<:NumTypes} = zeros(S, dims)
17042
#
17143

44+
# matrix multiplication
45+
46+
"""
47+
MatMulMode
48+
49+
Matrix multiplication mode type.
50+
51+
Available mode types:
52+
- `:slow` (default): generic algorithm.
53+
- `:fast` : Rump's algorithm.
54+
"""
55+
struct MatMulMode{T} end
56+
57+
#
58+
17259
LinearAlgebra.mul!(C::AbstractVecOrMat{<:RealOrComplexI}, A::AbstractMatrix{<:RealOrComplexI}, B::AbstractVecOrMat{<:RealOrComplexI}) =
17360
LinearAlgebra.mul!(C, A, B, interval(true), interval(false))
17461

0 commit comments

Comments
 (0)