Make a LinearAlgebra extension and use the improved fast matrix multiplication

Author: OlivierHnt

Project.toml

@@ -7,19 +7,22 @@ version = "0.22.23"
 CRlibm = "96374032-68de-5a5b-8d9e-752f78720389"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+OpenBLASConsistentFPCSR_jll = "6cdc7f73-28fd-5e50-80fb-958a8875b1af"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RoundingEmulator = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705"
 
 [weakdeps]
 DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 
 [extensions]
 IntervalArithmeticDiffRulesExt = "DiffRules"
 IntervalArithmeticForwardDiffExt = "ForwardDiff"
 IntervalArithmeticIntervalSetsExt = "IntervalSets"
+IntervalArithmeticLinearAlgebraExt = "LinearAlgebra"
 IntervalArithmeticRecipesBaseExt = "RecipesBase"
 
 [compat]
@@ -29,6 +32,7 @@ ForwardDiff = "0.10"
 IntervalSets = "0.7"
 LinearAlgebra = "1.10"
 MacroTools = "0.5"
+OpenBLASConsistentFPCSR_jll = "0.3.29"
 Printf = "1.10"
 RecipesBase = "1"
 RoundingEmulator = "0.2"

ext/IntervalArithmeticLinearAlgebraExt.jl

@@ -0,0 +1,114 @@
+module IntervalArithmeticLinearAlgebraExt
+
+using IntervalArithmetic
+import LinearAlgebra
+
+# contructor for `UniformScaling`
+
+IntervalArithmetic.interval(::Type{T}, J::LinearAlgebra.UniformScaling, d::IntervalArithmetic.Decoration = com; format::Symbol = :infsup) where {T} =
+    LinearAlgebra.UniformScaling(interval(T, J.λ, d; format = format))
+IntervalArithmetic.interval(J::LinearAlgebra.UniformScaling, d::IntervalArithmetic.Decoration = com; format::Symbol = :infsup) =
+    LinearAlgebra.UniformScaling(interval(J.λ, d; format = format))
+
+# by-pass generic `opnorm` from LinearAlgebra to prevent NG flag
+
+function LinearAlgebra.opnorm1(A::AbstractMatrix{T}) where {T<:RealOrComplexI}
+    LinearAlgebra.require_one_based_indexing(A)
+    m, n = size(A)
+    Tnorm = typeof(float(real(zero(T))))
+    Tsum = promote_type(Float64, Tnorm)
+    nrm = zero(Tsum)
+    @inbounds begin
+        for j = 1:n
+            nrmj = zero(Tsum)
+            for i = 1:m
+                nrmj += LinearAlgebra.norm(A[i,j])
+            end
+            nrm = max(nrm, nrmj)
+        end
+    end
+    return convert(Tnorm, nrm)
+end
+
+function LinearAlgebra.opnormInf(A::AbstractMatrix{T}) where {T<:RealOrComplexI}
+    LinearAlgebra.require_one_based_indexing(A)
+    m, n = size(A)
+    Tnorm = typeof(float(real(zero(T))))
+    Tsum = promote_type(Float64, Tnorm)
+    nrm = zero(Tsum)
+    @inbounds begin
+        for i = 1:m
+            nrmi = zero(Tsum)
+            for j = 1:n
+                nrmi += LinearAlgebra.norm(A[i,j])
+            end
+            nrm = max(nrm, nrmi)
+        end
+    end
+    return convert(Tnorm, nrm)
+end
+
+# matrix eigenvalues
+
+function LinearAlgebra.eigvals!(A::AbstractMatrix{<:Interval}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby)
+    # note: this function does not overwrite `A`
+    v = _eigvals(A, permute, scale, sortby)
+    isreal(v) && return v
+    _fold_conjugate!(v)
+    isreal(v) && return real(v)
+    return v
+end
+
+LinearAlgebra.eigvals!(A::AbstractMatrix{<:Complex{<:Interval}}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby) =
+    # note: this function does not overwrite `A`
+    _eigvals(A, permute, scale, sortby)
+
+function _eigvals(A, permute, scale, sortby)
+    # Gershgorin circle theorem
+    B = _similarity_transform(A, permute, scale, sortby)
+    v = LinearAlgebra.diag(B)
+    T = eltype(v)
+    for j ∈ axes(B, 1)
+        r = zero(T)
+        for i ∈ axes(B, 2)
+            if i ≠ j
+                r += abs(B[i,j])
+            end
+        end
+        v[j] = interval(v[j], r; format = :midpoint)
+    end
+    return v
+end
+
+function _similarity_transform(A, permute, scale, sortby)
+    mA = mid.(A)
+    mλ, mV = LinearAlgebra.eigen(mA; permute = permute, scale = scale, sortby = sortby)
+    mλ .+= LinearAlgebra.diag(mV \ (mA * mV - mV * LinearAlgebra.Diagonal(mλ)))
+    Λ = LinearAlgebra.Diagonal(interval(mλ))
+    V = interval(mV)
+    V .= Λ .+ inv(V) * (A * V - V * Λ)
+    return V
+end
+
+function _fold_conjugate!(v)
+    for i ∈ eachindex(v)
+        vᵢ = v[i]
+        idxs = findall(j -> (j ≠ i) & !isdisjoint_interval(conj(vᵢ), v[j]), eachindex(v))
+        if isempty(idxs)
+            v[i] = real(vᵢ)
+        else
+            w = view(v, idxs)
+            z = conj(intersect_interval(conj(vᵢ), reduce(intersect_interval, w)))
+            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))))
+            v[i] = z
+        end
+    end
+    return v
+end
+
+# matrix determinant
+
+LinearAlgebra.det(A::AbstractMatrix{<:Interval}) = real(reduce(*, LinearAlgebra.eigvals(A)))
+LinearAlgebra.det(A::AbstractMatrix{<:Complex{<:Interval}}) = reduce(*, LinearAlgebra.eigvals(A))
+
+end

src/IntervalArithmetic.jl

@@ -84,6 +84,14 @@ include("display.jl")
 #
 
 import LinearAlgebra
+import OpenBLASConsistentFPCSR_jll # 32-bit systems are not supported
+
+if Int != Int32
+    # use the same number of threads as the default BLAS library
+    ccall((:openblas_set_num_threads64_, OpenBLASConsistentFPCSR_jll.libopenblas),
+        Cint, (Cint,),
+        LinearAlgebra.BLAS.get_num_threads())
+end
 
 include("matmul.jl")