From 12957e0f93f9e2347e80f36a8cf80df20a2e82f7 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Fri, 2 Sep 2022 16:35:59 +0900 Subject: [PATCH 01/16] Implemented Theorem 7 in Miya14 --- docs/src/api/eigenvalues.md | 7 ++ docs/src/references.md | 28 ++++++ src/IntervalLinearAlgebra.jl | 1 + src/eigenvalues/miyajima_svd.jl | 105 +++++++++++++++++++++ test/runtests.jl | 1 + test/test_eigenvalues/test_miyajima_svd.jl | 12 +++ 6 files changed, 154 insertions(+) create mode 100644 src/eigenvalues/miyajima_svd.jl create mode 100644 test/test_eigenvalues/test_miyajima_svd.jl diff --git a/docs/src/api/eigenvalues.md b/docs/src/api/eigenvalues.md index 424dcadb..887a1250 100644 --- a/docs/src/api/eigenvalues.md +++ b/docs/src/api/eigenvalues.md @@ -20,3 +20,10 @@ Pages=["verify_eigs.jl"] Private=false ``` +## Singular Values Verification + +```@autodocs +Modules=[IntervalLinearAlgebra] +Pages=["miyajima_svd.jl"] +Private=false +``` \ No newline at end of file diff --git a/docs/src/references.md b/docs/src/references.md index 1415abab..69f9c5b6 100644 --- a/docs/src/references.md +++ b/docs/src/references.md @@ -104,6 +104,34 @@ L. Jaulin and B. Desrochers, [*Introduction to the algebra of separators with ap ``` --- + +#### [MIYA14] + +```@raw html + +``` +--- + + #### [NEU90] ```@raw html diff --git a/src/IntervalLinearAlgebra.jl b/src/IntervalLinearAlgebra.jl index a210653f..91de4ff2 100644 --- a/src/IntervalLinearAlgebra.jl +++ b/src/IntervalLinearAlgebra.jl @@ -42,6 +42,7 @@ include("pils/pils_solvers.jl") include("eigenvalues/interval_eigenvalues.jl") include("eigenvalues/verify_eigs.jl") +include("eigenvalues/miyajima_svd.jl") function __init__() @require IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153" include("linear_systems/oettli_nonlinear.jl") diff --git a/src/eigenvalues/miyajima_svd.jl b/src/eigenvalues/miyajima_svd.jl new file mode 100644 index 00000000..7dd6e333 --- /dev/null +++ b/src/eigenvalues/miyajima_svd.jl @@ -0,0 +1,105 @@ +# In this file we implement some of the algorithms from +# +# Shinya Miyajima +# +# Verified bounds for all the singular values of matrix +# Japan J. Indust. Appl. Math. +# DOI 10.1007/s13160-014-0145-5 + +abstract type AbstractIntervalSVDSolver end + +struct M1 <: AbstractIntervalSVDSolver end + +# The method is called svdbox to mirror eigenbox, even if it would be best to call it +# svd interval + +using LinearAlgebra + +function _power_iteration_singularvalue(A, max_iter) + n = size(A, 2) + xp = rand(n) + @inbounds for _ in 1:max_iter + xp .= A'*(A*xp) + xp ./= norm(xp) + end + return xp +end + + +function _bound_perron_frobenius_singularvalue(M, max_iter=10) + + size(M) == (1, 1) && return M[1] + xpf = IA.Interval.(_power_iteration_singularvalue(mid.(M), max_iter)) + Mxpf = M'* (M * xpf) + ρ = zero(eltype(M)) + @inbounds for (i, xi) in enumerate(xpf) + iszero(xi) && continue + tmp = Mxpf[i] / xi + ρ = max(ρ, tmp.hi) + end + return ρ +end + + +""" + svdbox(A[, method = M1()]) + +Return an enclosure for all the singular values of `A`. + +# Input +- `A` -- interval matrix +- `method` -- the method used to compute the enclosure + Possible values are + - M1 Algorithm from Theorem 7 in [[MIYA14]](@ref) + +### Algorithm + +The algorithms used by the function are described in [[MIYA14]](@ref). + +### Example +The matrix encloses a matrix that has singular values +``3, √5, 2, 0``. + +```julia +julia> A = [0.9..1.1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0] +4×5 Matrix{Interval{Float64}}: + [0.899999, 1.10001] [0, 0] [0, 0] [0, 0] [2, 2] + [0, 0] [0, 0] [3, 3] [0, 0] [0, 0] + [0, 0] [0, 0] [0, 0] [0, 0] [0, 0] + [0, 0] [2, 2] [0, 0] [0, 0] [0, 0] + + julia> IntervalLinearAlgebra.svdbox(A, IntervalLinearAlgebra.M1()) + 4-element Vector{Interval{Float64}}: + [2.98999, 3.01001] + [2.22606, 2.24607] + [1.98999, 2.01001] + [-0.0100001, 0.0100001] +``` + +""" +function svdbox(A::AbstractMatrix{Interval{T}}, ::M1) where T + mA = mid.(A) + svdA = svd(mA) + U = Interval{T}.(svdA.U) + Vt = Interval{T}.(svdA.Vt) + Σ = diagm(Interval{T}.(svdA.S)) + V = Interval{T}.(svdA.V) + + E = U*Σ*Vt - A + normE = _bound_perron_frobenius_singularvalue(E) + + F = Vt*V-I + normF = _bound_perron_frobenius_singularvalue(F) + + G = U'*U-I + normG = _bound_perron_frobenius_singularvalue(G) + + @assert normF < 1 "It is not possible to verify the singular values with this precision" + @assert normG < 1 "It is not possible to verify the singular values with this precision" + + svdbounds = [hull(σ*sqrt((1-normF)*(1-normG))-normE, + σ*sqrt((1+normF)*(1+normG))+normE) + for σ in diag(Σ)] + + return svdbounds +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index de035dc7..52f0974d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,6 +10,7 @@ include("test_utils.jl") include("test_eigenvalues/test_interval_eigenvalues.jl") include("test_eigenvalues/test_verify_eigs.jl") +include("test_eigenvalues/test_miyajima_svd.jl") include("test_solvers/test_enclosures.jl") include("test_solvers/test_epsilon_inflation.jl") include("test_solvers/test_precondition.jl") diff --git a/test/test_eigenvalues/test_miyajima_svd.jl b/test/test_eigenvalues/test_miyajima_svd.jl new file mode 100644 index 00000000..d2df7905 --- /dev/null +++ b/test/test_eigenvalues/test_miyajima_svd.jl @@ -0,0 +1,12 @@ +@testset "verify singular value perron frobenius " begin + A = [0.5 0.5;0.3 0.7] + ρ = IntervalLinearAlgebra._bound_perron_frobenius_singularvalue(A) + @test 1 ≤ ρ +end + +@testset "verified svd" begin + A = [0.9..1.1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0] + Σ = IntervalLinearAlgebra.svdbox(A, IntervalLinearAlgebra.M1()) + + @test all([3 ∈ Σ[1], 5 ∈ Σ[2]^2, 2 ∈ Σ[3], 0 ∈ Σ[4]]) +end From 6a863c8ba71a99ac69fc8eacead2edc6dc261327 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Fri, 2 Sep 2022 17:01:06 +0900 Subject: [PATCH 02/16] Corrected bound on the 2-norm of E, F and G --- src/eigenvalues/miyajima_svd.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/eigenvalues/miyajima_svd.jl b/src/eigenvalues/miyajima_svd.jl index 7dd6e333..06d558b5 100644 --- a/src/eigenvalues/miyajima_svd.jl +++ b/src/eigenvalues/miyajima_svd.jl @@ -25,7 +25,9 @@ function _power_iteration_singularvalue(A, max_iter) return xp end - +# we use this function to bound the $2$ norm of a matrix, +# as remarked in Rump, Verified bounds for singular values, +# Equation 3.3 function _bound_perron_frobenius_singularvalue(M, max_iter=10) size(M) == (1, 1) && return M[1] @@ -77,7 +79,7 @@ julia> A = [0.9..1.1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0] ``` """ -function svdbox(A::AbstractMatrix{Interval{T}}, ::M1) where T +function svdbox(A::AbstractMatrix{Interval{T}}, method = ::M1) where T mA = mid.(A) svdA = svd(mA) U = Interval{T}.(svdA.U) @@ -86,13 +88,13 @@ function svdbox(A::AbstractMatrix{Interval{T}}, ::M1) where T V = Interval{T}.(svdA.V) E = U*Σ*Vt - A - normE = _bound_perron_frobenius_singularvalue(E) + normE = sqrt(_bound_perron_frobenius_singularvalue(E)) F = Vt*V-I - normF = _bound_perron_frobenius_singularvalue(F) + normF = sqrt(_bound_perron_frobenius_singularvalue(F)) G = U'*U-I - normG = _bound_perron_frobenius_singularvalue(G) + normG = sqrt(_bound_perron_frobenius_singularvalue(G)) @assert normF < 1 "It is not possible to verify the singular values with this precision" @assert normG < 1 "It is not possible to verify the singular values with this precision" From dc90cd9450359651dd3b09f3c11227842a740e5d Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Fri, 2 Sep 2022 18:14:43 +0900 Subject: [PATCH 03/16] Added Algorithm R1 from Miya14 --- src/eigenvalues/miyajima_svd.jl | 50 ++++++++++++++++++---- test/test_eigenvalues/test_miyajima_svd.jl | 4 ++ 2 files changed, 46 insertions(+), 8 deletions(-) diff --git a/src/eigenvalues/miyajima_svd.jl b/src/eigenvalues/miyajima_svd.jl index 06d558b5..8410905f 100644 --- a/src/eigenvalues/miyajima_svd.jl +++ b/src/eigenvalues/miyajima_svd.jl @@ -9,6 +9,7 @@ abstract type AbstractIntervalSVDSolver end struct M1 <: AbstractIntervalSVDSolver end +struct R1 <: AbstractIntervalSVDSolver end # The method is called svdbox to mirror eigenbox, even if it would be best to call it # svd interval @@ -53,6 +54,7 @@ Return an enclosure for all the singular values of `A`. - `method` -- the method used to compute the enclosure Possible values are - M1 Algorithm from Theorem 7 in [[MIYA14]](@ref) + - R1 Rump Algorithm from Theorem 5 in [[MIYA14]](@ref) ### Algorithm @@ -72,14 +74,14 @@ julia> A = [0.9..1.1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0] julia> IntervalLinearAlgebra.svdbox(A, IntervalLinearAlgebra.M1()) 4-element Vector{Interval{Float64}}: - [2.98999, 3.01001] - [2.22606, 2.24607] - [1.98999, 2.01001] - [-0.0100001, 0.0100001] + [2.89999, 3.10001] + [2.13606, 2.33607] + [1.89999, 2.10001] + [-0.100001, 0.100001] ``` """ -function svdbox(A::AbstractMatrix{Interval{T}}, method = ::M1) where T +function svdbox(A::AbstractMatrix{Interval{T}}, ::M1) where T mA = mid.(A) svdA = svd(mA) U = Interval{T}.(svdA.U) @@ -88,13 +90,13 @@ function svdbox(A::AbstractMatrix{Interval{T}}, method = ::M1) where T V = Interval{T}.(svdA.V) E = U*Σ*Vt - A - normE = sqrt(_bound_perron_frobenius_singularvalue(E)) + normE = sqrt(_bound_perron_frobenius_singularvalue(abs.(E))) F = Vt*V-I - normF = sqrt(_bound_perron_frobenius_singularvalue(F)) + normF = sqrt(_bound_perron_frobenius_singularvalue(abs.(F))) G = U'*U-I - normG = sqrt(_bound_perron_frobenius_singularvalue(G)) + normG = sqrt(_bound_perron_frobenius_singularvalue(abs.(G))) @assert normF < 1 "It is not possible to verify the singular values with this precision" @assert normG < 1 "It is not possible to verify the singular values with this precision" @@ -104,4 +106,36 @@ function svdbox(A::AbstractMatrix{Interval{T}}, method = ::M1) where T for σ in diag(Σ)] return svdbounds +end + +function certifysvd(A::AbstractMatrix{Interval{T}}, U, V, Vt) where {T} + IU = Interval{T}.(U) + IV = Interval{T}.(V) + IVt = Interval{T}.(Vt) + + F = IVt*IV-I + G = (IU)'*IU-I + normF = sqrt(_bound_perron_frobenius_singularvalue(abs.(F))) + normG = sqrt(_bound_perron_frobenius_singularvalue(abs.(G))) + + @assert normF < 1 "It is not possible to verify the singular values with this precision" + @assert normG < 1 "It is not possible to verify the singular values with this precision" + + M = IU'*A*IV + D = diag(M) + E = M-diagm(D) + + normE = sqrt(_bound_perron_frobenius_singularvalue(abs.(E))) + + svdbounds = [hull((σ-normE)/sqrt((1+normF)*(1+normG)), + (σ+normE)/sqrt((1-normF)*(1-normG))) + for σ in D] + + return sort!(svdbounds, rev = true) +end + +function svdbox(A::AbstractMatrix{Interval{T}}, ::R1) where T + mA = mid.(A) + svdA = svd(mA) + return certifysvd(A, svdA.U, svdA.V, svdA.Vt) end \ No newline at end of file diff --git a/test/test_eigenvalues/test_miyajima_svd.jl b/test/test_eigenvalues/test_miyajima_svd.jl index d2df7905..5a166660 100644 --- a/test/test_eigenvalues/test_miyajima_svd.jl +++ b/test/test_eigenvalues/test_miyajima_svd.jl @@ -9,4 +9,8 @@ end Σ = IntervalLinearAlgebra.svdbox(A, IntervalLinearAlgebra.M1()) @test all([3 ∈ Σ[1], 5 ∈ Σ[2]^2, 2 ∈ Σ[3], 0 ∈ Σ[4]]) + + Σ = IntervalLinearAlgebra.svdbox(A, IntervalLinearAlgebra.R1()) + @test all([3 ∈ Σ[1], 5 ∈ Σ[2]^2, 2 ∈ Σ[3], 0 ∈ Σ[4]]) + end From 5fc049f2fa76a945a08ea2732a20ddf5e675e5d5 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Fri, 2 Sep 2022 18:47:04 +0900 Subject: [PATCH 04/16] Missing absolute value in R1 --- src/eigenvalues/miyajima_svd.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/eigenvalues/miyajima_svd.jl b/src/eigenvalues/miyajima_svd.jl index 8410905f..70fd58d0 100644 --- a/src/eigenvalues/miyajima_svd.jl +++ b/src/eigenvalues/miyajima_svd.jl @@ -127,8 +127,8 @@ function certifysvd(A::AbstractMatrix{Interval{T}}, U, V, Vt) where {T} normE = sqrt(_bound_perron_frobenius_singularvalue(abs.(E))) - svdbounds = [hull((σ-normE)/sqrt((1+normF)*(1+normG)), - (σ+normE)/sqrt((1-normF)*(1-normG))) + svdbounds = [hull((abs(σ)-normE)/sqrt((1+normF)*(1+normG)), + (abs(σ)+normE)/sqrt((1-normF)*(1-normG))) for σ in D] return sort!(svdbounds, rev = true) From 76299e4625a88f735e3964c33eefc7774d1c873a Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Sun, 5 Feb 2023 17:04:17 +0900 Subject: [PATCH 05/16] Fixed complex matrix multiplication --- src/IntervalLinearAlgebra.jl | 5 ++++ src/eigenvalues/miyajima_svd.jl | 43 +++++++++++++++++++++++++++++++-- src/multiplication.jl | 22 +++++++++++++++-- src/utils.jl | 2 +- 4 files changed, 67 insertions(+), 5 deletions(-) diff --git a/src/IntervalLinearAlgebra.jl b/src/IntervalLinearAlgebra.jl index 91de4ff2..14e116ee 100644 --- a/src/IntervalLinearAlgebra.jl +++ b/src/IntervalLinearAlgebra.jl @@ -51,4 +51,9 @@ end set_multiplication_mode(config[:multiplication]) +using LinearAlgebra +BLAS.set_num_threads(1) +@warn "The number of BLAS threads was set to 1 to ensure rounding mode is consistent" + + end diff --git a/src/eigenvalues/miyajima_svd.jl b/src/eigenvalues/miyajima_svd.jl index 70fd58d0..9b7aaa16 100644 --- a/src/eigenvalues/miyajima_svd.jl +++ b/src/eigenvalues/miyajima_svd.jl @@ -31,6 +31,10 @@ end # Equation 3.3 function _bound_perron_frobenius_singularvalue(M, max_iter=10) + if any(M.<0) + @warn "The matrix has nonpositive entries; M'*M could be nonnegative, but careful" + end + size(M) == (1, 1) && return M[1] xpf = IA.Interval.(_power_iteration_singularvalue(mid.(M), max_iter)) Mxpf = M'* (M * xpf) @@ -43,9 +47,9 @@ function _bound_perron_frobenius_singularvalue(M, max_iter=10) return ρ end - +export svdbox """ - svdbox(A[, method = M1()]) + svdbox(A[, method = R1()]) Return an enclosure for all the singular values of `A`. @@ -81,6 +85,8 @@ julia> A = [0.9..1.1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0] ``` """ +svdbox(A) = svdbox(A, R1()) + function svdbox(A::AbstractMatrix{Interval{T}}, ::M1) where T mA = mid.(A) svdA = svd(mA) @@ -134,8 +140,41 @@ function certifysvd(A::AbstractMatrix{Interval{T}}, U, V, Vt) where {T} return sort!(svdbounds, rev = true) end +function certifysvd(A::AbstractMatrix{Complex{Interval{T}}}, U, V, Vt) where {T} + IU = Interval{T}.(real(U))+im*Interval{T}.(imag(U)) + IV = Interval{T}.(real(V))+im*Interval{T}.(imag(V)) + IVt = Interval{T}.(real(Vt))+im*Interval{T}.(imag(Vt)) + + F = IVt*IV-I + G = (IU)'*IU-I + normF = sqrt(_bound_perron_frobenius_singularvalue(abs.(F))) + normG = sqrt(_bound_perron_frobenius_singularvalue(abs.(G))) + + @assert normF < 1 "It is not possible to verify the singular values with this precision" + @assert normG < 1 "It is not possible to verify the singular values with this precision" + + M = IU'*A*IV + D = diag(M) + E = M-diagm(D) + + normE = sqrt(_bound_perron_frobenius_singularvalue(abs.(E))) + + svdbounds = [hull((abs(σ)-normE)/sqrt((1+normF)*(1+normG)), + (abs(σ)+normE)/sqrt((1-normF)*(1-normG))) + for σ in D] + + return sort!(svdbounds, rev = true) +end + + function svdbox(A::AbstractMatrix{Interval{T}}, ::R1) where T mA = mid.(A) svdA = svd(mA) return certifysvd(A, svdA.U, svdA.V, svdA.Vt) +end + +function svdbox(A::AbstractMatrix{Complex{Interval{T}}}, ::R1) where T + mA = mid.(A) + svdA = svd(mA) + return certifysvd(A, svdA.U, svdA.V, svdA.Vt) end \ No newline at end of file diff --git a/src/multiplication.jl b/src/multiplication.jl index b4342be6..4233e32e 100644 --- a/src/multiplication.jl +++ b/src/multiplication.jl @@ -26,18 +26,36 @@ Sets the algorithm used to perform matrix multiplication with interval matrices. """ function set_multiplication_mode(multype) type = MultiplicationType{multype}() + + # real matrices @eval *(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{Interval{T}}) where T = *($type, A, B) @eval *(A::AbstractMatrix{T}, B::AbstractMatrix{Interval{T}}) where T = *($type, A, B) - + @eval *(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{T}) where T = *($type, A, B) - + @eval *(A::Diagonal, B::AbstractMatrix{Interval{T}}) where T = *($type, A, B) + @eval *(A::AbstractMatrix{Interval{T}}, B::Diagonal) where T = *($type, A, B) + config[:multiplication] = multype end +function *(A::AbstractMatrix{Complex{Interval{T}}}, B::AbstractMatrix) where T + return real(A)*B+im*imag(A)*B +end + +function *(A::AbstractMatrix, B::AbstractMatrix{Complex{Interval{T}}}) where T + return A*real(B)+im*A*imag(B) +end + +function *(A::AbstractMatrix{Complex{Interval{T}}}, B::AbstractMatrix{Complex{Interval{T}}}) where T + rA, iA = real(A), imag(A) + rB, iB = real(B), imag(B) + return rA*rB-iA*iB+im*(iA*rB+rA*iB) +end + function *(::MultiplicationType{:slow}, A, B) TS = promote_type(eltype(A), eltype(B)) return mul!(similar(B, TS, (size(A,1), size(B,2))), A, B) diff --git a/src/utils.jl b/src/utils.jl index 6bc44e4c..4fd13d96 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,4 +1,4 @@ -mid(x::Real) = x +# mid(x::Real) = x #now exported by IntervalArithmetic mid(x::Complex) = x """ From 70dd6745a61199ab54b58082b8a9d9edca2a0578 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Thu, 9 Feb 2023 15:06:14 +0900 Subject: [PATCH 06/16] Added complex test --- test/test_eigenvalues/test_miyajima_svd.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/test_eigenvalues/test_miyajima_svd.jl b/test/test_eigenvalues/test_miyajima_svd.jl index 5a166660..daf2535d 100644 --- a/test/test_eigenvalues/test_miyajima_svd.jl +++ b/test/test_eigenvalues/test_miyajima_svd.jl @@ -13,4 +13,8 @@ end Σ = IntervalLinearAlgebra.svdbox(A, IntervalLinearAlgebra.R1()) @test all([3 ∈ Σ[1], 5 ∈ Σ[2]^2, 2 ∈ Σ[3], 0 ∈ Σ[4]]) + A = im*A + + Σ = IntervalLinearAlgebra.svdbox(A, IntervalLinearAlgebra.R1()) + @test all([3 ∈ Σ[1], 5 ∈ Σ[2]^2, 2 ∈ Σ[3], 0 ∈ Σ[4]]) end From 561d1c396fd99f3d486359ba8613882b3ca20998 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Thu, 9 Feb 2023 15:34:36 +0900 Subject: [PATCH 07/16] Cheery picked right commits --- Project.toml | 1 + src/IntervalLinearAlgebra.jl | 14 ++++++++++++++ 2 files changed, 15 insertions(+) diff --git a/Project.toml b/Project.toml index 5e76dfaf..b0c83124 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.1.5" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +OpenBLASConsistentFPCSR_jll = "6cdc7f73-28fd-5e50-80fb-958a8875b1af" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/src/IntervalLinearAlgebra.jl b/src/IntervalLinearAlgebra.jl index a210653f..4eaa2054 100644 --- a/src/IntervalLinearAlgebra.jl +++ b/src/IntervalLinearAlgebra.jl @@ -43,9 +43,23 @@ include("pils/pils_solvers.jl") include("eigenvalues/interval_eigenvalues.jl") include("eigenvalues/verify_eigs.jl") +using LinearAlgebra + +if Sys.ARCH == :x86_64 + using OpenBLASConsistentFPCSR_jll +end + function __init__() @require IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153" include("linear_systems/oettli_nonlinear.jl") @require LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" include("linear_systems/oettli_linear.jl") + if Sys.ARCH == :x86_64 + @info "Switching to OpenBLAS with ConsistentFPCSR = 1 flag enabled, guarantees + coherent floating point rounding mode over all threads" + BLAS.lbt_forward(OpenBLASConsistentFPCSR_jll.libopenblas_path; verbose = true) + else + BLAS.set_num_threads(1) + @warn "The number of BLAS threads was set to 1 to ensure rounding mode is consistent" + end end set_multiplication_mode(config[:multiplication]) From 579dcc6bfaac11ce3ee7033ff9ed98228f0367f0 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Thu, 9 Feb 2023 16:29:58 +0900 Subject: [PATCH 08/16] Added methods to currectly dispatch complex mul --- src/multiplication.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/multiplication.jl b/src/multiplication.jl index b4342be6..36c1f545 100644 --- a/src/multiplication.jl +++ b/src/multiplication.jl @@ -38,6 +38,21 @@ function set_multiplication_mode(multype) config[:multiplication] = multype end +function *(A::AbstractMatrix{Complex{Interval{T}}}, B::AbstractMatrix) where T + return real(A)*B+im*imag(A)*B +end + +function *(A::AbstractMatrix, B::AbstractMatrix{Complex{Interval{T}}}) where T + return A*real(B)+im*A*imag(B) +end + +function *(A::AbstractMatrix{Complex{Interval{T}}}, B::AbstractMatrix{Complex{Interval{T}}}) where T + rA, iA = real(A), imag(A) + rB, iB = real(B), imag(B) + return rA*rB-iA*iB+im*(iA*rB+rA*iB) +end + + function *(::MultiplicationType{:slow}, A, B) TS = promote_type(eltype(A), eltype(B)) return mul!(similar(B, TS, (size(A,1), size(B,2))), A, B) From fe5b3883f11598618ca2a8fc610e80c98ca84160 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Thu, 9 Feb 2023 16:46:54 +0900 Subject: [PATCH 09/16] Added methods to dispatch correctly cpx matrix mul --- src/multiplication.jl | 8 ++++++++ test/test_multiplication.jl | 8 ++++++++ 2 files changed, 16 insertions(+) diff --git a/src/multiplication.jl b/src/multiplication.jl index 36c1f545..21b9c584 100644 --- a/src/multiplication.jl +++ b/src/multiplication.jl @@ -52,6 +52,14 @@ function *(A::AbstractMatrix{Complex{Interval{T}}}, B::AbstractMatrix{Complex{In return rA*rB-iA*iB+im*(iA*rB+rA*iB) end +function *(A::AbstractMatrix{Complex{T}}, B::AbstractMatrix{Interval{T}}) where {T} + return real(A)*B+im*imag(A)*B +end + +function *(A::AbstractMatrix{Interval{T}}, B::AbstractMatrix{Complex{T}}) where {T} + return A*real(B)+im*A*imag(B) +end + function *(::MultiplicationType{:slow}, A, B) TS = promote_type(eltype(A), eltype(B)) diff --git a/test/test_multiplication.jl b/test/test_multiplication.jl index d7c9b9ee..caea00ae 100644 --- a/test/test_multiplication.jl +++ b/test/test_multiplication.jl @@ -4,14 +4,22 @@ @test get_multiplication_mode() == Dict(:multiplication => :fast) A = [2..4 -2..1; -1..2 2..4] + imA = im*A + set_multiplication_mode(:slow) @test A * A == [0..18 -16..8; -8..16 0..18] + @test imA * imA == -1*[0..18 -16..8; -8..16 0..18] set_multiplication_mode(:rank1) @test A * A == [0..18 -16..8; -8..16 0..18] + @test imA * imA == -1*[0..18 -16..8; -8..16 0..18] set_multiplication_mode(:fast) @test A * A == [-2..19.5 -16..10; -10..16 -2..19.5] @test A * mid.(A) == [5..12.5 -8..2; -2..8 5..12.5] @test mid.(A) * A == [5..12.5 -8..2; -2..8 5..12.5] + + @test imA * imA == -1*[-2..19.5 -16..10; -10..16 -2..19.5] + @test mid.(A) * imA == im*[5..12.5 -8..2; -2..8 5..12.5] + @test imA * mid.(A) == im*[5..12.5 -8..2; -2..8 5..12.5] end From 0681fde3ab33974a2d1f2b71f533ebbf42828bd9 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Thu, 16 Feb 2023 20:10:25 +0900 Subject: [PATCH 10/16] Update src/IntervalLinearAlgebra.jl Co-authored-by: Luca Ferranti <49938764+lucaferranti@users.noreply.github.com> --- src/IntervalLinearAlgebra.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/IntervalLinearAlgebra.jl b/src/IntervalLinearAlgebra.jl index 7e9ac63f..d21c73e4 100644 --- a/src/IntervalLinearAlgebra.jl +++ b/src/IntervalLinearAlgebra.jl @@ -53,7 +53,7 @@ function __init__() @require LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" include("linear_systems/oettli_linear.jl") if Sys.ARCH == :x86_64 @info "Switching to OpenBLAS with ConsistentFPCSR = 1 flag enabled, guarantees - coherent floating point rounding mode over all threads" + correct floating point rounding mode over all threads." BLAS.lbt_forward(OpenBLASConsistentFPCSR_jll.libopenblas_path; verbose = true) else BLAS.set_num_threads(1) From 50a8aba601f185b1305d81d0dc585d98c8ae4ecc Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Sat, 18 Feb 2023 13:17:25 +0900 Subject: [PATCH 11/16] Added function that tests that OpenBLAS is working --- Project.toml | 6 ++- src/IntervalLinearAlgebra.jl | 20 ++++++++++ src/numerical_test/multithread.jl | 38 +++++++++++++++++++ test/runtests.jl | 2 +- .../test_numerical_test.jl | 9 +++++ 5 files changed, 73 insertions(+), 2 deletions(-) create mode 100644 src/numerical_test/multithread.jl create mode 100644 test/test_numerical_test/test_numerical_test.jl diff --git a/Project.toml b/Project.toml index eabc4382..88e8e079 100644 --- a/Project.toml +++ b/Project.toml @@ -4,12 +4,16 @@ authors = ["Luca Ferranti"] version = "0.1.6" [deps] +BLISBLAS = "6f275bd8-fec0-4d39-945b-7e95a765fa1e" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" +Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OpenBLASConsistentFPCSR_jll = "6cdc7f73-28fd-5e50-80fb-958a8875b1af" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +SetRoundingLLVM = "259bb760-bc25-474b-814a-89aa0a9f335c" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] @@ -24,4 +28,4 @@ julia = "1.6" [extras] IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153" -LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" \ No newline at end of file +LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" diff --git a/src/IntervalLinearAlgebra.jl b/src/IntervalLinearAlgebra.jl index d21c73e4..3d4dc0c0 100644 --- a/src/IntervalLinearAlgebra.jl +++ b/src/IntervalLinearAlgebra.jl @@ -42,10 +42,15 @@ include("pils/pils_solvers.jl") include("eigenvalues/interval_eigenvalues.jl") include("eigenvalues/verify_eigs.jl") +include("numerical_test/multithread.jl") + using LinearAlgebra if Sys.ARCH == :x86_64 using OpenBLASConsistentFPCSR_jll +else + @warn "The behaviour of multithreaded OpenBlas on this architecture is unclear, + we will import MKL" end function __init__() @@ -55,9 +60,24 @@ function __init__() @info "Switching to OpenBLAS with ConsistentFPCSR = 1 flag enabled, guarantees correct floating point rounding mode over all threads." BLAS.lbt_forward(OpenBLASConsistentFPCSR_jll.libopenblas_path; verbose = true) + + N = BLAS.get_num_threads() + K = 1024 + if NumericalTest.rounding_test(N, K) + @info "OpenBLAS is giving correct rounding on a ($K,$K) test matrix on $N threads" + else + @warn "OpenBLAS is not rounding correctly on the test matrix" + @warn "The number of BLAS threads was set to 1 to ensure rounding mode is consistent" + if !NumericalTest.rounding_test(1, K) + @warn "The rounding test failed on 1 thread" + end + end else BLAS.set_num_threads(1) @warn "The number of BLAS threads was set to 1 to ensure rounding mode is consistent" + if !NumericalTest.rounding_test(1, K) + @warn "The rounding test failed on 1 thread" + end end end diff --git a/src/numerical_test/multithread.jl b/src/numerical_test/multithread.jl new file mode 100644 index 00000000..b80e2125 --- /dev/null +++ b/src/numerical_test/multithread.jl @@ -0,0 +1,38 @@ +module NumericalTest + +export rounding_test + +function test_matrix(k) + A = zeros(Float64, (k, k)) + A[:, end] = fill(2^(-53), k) + for i in 1:k-1 + A[i,i] = 1.0 + end + return A +end + +using LinearAlgebra, SetRoundingLLVM +""" + rounding_test(n, k) + +Let `u=fill(2^(-53), k-1)` and let A be the matrix +[I u; +0 2^(-53)] + +This test checks the result of A*A' in different rounding modes, +running BLAS on `n` threads +""" +function rounding_test(n,k) + + + BLAS.set_num_threads( n ) + A = test_matrix( k ) + llvm_setrounding(RoundUp) + B = BLAS.gemm('N', 'T', 1.0, A, A) + llvm_setrounding(RoundNearest) + + return all([B[i,i]==nextfloat(1.0) for i in 1:k-1]) +end + + +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index de035dc7..4d635af9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,7 @@ include("test_classify.jl") include("test_multiplication.jl") include("test_utils.jl") - +include("test_numerical_test/test_numerical_test.jl") include("test_eigenvalues/test_interval_eigenvalues.jl") include("test_eigenvalues/test_verify_eigs.jl") include("test_solvers/test_enclosures.jl") diff --git a/test/test_numerical_test/test_numerical_test.jl b/test/test_numerical_test/test_numerical_test.jl new file mode 100644 index 00000000..70710e3d --- /dev/null +++ b/test/test_numerical_test/test_numerical_test.jl @@ -0,0 +1,9 @@ +@testset "Test Numerical Test" + A = IntervalLinearAlgebra.NumericalTest.test_matrix(4) + @test A == [1.0 0 0 2^(-53); + 0 1.0 0 2^(-53); + 0 0 1.0 2^(-53); + 0 0 0 2^(-53)] + + @test IntervalLinearAlgebra.NumericalTest.rounding_test(4) == true +end \ No newline at end of file From b155e357195fc4b80a47bd38315f14f534bcf177 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Sat, 18 Feb 2023 13:52:46 +0900 Subject: [PATCH 12/16] Removed a wrong BLAS set number of threads --- src/IntervalLinearAlgebra.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/IntervalLinearAlgebra.jl b/src/IntervalLinearAlgebra.jl index e9d75b46..90105653 100644 --- a/src/IntervalLinearAlgebra.jl +++ b/src/IntervalLinearAlgebra.jl @@ -85,8 +85,5 @@ end set_multiplication_mode(config[:multiplication]) using LinearAlgebra -BLAS.set_num_threads(1) -@warn "The number of BLAS threads was set to 1 to ensure rounding mode is consistent" - end From 9fbcb4c8151d4159fe804ad17f2ea4b21857f22b Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Sat, 18 Feb 2023 14:07:47 +0900 Subject: [PATCH 13/16] Cleaned dependencies --- Project.toml | 3 --- src/IntervalLinearAlgebra.jl | 2 +- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 88e8e079..3a256f42 100644 --- a/Project.toml +++ b/Project.toml @@ -4,12 +4,9 @@ authors = ["Luca Ferranti"] version = "0.1.6" [deps] -BLISBLAS = "6f275bd8-fec0-4d39-945b-7e95a765fa1e" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" -Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OpenBLASConsistentFPCSR_jll = "6cdc7f73-28fd-5e50-80fb-958a8875b1af" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" diff --git a/src/IntervalLinearAlgebra.jl b/src/IntervalLinearAlgebra.jl index 90105653..3b1d1d2c 100644 --- a/src/IntervalLinearAlgebra.jl +++ b/src/IntervalLinearAlgebra.jl @@ -51,7 +51,7 @@ if Sys.ARCH == :x86_64 using OpenBLASConsistentFPCSR_jll else @warn "The behaviour of multithreaded OpenBlas on this architecture is unclear, - we will import MKL" + we will fallback to single threaded OpenBLAS" end function __init__() From aa7a07f2aea9fef42fe5e0b9e56785f5af468cd2 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Sat, 18 Feb 2023 16:12:38 +0900 Subject: [PATCH 14/16] Removed dependencies --- Project.toml | 4 ---- src/numerical_test/multithread.jl | 10 +++++----- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 88e8e079..bb5d28d4 100644 --- a/Project.toml +++ b/Project.toml @@ -4,16 +4,12 @@ authors = ["Luca Ferranti"] version = "0.1.6" [deps] -BLISBLAS = "6f275bd8-fec0-4d39-945b-7e95a765fa1e" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" -Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OpenBLASConsistentFPCSR_jll = "6cdc7f73-28fd-5e50-80fb-958a8875b1af" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" -SetRoundingLLVM = "259bb760-bc25-474b-814a-89aa0a9f335c" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] diff --git a/src/numerical_test/multithread.jl b/src/numerical_test/multithread.jl index b80e2125..5c01e024 100644 --- a/src/numerical_test/multithread.jl +++ b/src/numerical_test/multithread.jl @@ -11,7 +11,7 @@ function test_matrix(k) return A end -using LinearAlgebra, SetRoundingLLVM +using LinearAlgebra """ rounding_test(n, k) @@ -27,10 +27,10 @@ function rounding_test(n,k) BLAS.set_num_threads( n ) A = test_matrix( k ) - llvm_setrounding(RoundUp) - B = BLAS.gemm('N', 'T', 1.0, A, A) - llvm_setrounding(RoundNearest) - + B = setrounding(Float64, RoundUp) do + BLAS.gemm('N', 'T', 1.0, A, A) + end + return all([B[i,i]==nextfloat(1.0) for i in 1:k-1]) end From 7d8c333ae7118c5fee347b0c8aa2d0f145e219f7 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Sat, 18 Feb 2023 16:40:03 +0900 Subject: [PATCH 15/16] Better test in Numerical Test testset --- test/test_numerical_test/test_numerical_test.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/test_numerical_test/test_numerical_test.jl b/test/test_numerical_test/test_numerical_test.jl index 70710e3d..d5e79107 100644 --- a/test/test_numerical_test/test_numerical_test.jl +++ b/test/test_numerical_test/test_numerical_test.jl @@ -1,9 +1,11 @@ -@testset "Test Numerical Test" +@testset "Test Numerical Test" begin A = IntervalLinearAlgebra.NumericalTest.test_matrix(4) @test A == [1.0 0 0 2^(-53); 0 1.0 0 2^(-53); 0 0 1.0 2^(-53); 0 0 0 2^(-53)] - @test IntervalLinearAlgebra.NumericalTest.rounding_test(4) == true + # we test the singlethread version, so that CI points + # out if setting rounding modes is broken + @test IntervalLinearAlgebra.NumericalTest.rounding_test(1, 2) == true end \ No newline at end of file From 2ea3de1cc5301b75f92689088494e29deeb102e8 Mon Sep 17 00:00:00 2001 From: Isaia Nisoli Date: Sat, 18 Feb 2023 16:50:21 +0900 Subject: [PATCH 16/16] Corrected missing variable --- src/IntervalLinearAlgebra.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/IntervalLinearAlgebra.jl b/src/IntervalLinearAlgebra.jl index 3d4dc0c0..975af56d 100644 --- a/src/IntervalLinearAlgebra.jl +++ b/src/IntervalLinearAlgebra.jl @@ -75,7 +75,7 @@ function __init__() else BLAS.set_num_threads(1) @warn "The number of BLAS threads was set to 1 to ensure rounding mode is consistent" - if !NumericalTest.rounding_test(1, K) + if !NumericalTest.rounding_test(1, 1024) @warn "The rounding test failed on 1 thread" end end