mixed float interval matrix multiplication (#83)

Author: lucaferranti

src/multiplication.jl

@@ -26,9 +26,15 @@ Sets the algorithm used to perform matrix multiplication with interval matrices.
 """
 function set_multiplication_mode(multype)
     type = MultiplicationType{multype}()
-    @eval *(A::AbstractMatrix{Interval{T}} where T, B::AbstractMatrix{Interval{T}} where T) =
+    @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
 
@@ -66,6 +72,58 @@ function *(::MultiplicationType{:fast},
     return Interval.(Cinf, Csup)
 end
 
+
+function *(::MultiplicationType{:fast},
+    A::AbstractMatrix{T},
+    B::AbstractMatrix{Interval{T}}) where {T<:Real}
+
+    Binf = inf.(B)
+    Bsup = sup.(B)
+
+    mB, R, Csup = setrounding(T, RoundUp) do
+        mB = Binf + 0.5 * (Bsup - Binf)
+
+        rB = mB - Binf
+
+        R = abs.(A) * rB
+        Csup = A * mB + R
+
+        return mB, R, Csup
+    end
+
+    Cinf = setrounding(T, RoundDown) do
+        A * mB - R
+    end
+
+    return Interval.(Cinf, Csup)
+end
+
+function *(::MultiplicationType{:fast},
+    A::AbstractMatrix{Interval{T}},
+    B::AbstractMatrix{T}) where {T<:Real}
+
+    Ainf = inf.(A)
+    Asup = sup.(A)
+
+    mA, R, Csup = setrounding(T, RoundUp) do
+        mA = Ainf + 0.5 * (Asup - Ainf)
+
+        rA = mA - Ainf
+
+        R = rA * abs.(B)
+        Csup = mA * B + R
+
+        return mA, R, Csup
+    end
+
+    Cinf = setrounding(T, RoundDown) do
+        mA * B - R
+    end
+
+    return Interval.(Cinf, Csup)
+end
+
+
 function *(::MultiplicationType{:rank1},
            A::AbstractMatrix{Interval{T}},
            B::AbstractMatrix{Interval{T}}) where {T<:Real}

test/test_multiplication.jl

@@ -5,11 +5,13 @@
 
     A = [2..4 -2..1; -1..2 2..4]
     set_multiplication_mode(:slow)
-    @test A*A == [0..18 -16..8; -8..16 0..18]
+    @test A * A == [0..18 -16..8; -8..16 0..18]
 
     set_multiplication_mode(:rank1)
-    @test A*A == [0..18 -16..8; -8..16 0..18]
+    @test A * A == [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 * 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]
 end