From d9dcd1ed8d7160c379dd7aae914b5f5a678dff73 Mon Sep 17 00:00:00 2001 From: OlivierHnt <38465572+OlivierHnt@users.noreply.github.com> Date: Tue, 22 Oct 2024 16:53:55 +0900 Subject: [PATCH] Make LinearAlgebra a dependency and add `MatMulMode` --- Project.toml | 3 +- src/IntervalArithmetic.jl | 28 ++- src/intervals/intervals.jl | 2 +- .../matmul.jl | 233 ++++++++++-------- test/interval_tests/linearalgebra.jl | 2 + 5 files changed, 162 insertions(+), 106 deletions(-) rename ext/IntervalArithmeticLinearAlgebraExt.jl => src/matmul.jl (63%) diff --git a/Project.toml b/Project.toml index 3d0d0dcb9..1fe0ce1ea 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.22.17" [deps] CRlibm_jll = "4e9b3aee-d8a1-5a3d-ad8b-7d824db253f0" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" RoundingEmulator = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705" @@ -13,14 +14,12 @@ DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [extensions] IntervalArithmeticDiffRulesExt = "DiffRules" IntervalArithmeticForwardDiffExt = "ForwardDiff" IntervalArithmeticRecipesBaseExt = "RecipesBase" IntervalArithmeticIntervalSetsExt = "IntervalSets" -IntervalArithmeticLinearAlgebraExt = "LinearAlgebra" [compat] CRlibm_jll = "1" diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl index 90ee99610..11abb953d 100644 --- a/src/IntervalArithmetic.jl +++ b/src/IntervalArithmetic.jl @@ -1,3 +1,14 @@ +""" + IntervalArithmetic + +Library for validated numerics using interval arithmetic. + +Default settings: +- power mode is set to `:fast` +- matrix multiplication mode is set to `:slow` + +Learn more: https://github.com/JuliaIntervals/IntervalArithmetic.jl +""" module IntervalArithmetic import CRlibm_jll @@ -8,16 +19,23 @@ using MacroTools: MacroTools, prewalk, postwalk, @capture # include("intervals/intervals.jl") +# convenient alias +const RealOrComplexI{T} = Union{Interval{T},Complex{Interval{T}}} +const ComplexI{T} = Complex{Interval{T}} +const RealIntervalType{T} = Union{BareInterval{T},Interval{T}} + export RealOrComplexI, ComplexI, RealIntervalType + +# include("display.jl") export setdisplay +# + include("symbols.jl") -# convenient alias -const RealOrComplexI{T} = Union{Interval{T},Complex{Interval{T}}} -const ComplexI{T} = Complex{Interval{T}} -const RealIntervalType{T} = Union{BareInterval{T},Interval{T}} - export RealOrComplexI, ComplexI, RealIntervalType +# + +include("matmul.jl") end diff --git a/src/intervals/intervals.jl b/src/intervals/intervals.jl index 4daf7142d..cb4685e92 100644 --- a/src/intervals/intervals.jl +++ b/src/intervals/intervals.jl @@ -1,7 +1,7 @@ # Construction and composability with numbers include("construction.jl") export BareInterval, bareinterval, decoration, ill, trv, def, dac, com, - Interval, interval, isguaranteed, @interval, @I_str + Interval, interval, isguaranteed, @interval, @I_str include("parsing.jl") include("real_interface.jl") include("exact_literals.jl") diff --git a/ext/IntervalArithmeticLinearAlgebraExt.jl b/src/matmul.jl similarity index 63% rename from ext/IntervalArithmeticLinearAlgebraExt.jl rename to src/matmul.jl index 15d9d598d..3ad84cc10 100644 --- a/ext/IntervalArithmeticLinearAlgebraExt.jl +++ b/src/matmul.jl @@ -1,11 +1,15 @@ -module IntervalArithmeticLinearAlgebraExt +import LinearAlgebra + + + +# + +interval(::Type{T}, J::LinearAlgebra.UniformScaling, d::Decoration = com; format::Symbol = :infsup) where {T} = + LinearAlgebra.UniformScaling(interval(T, J.λ, d; format = format)) +interval(J::LinearAlgebra.UniformScaling, d::Decoration = com; format::Symbol = :infsup) = + LinearAlgebra.UniformScaling(interval(J.λ, d; format = format)) -using IntervalArithmetic, LinearAlgebra -IntervalArithmetic.interval(::Type{T}, J::UniformScaling, d::IntervalArithmetic.Decoration = com; format::Symbol = :infsup) where {T} = - UniformScaling(interval(T, J.λ, d; format = format)) -IntervalArithmetic.interval(J::UniformScaling, d::IntervalArithmetic.Decoration = com; format::Symbol = :infsup) = - UniformScaling(interval(J.λ, d; format = format)) # by-pass generic `opnorm` from LinearAlgebra to prevent NG flag @@ -19,7 +23,7 @@ function LinearAlgebra.opnorm1(A::AbstractMatrix{T}) where {T<:RealOrComplexI} for j = 1:n nrmj = zero(Tsum) for i = 1:m - nrmj += norm(A[i,j]) + nrmj += LinearAlgebra.norm(A[i,j]) end nrm = max(nrm, nrmj) end @@ -37,7 +41,7 @@ function LinearAlgebra.opnormInf(A::AbstractMatrix{T}) where {T<:RealOrComplexI} for i = 1:m nrmi = zero(Tsum) for j = 1:n - nrmi += norm(A[i,j]) + nrmi += LinearAlgebra.norm(A[i,j]) end nrm = max(nrm, nrmi) end @@ -45,20 +49,7 @@ function LinearAlgebra.opnormInf(A::AbstractMatrix{T}) where {T<:RealOrComplexI} return convert(Tnorm, nrm) end -# convenient function to propagate NG flag - -function _ensure_ng_flag!(C::AbstractVecOrMat{<:Interval}, ng_flag::Bool) - C .= IntervalArithmetic._unsafe_interval.(getfield.(C, :bareinterval), decoration.(C), ng_flag) - return C -end -function _ensure_ng_flag!(C::AbstractVecOrMat{<:Complex{<:Interval}}, ng_flag::Bool) - C .= complex.( - IntervalArithmetic._unsafe_interval.(getfield.(real.(C), :bareinterval), decoration.(C), ng_flag), - IntervalArithmetic._unsafe_interval.(getfield.(imag.(C), :bareinterval), decoration.(C), ng_flag) - ) - return C -end # matrix inversion # Note: use the contraction mapping theorem, only works when the entries of A have small radii @@ -66,9 +57,9 @@ end function Base.inv(A::Matrix{<:RealOrComplexI}) mid_A = mid.(A) approx_A⁻¹ = interval(inv(mid_A)) - F = A * approx_A⁻¹ - I - Y = opnorm(approx_A⁻¹ * F, Inf) - Z₁ = opnorm(F, Inf) + F = A * approx_A⁻¹ - interval(LinearAlgebra.I) + Y = LinearAlgebra.opnorm(approx_A⁻¹ * F, Inf) + Z₁ = LinearAlgebra.opnorm(F, Inf) if isbounded(Y) & strictprecedes(Z₁, one(one(Z₁))) A⁻¹ = interval.(approx_A⁻¹, inf(interval(mag(Y)) / (one(Z₁) - interval(mag(Z₁)))); format = :midpoint) else @@ -78,53 +69,88 @@ function Base.inv(A::Matrix{<:RealOrComplexI}) return A⁻¹ end -# fast matrix multiplication -# Note: Rump's algorithm + + +# matrix multiplication + +""" + MatMulMode + +Matrix multiplication type. + +Available mode types: +- `:slow` (default): generic algorithm. +- `:fast` : Rump's algorithm. +""" +struct MatMulMode{T} end + +matmul_mode() = MatMulMode{:slow}() + +# function Base.:*(A::AbstractMatrix{<:RealOrComplexI}, B::AbstractMatrix{<:RealOrComplexI}) T = promote_type(eltype(A), eltype(B)) C = zeros(T, size(A, 1), size(B, 2)) - return mul!(C, A, B, one(T), zero(T)) + return LinearAlgebra.mul!(C, A, B, one(T), zero(T)) end function Base.:*(A::AbstractMatrix{<:RealOrComplexI}, B::AbstractVector{<:RealOrComplexI}) T = promote_type(eltype(A), eltype(B)) C = zeros(T, size(A, 1)) - return mul!(C, A, B, one(T), zero(T)) + return LinearAlgebra.mul!(C, A, B, one(T), zero(T)) end function Base.:*(A::AbstractMatrix{<:RealOrComplexI}, B::AbstractMatrix) T = promote_type(eltype(A), eltype(B)) C = zeros(T, size(A, 1), size(B, 2)) - return mul!(C, A, B, one(T), zero(T)) + return LinearAlgebra.mul!(C, A, B, one(T), zero(T)) end function Base.:*(A::AbstractMatrix{<:RealOrComplexI}, B::AbstractVector) T = promote_type(eltype(A), eltype(B)) C = zeros(T, size(A, 1)) - return mul!(C, A, B, one(T), zero(T)) + return LinearAlgebra.mul!(C, A, B, one(T), zero(T)) end function Base.:*(A::AbstractMatrix, B::AbstractMatrix{<:RealOrComplexI}) T = promote_type(eltype(A), eltype(B)) C = zeros(T, size(A, 1), size(B, 2)) - return mul!(C, A, B, one(T), zero(T)) + return LinearAlgebra.mul!(C, A, B, one(T), zero(T)) end function Base.:*(A::AbstractMatrix, B::AbstractVector{<:RealOrComplexI}) T = promote_type(eltype(A), eltype(B)) C = zeros(T, size(A, 1)) - return mul!(C, A, B, one(T), zero(T)) + return LinearAlgebra.mul!(C, A, B, one(T), zero(T)) end LinearAlgebra.mul!(C::AbstractVecOrMat{<:RealOrComplexI}, A::AbstractMatrix{<:RealOrComplexI}, B::AbstractVecOrMat{<:RealOrComplexI}, α::Number, β::Number) = - _mul!(C, A, B, α, β) + _mul!(matmul_mode(), C, A, B, α, β) LinearAlgebra.mul!(C::AbstractVecOrMat{<:RealOrComplexI}, A::AbstractMatrix, B::AbstractVecOrMat{<:RealOrComplexI}, α::Number, β::Number) = - _mul!(C, A, B, α, β) + _mul!(matmul_mode(), C, A, B, α, β) LinearAlgebra.mul!(C::AbstractVecOrMat{<:RealOrComplexI}, A::AbstractMatrix{<:RealOrComplexI}, B::AbstractVecOrMat, α::Number, β::Number) = - _mul!(C, A, B, α, β) + _mul!(matmul_mode(), C, A, B, α, β) + +_mul!(::MatMulMode{:slow}, C, A, B, α, β) = LinearAlgebra._mul!(C, A, B, α, β) + +# fast matrix multiplication +# Note: Rump's algorithm + +_mul!(::MatMulMode{:fast}, C, A::AbstractMatrix{<:Interval{<:Rational}}, B::AbstractVecOrMat{<:Interval{<:Rational}}, α, β) = + LinearAlgebra._mul!(C, A, B, α, β) +_mul!(::MatMulMode{:fast}, C, A::AbstractMatrix{<:Interval{<:Rational}}, B::AbstractVecOrMat, α, β) = + LinearAlgebra._mul!(C, A, B, α, β) +_mul!(::MatMulMode{:fast}, C, A::AbstractMatrix, B::AbstractVecOrMat{<:Interval{<:Rational}}, α, β) = + LinearAlgebra._mul!(C, A, B, α, β) + +_mul!(::MatMulMode{:fast}, C, A::AbstractMatrix{<:Complex{<:Interval{<:Rational}}}, B::AbstractVecOrMat{<:Complex{<:Interval{<:Rational}}}, α, β) = + LinearAlgebra._mul!(C, A, B, α, β) +_mul!(::MatMulMode{:fast}, C, A::AbstractMatrix{<:Complex{<:Interval{<:Rational}}}, B::AbstractVecOrMat{<:Interval{<:Rational}}, α, β) = + LinearAlgebra._mul!(C, A, B, α, β) +_mul!(::MatMulMode{:fast}, C, A::AbstractMatrix{<:Interval{<:Rational}}, B::AbstractVecOrMat{<:Complex{<:Interval{<:Rational}}}, α, β) = + LinearAlgebra._mul!(C, A, B, α, β) for (T, S) ∈ ((:Interval, :Interval), (:Interval, :Any), (:Any, :Interval)) - @eval function _mul!(C, A::AbstractMatrix{<:$T}, B::AbstractVecOrMat{<:$S}, α, β) + @eval function _mul!(::MatMulMode{:fast}, C, A::AbstractMatrix{<:$T}, B::AbstractVecOrMat{<:$S}, α, β) CoefType = eltype(C) if iszero(α) if iszero(β) @@ -133,7 +159,7 @@ for (T, S) ∈ ((:Interval, :Interval), (:Interval, :Any), (:Any, :Interval)) C .*= β end else - BoundType = IntervalArithmetic.numtype(CoefType) + BoundType = numtype(CoefType) ABinf, ABsup = __mul(A, B) if isone(α) if iszero(β) @@ -161,7 +187,7 @@ end for (T, S) ∈ ((:(Complex{<:Interval}), :(Complex{<:Interval})), (:(Complex{<:Interval}), :Complex), (:Complex, :(Complex{<:Interval}))) - @eval function _mul!(C, A::AbstractMatrix{<:$T}, B::AbstractVecOrMat{<:$S}, α, β) + @eval function _mul!(::MatMulMode{:fast}, C, A::AbstractMatrix{<:$T}, B::AbstractVecOrMat{<:$S}, α, β) CoefType = eltype(C) if iszero(α) if iszero(β) @@ -170,7 +196,7 @@ for (T, S) ∈ ((:(Complex{<:Interval}), :(Complex{<:Interval})), C .*= β end else - BoundType = IntervalArithmetic.numtype(CoefType) + BoundType = numtype(CoefType) A_real, A_imag = reim(A) B_real, B_imag = reim(B) ABinf_1, ABsup_1 = __mul(A_real, B_real) @@ -209,7 +235,7 @@ end for (T, S) ∈ ((:(Complex{<:Interval}), :Interval), (:(Complex{<:Interval}), :Any), (:Complex, :Interval)) @eval begin - function _mul!(C, A::AbstractMatrix{<:$T}, B::AbstractVecOrMat{<:$S}, α, β) + function _mul!(::MatMulMode{:fast}, C, A::AbstractMatrix{<:$T}, B::AbstractVecOrMat{<:$S}, α, β) CoefType = eltype(C) if iszero(α) if iszero(β) @@ -218,7 +244,7 @@ for (T, S) ∈ ((:(Complex{<:Interval}), :Interval), (:(Complex{<:Interval}), :A C .*= β end else - BoundType = IntervalArithmetic.numtype(CoefType) + BoundType = numtype(CoefType) A_real, A_imag = reim(A) ABinf_real, ABsup_real = __mul(A_real, B) ABinf_imag, ABsup_imag = __mul(A_imag, B) @@ -245,7 +271,7 @@ for (T, S) ∈ ((:(Complex{<:Interval}), :Interval), (:(Complex{<:Interval}), :A return C end - function _mul!(C, A::AbstractMatrix{<:$S}, B::AbstractVecOrMat{<:$T}, α, β) + function _mul!(::MatMulMode{:fast}, C, A::AbstractMatrix{<:$S}, B::AbstractVecOrMat{<:$T}, α, β) CoefType = eltype(C) if iszero(α) if iszero(β) @@ -254,7 +280,7 @@ for (T, S) ∈ ((:(Complex{<:Interval}), :Interval), (:(Complex{<:Interval}), :A C .*= β end else - BoundType = IntervalArithmetic.numtype(CoefType) + BoundType = numtype(CoefType) B_real, B_imag = reim(B) ABinf_real, ABsup_real = __mul(A, B_real) ABinf_imag, ABsup_imag = __mul(A, B_imag) @@ -283,18 +309,16 @@ for (T, S) ∈ ((:(Complex{<:Interval}), :Interval), (:(Complex{<:Interval}), :A end end -__mul(A::AbstractMatrix{<:Interval{<:Rational}}, B::AbstractVecOrMat{<:Interval{<:Rational}}) = A * B - function __mul(A::AbstractMatrix{T}, B::AbstractVecOrMat{S}) where {T,S} - NewType = IntervalArithmetic.promote_numtype(T, S) + NewType = promote_numtype(T, S) return __mul(interval.(NewType, A), interval.(NewType, B)) end function __mul(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{Interval{T}}) where {T<:AbstractFloat} - mA = IntervalArithmetic._div_round.(IntervalArithmetic._add_round.(inf.(A), sup.(A), RoundUp), convert(T, 2), RoundUp) # (inf.(A) .+ sup.(A)) ./ 2 - rA = IntervalArithmetic._sub_round.(mA, inf.(A), RoundUp) - mB = IntervalArithmetic._div_round.(IntervalArithmetic._add_round.(inf.(B), sup.(B), RoundUp), convert(T, 2), RoundUp) # (inf.(B) .+ sup.(B)) ./ 2 - rB = IntervalArithmetic._sub_round.(mB, inf.(B), RoundUp) + mA = _div_round.(_add_round.(inf.(A), sup.(A), RoundUp), convert(T, 2), RoundUp) # (inf.(A) .+ sup.(A)) ./ 2 + rA = _sub_round.(mA, inf.(A), RoundUp) + mB = _div_round.(_add_round.(inf.(B), sup.(B), RoundUp), convert(T, 2), RoundUp) # (inf.(B) .+ sup.(B)) ./ 2 + rB = _sub_round.(mB, inf.(B), RoundUp) Cinf = zeros(T, size(A, 1), size(B, 2)) Csup = zeros(T, size(A, 1), size(B, 2)) @@ -302,14 +326,14 @@ function __mul(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{Interval{T}}) w Threads.@threads for j ∈ axes(B, 2) for l ∈ axes(A, 2) @inbounds for i ∈ axes(A, 1) - U_ij = IntervalArithmetic._mul_round(abs(mA[i,l]), rB[l,j], RoundUp) - V_ij = IntervalArithmetic._mul_round(rA[i,l], IntervalArithmetic._add_round(abs(mB[l,j]), rB[l,j], RoundUp), RoundUp) - rC_ij = IntervalArithmetic._add_round(U_ij, V_ij, RoundUp) - mAmB_up_ij = IntervalArithmetic._mul_round(mA[i,l], mB[l,j], RoundUp) - mAmB_down_ij = IntervalArithmetic._mul_round(mA[i,l], mB[l,j], RoundDown) - - Cinf[i,j] = IntervalArithmetic._add_round(IntervalArithmetic._sub_round(mAmB_down_ij, rC_ij, RoundDown), Cinf[i,j], RoundDown) - Csup[i,j] = IntervalArithmetic._add_round(IntervalArithmetic._add_round(mAmB_up_ij, rC_ij, RoundUp), Csup[i,j], RoundUp) + U_ij = _mul_round(abs(mA[i,l]), rB[l,j], RoundUp) + V_ij = _mul_round(rA[i,l], _add_round(abs(mB[l,j]), rB[l,j], RoundUp), RoundUp) + rC_ij = _add_round(U_ij, V_ij, RoundUp) + mAmB_up_ij = _mul_round(mA[i,l], mB[l,j], RoundUp) + mAmB_down_ij = _mul_round(mA[i,l], mB[l,j], RoundDown) + + Cinf[i,j] = _add_round(_sub_round(mAmB_down_ij, rC_ij, RoundDown), Cinf[i,j], RoundDown) + Csup[i,j] = _add_round(_add_round(mAmB_up_ij, rC_ij, RoundUp), Csup[i,j], RoundUp) end end end @@ -318,8 +342,8 @@ function __mul(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{Interval{T}}) w end function __mul(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{T}) where {T<:AbstractFloat} - mA = IntervalArithmetic._div_round.(IntervalArithmetic._add_round.(inf.(A), sup.(A), RoundUp), convert(T, 2), RoundUp) # (inf.(A) .+ sup.(A)) ./ 2 - rA = IntervalArithmetic._sub_round.(mA, inf.(A), RoundUp) + mA = _div_round.(_add_round.(inf.(A), sup.(A), RoundUp), convert(T, 2), RoundUp) # (inf.(A) .+ sup.(A)) ./ 2 + rA = _sub_round.(mA, inf.(A), RoundUp) Cinf = zeros(T, size(A, 1), size(B, 2)) Csup = zeros(T, size(A, 1), size(B, 2)) @@ -327,12 +351,12 @@ function __mul(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{T}) where {T<:A Threads.@threads for j ∈ axes(B, 2) for l ∈ axes(A, 2) @inbounds for i ∈ axes(A, 1) - rC_ij = IntervalArithmetic._mul_round(rA[i,l], abs(B[l,j]), RoundUp) - mAmB_up_ij = IntervalArithmetic._mul_round(mA[i,l], B[l,j], RoundUp) - mAmB_down_ij = IntervalArithmetic._mul_round(mA[i,l], B[l,j], RoundDown) + rC_ij = _mul_round(rA[i,l], abs(B[l,j]), RoundUp) + mAmB_up_ij = _mul_round(mA[i,l], B[l,j], RoundUp) + mAmB_down_ij = _mul_round(mA[i,l], B[l,j], RoundDown) - Cinf[i,j] = IntervalArithmetic._add_round(IntervalArithmetic._sub_round(mAmB_down_ij, rC_ij, RoundDown), Cinf[i,j], RoundDown) - Csup[i,j] = IntervalArithmetic._add_round(IntervalArithmetic._add_round(mAmB_up_ij, rC_ij, RoundUp), Csup[i,j], RoundUp) + Cinf[i,j] = _add_round(_sub_round(mAmB_down_ij, rC_ij, RoundDown), Cinf[i,j], RoundDown) + Csup[i,j] = _add_round(_add_round(mAmB_up_ij, rC_ij, RoundUp), Csup[i,j], RoundUp) end end end @@ -341,8 +365,8 @@ function __mul(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{T}) where {T<:A end function __mul(A::AbstractMatrix{T}, B::AbstractMatrix{Interval{T}}) where {T<:AbstractFloat} - mB = IntervalArithmetic._div_round.(IntervalArithmetic._add_round.(inf.(B), sup.(B), RoundUp), convert(T, 2), RoundUp) # (inf.(B) .+ sup.(B)) ./ 2 - rB = IntervalArithmetic._sub_round.(mB, inf.(B), RoundUp) + mB = _div_round.(_add_round.(inf.(B), sup.(B), RoundUp), convert(T, 2), RoundUp) # (inf.(B) .+ sup.(B)) ./ 2 + rB = _sub_round.(mB, inf.(B), RoundUp) Cinf = zeros(T, size(A, 1), size(B, 2)) Csup = zeros(T, size(A, 1), size(B, 2)) @@ -350,12 +374,12 @@ function __mul(A::AbstractMatrix{T}, B::AbstractMatrix{Interval{T}}) where {T<:A Threads.@threads for j ∈ axes(B, 2) for l ∈ axes(A, 2) @inbounds for i ∈ axes(A, 1) - rC_ij = IntervalArithmetic._mul_round(abs(A[i,l]), rB[l,j], RoundUp) - mAmB_up_ij = IntervalArithmetic._mul_round(A[i,l], mB[l,j], RoundUp) - mAmB_down_ij = IntervalArithmetic._mul_round(A[i,l], mB[l,j], RoundDown) + rC_ij = _mul_round(abs(A[i,l]), rB[l,j], RoundUp) + mAmB_up_ij = _mul_round(A[i,l], mB[l,j], RoundUp) + mAmB_down_ij = _mul_round(A[i,l], mB[l,j], RoundDown) - Cinf[i,j] = IntervalArithmetic._add_round(IntervalArithmetic._sub_round(mAmB_down_ij, rC_ij, RoundDown), Cinf[i,j], RoundDown) - Csup[i,j] = IntervalArithmetic._add_round(IntervalArithmetic._add_round(mAmB_up_ij, rC_ij, RoundUp), Csup[i,j], RoundUp) + Cinf[i,j] = _add_round(_sub_round(mAmB_down_ij, rC_ij, RoundDown), Cinf[i,j], RoundDown) + Csup[i,j] = _add_round(_add_round(mAmB_up_ij, rC_ij, RoundUp), Csup[i,j], RoundUp) end end end @@ -364,24 +388,24 @@ function __mul(A::AbstractMatrix{T}, B::AbstractMatrix{Interval{T}}) where {T<:A end function __mul(A::AbstractMatrix{Interval{T}}, B::AbstractVector{Interval{T}}) where {T<:AbstractFloat} - mA = IntervalArithmetic._div_round.(IntervalArithmetic._add_round.(inf.(A), sup.(A), RoundUp), convert(T, 2), RoundUp) # (inf.(A) .+ sup.(A)) ./ 2 - rA = IntervalArithmetic._sub_round.(mA, inf.(A), RoundUp) - mB = IntervalArithmetic._div_round.(IntervalArithmetic._add_round.(inf.(B), sup.(B), RoundUp), convert(T, 2), RoundUp) # (inf.(B) .+ sup.(B)) ./ 2 - rB = IntervalArithmetic._sub_round.(mB, inf.(B), RoundUp) + mA = _div_round.(_add_round.(inf.(A), sup.(A), RoundUp), convert(T, 2), RoundUp) # (inf.(A) .+ sup.(A)) ./ 2 + rA = _sub_round.(mA, inf.(A), RoundUp) + mB = _div_round.(_add_round.(inf.(B), sup.(B), RoundUp), convert(T, 2), RoundUp) # (inf.(B) .+ sup.(B)) ./ 2 + rB = _sub_round.(mB, inf.(B), RoundUp) Cinf = zeros(T, size(A, 1)) Csup = zeros(T, size(A, 1)) Threads.@threads for i ∈ axes(A, 1) @inbounds for l ∈ axes(A, 2) - U_il = IntervalArithmetic._mul_round(abs(mA[i,l]), rB[l], RoundUp) - V_il = IntervalArithmetic._mul_round(rA[i,l], IntervalArithmetic._add_round(abs(mB[l]), rB[l], RoundUp), RoundUp) - rC_il = IntervalArithmetic._add_round(U_il, V_il, RoundUp) - mAmB_up_il = IntervalArithmetic._mul_round(mA[i,l], mB[l], RoundUp) - mAmB_down_il = IntervalArithmetic._mul_round(mA[i,l], mB[l], RoundDown) - - Cinf[i] = IntervalArithmetic._add_round(IntervalArithmetic._sub_round(mAmB_down_il, rC_il, RoundDown), Cinf[i], RoundDown) - Csup[i] = IntervalArithmetic._add_round(IntervalArithmetic._add_round(mAmB_up_il, rC_il, RoundUp), Csup[i], RoundUp) + U_il = _mul_round(abs(mA[i,l]), rB[l], RoundUp) + V_il = _mul_round(rA[i,l], _add_round(abs(mB[l]), rB[l], RoundUp), RoundUp) + rC_il = _add_round(U_il, V_il, RoundUp) + mAmB_up_il = _mul_round(mA[i,l], mB[l], RoundUp) + mAmB_down_il = _mul_round(mA[i,l], mB[l], RoundDown) + + Cinf[i] = _add_round(_sub_round(mAmB_down_il, rC_il, RoundDown), Cinf[i], RoundDown) + Csup[i] = _add_round(_add_round(mAmB_up_il, rC_il, RoundUp), Csup[i], RoundUp) end end @@ -389,20 +413,20 @@ function __mul(A::AbstractMatrix{Interval{T}}, B::AbstractVector{Interval{T}}) w end function __mul(A::AbstractMatrix{Interval{T}}, B::AbstractVector{T}) where {T<:AbstractFloat} - mA = IntervalArithmetic._div_round.(IntervalArithmetic._add_round.(inf.(A), sup.(A), RoundUp), convert(T, 2), RoundUp) # (inf.(A) .+ sup.(A)) ./ 2 - rA = IntervalArithmetic._sub_round.(mA, inf.(A), RoundUp) + mA = _div_round.(_add_round.(inf.(A), sup.(A), RoundUp), convert(T, 2), RoundUp) # (inf.(A) .+ sup.(A)) ./ 2 + rA = _sub_round.(mA, inf.(A), RoundUp) Cinf = zeros(T, size(A, 1)) Csup = zeros(T, size(A, 1)) Threads.@threads for i ∈ axes(A, 1) @inbounds for l ∈ axes(A, 2) - rC_il = IntervalArithmetic._mul_round(rA[i,l], abs(B[l]), RoundUp) - mAB_up_il = IntervalArithmetic._mul_round(mA[i,l], B[l], RoundUp) - mAB_down_il = IntervalArithmetic._mul_round(mA[i,l], B[l], RoundDown) + rC_il = _mul_round(rA[i,l], abs(B[l]), RoundUp) + mAB_up_il = _mul_round(mA[i,l], B[l], RoundUp) + mAB_down_il = _mul_round(mA[i,l], B[l], RoundDown) - Cinf[i] = IntervalArithmetic._add_round(IntervalArithmetic._sub_round(mAB_down_il, rC_il, RoundDown), Cinf[i], RoundDown) - Csup[i] = IntervalArithmetic._add_round(IntervalArithmetic._add_round(mAB_up_il, rC_il, RoundUp), Csup[i], RoundUp) + Cinf[i] = _add_round(_sub_round(mAB_down_il, rC_il, RoundDown), Cinf[i], RoundDown) + Csup[i] = _add_round(_add_round(mAB_up_il, rC_il, RoundUp), Csup[i], RoundUp) end end @@ -410,24 +434,37 @@ function __mul(A::AbstractMatrix{Interval{T}}, B::AbstractVector{T}) where {T<:A end function __mul(A::AbstractMatrix{T}, B::AbstractVector{Interval{T}}) where {T<:AbstractFloat} - mB = IntervalArithmetic._div_round.(IntervalArithmetic._add_round.(inf.(B), sup.(B), RoundUp), convert(T, 2), RoundUp) # (inf.(B) .+ sup.(B)) ./ 2 - rB = IntervalArithmetic._sub_round.(mB, inf.(B), RoundUp) + mB = _div_round.(_add_round.(inf.(B), sup.(B), RoundUp), convert(T, 2), RoundUp) # (inf.(B) .+ sup.(B)) ./ 2 + rB = _sub_round.(mB, inf.(B), RoundUp) Cinf = zeros(T, size(A, 1)) Csup = zeros(T, size(A, 1)) Threads.@threads for i ∈ axes(A, 1) @inbounds for l ∈ axes(A, 2) - rC_il = IntervalArithmetic._mul_round(abs(A[i,l]), rB[l], RoundUp) - AmB_up_il = IntervalArithmetic._mul_round(A[i,l], mB[l], RoundUp) - AmB_down_il = IntervalArithmetic._mul_round(A[i,l], mB[l], RoundDown) + rC_il = _mul_round(abs(A[i,l]), rB[l], RoundUp) + AmB_up_il = _mul_round(A[i,l], mB[l], RoundUp) + AmB_down_il = _mul_round(A[i,l], mB[l], RoundDown) - Cinf[i] = IntervalArithmetic._add_round(IntervalArithmetic._sub_round(AmB_down_il, rC_il, RoundDown), Cinf[i], RoundDown) - Csup[i] = IntervalArithmetic._add_round(IntervalArithmetic._add_round(AmB_up_il, rC_il, RoundUp), Csup[i], RoundUp) + Cinf[i] = _add_round(_sub_round(AmB_down_il, rC_il, RoundDown), Cinf[i], RoundDown) + Csup[i] = _add_round(_add_round(AmB_up_il, rC_il, RoundUp), Csup[i], RoundUp) end end return Cinf, Csup end +# convenient function to propagate NG flag + +function _ensure_ng_flag!(C::AbstractVecOrMat{<:Interval}, ng_flag::Bool) + C .= _unsafe_interval.(getfield.(C, :bareinterval), decoration.(C), ng_flag) + return C +end + +function _ensure_ng_flag!(C::AbstractVecOrMat{<:Complex{<:Interval}}, ng_flag::Bool) + C .= complex.( + _unsafe_interval.(getfield.(real.(C), :bareinterval), decoration.(C), ng_flag), + _unsafe_interval.(getfield.(imag.(C), :bareinterval), decoration.(C), ng_flag) + ) + return C end diff --git a/test/interval_tests/linearalgebra.jl b/test/interval_tests/linearalgebra.jl index 3ed90c0e0..e0152fc53 100644 --- a/test/interval_tests/linearalgebra.jl +++ b/test/interval_tests/linearalgebra.jl @@ -8,6 +8,7 @@ import LinearAlgebra end @testset "Matrix multiplication" begin + IntervalArithmetic.matmul_mode() = IntervalArithmetic.MatMulMode{:fast}() A = [interval(2, 4) interval(-2, 1) ; interval(-1, 2) interval(2, 4)] imA = interval(im) * A @@ -18,4 +19,5 @@ end @test all(isequal_interval.(imA * imA, -[interval(-2, 19.5) interval(-16, 10) ; interval(-10, 16) interval(-2, 19.5)])) @test all(isequal_interval.(mid.(A) * imA, interval(im)*[interval(5, 12.5) interval(-8, 2) ; interval(-2, 8) interval(5, 12.5)])) @test all(isequal_interval.(imA * mid.(A), interval(im)*[interval(5, 12.5) interval(-8, 2) ; interval(-2, 8) interval(5, 12.5)])) + IntervalArithmetic.matmul_mode() = IntervalArithmetic.MatMulMode{:slow}() end