From 17da336ff6437c9b743fcf7d9e4fe243f13b3f82 Mon Sep 17 00:00:00 2001
From: Marcelo Forets <mforets@gmail.com>
Date: Sun, 10 Oct 2021 09:02:35 -0400
Subject: [PATCH 1/5] add pencil types

---
 docs/src/lib/types.md   |  12 +++++
 src/IntervalMatrices.jl |   2 +
 src/pencils.jl          | 114 ++++++++++++++++++++++++++++++++++++++++
 test/pencils.jl         |  25 +++++++++
 test/runtests.jl        |   1 +
 5 files changed, 154 insertions(+)
 create mode 100644 src/pencils.jl
 create mode 100644 test/pencils.jl

diff --git a/docs/src/lib/types.md b/docs/src/lib/types.md
index 64348cce..8b69a9d3 100644
--- a/docs/src/lib/types.md
+++ b/docs/src/lib/types.md
@@ -28,3 +28,15 @@ IntervalMatrix
 ```@docs
 IntervalMatrixPower
 ```
+
+## Matrix pencil
+
+```@docs
+IntervalMatrixPencil
+```
+
+## Affine interval matrix
+
+```@docs
+AffineIntervalMatrix
+```
diff --git a/src/IntervalMatrices.jl b/src/IntervalMatrices.jl
index 53dd7f08..00430fc2 100644
--- a/src/IntervalMatrices.jl
+++ b/src/IntervalMatrices.jl
@@ -13,6 +13,7 @@ using Reexport
 # Type defining an interval matrix
 # ================================
 include("matrix.jl")
+include("pencils.jl")
 
 # =================================
 # Operations for interval matrices
@@ -61,6 +62,7 @@ export IntervalMatrixPower,
        base,
        index
 
+export IntervalMatrixPencil, AffineIntervalMatrix
 
 set_multiplication_mode(config[:multiplication])
 
diff --git a/src/pencils.jl b/src/pencils.jl
new file mode 100644
index 00000000..4c47d904
--- /dev/null
+++ b/src/pencils.jl
@@ -0,0 +1,114 @@
+"""
+    IntervalMatrixPencil{T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}} <: AbstractIntervalMatrix{IT}
+
+Interval matrix representing the matrix
+
+```math
+A₀ + λA₁,
+```
+where ``A₀`` and ``A₁`` are real (or complex) matrices, and ``λ`` is an interval.
+
+### Fields
+
+- `A0` -- matrix
+- `A1` -- matrix
+- `λ`  -- interval
+
+### Examples
+
+The matrix pencil ``I + [1 1; -1 1] * (0 .. 1)`` is:
+
+```jldoctest
+julia> using LinearAlgebra
+
+julia> P = IntervalMatrixPencil(Matrix(1.0I, 2, 2), [1 1; -1 1.], 0 .. 1);
+
+julia> P
+2×2 IntervalMatrixPencil{Float64, Interval{Float64}, Matrix{Float64}, Matrix{Float64}}:
+  [1, 2]  [0, 1]
+ [-1, 0]  [1, 2]
+```
+"""
+struct IntervalMatrixPencil{T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}} <: AbstractIntervalMatrix{IT}
+    A0::MT0
+    A1::MT1
+    λ::IT
+
+    # inner constructor with dimension check
+    function IntervalMatrixPencil(A0::MT0, A1::MT1, λ::IT) where {T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}}
+        @assert checksquare(A0) == checksquare(A1) "the size of `A0` and `A1` should match, got $(size(A0)) and $(size(A1)) respectively"
+        return new{T, IT, MT0, MT1}(A0, A1, λ)
+    end
+end
+
+IndexStyle(::Type{<:IntervalMatrixPencil}) = IndexLinear()
+size(M::IntervalMatrixPencil) = size(M.A0)
+getindex(M::IntervalMatrixPencil, i::Int) = getindex(M.A0, i) + M.λ * getindex(M.A1, i)
+function setindex!(M::IntervalMatrixPencil{T}, X::T, inds...) where {T}
+    setindex!(M.A0, X, inds...)
+    if !iszero(getindex(M.A1, inds...)) || !iszero(M.λ)
+        setindex!(M.A1, zero(T), inds...)
+    end
+end
+copy(M::IntervalMatrixPencil) = IntervalMatrixPencil(copy(M.A0), copy(M.A1), M.λ)
+
+"""
+    AffineIntervalMatrix{T, IT, MT0<:AbstractMatrix{T}, MT<:AbstractMatrix{T}, MTA<:AbstractVector{MT}} <: AbstractIntervalMatrix{IT}
+
+Interval matrix representing the matrix
+
+```math
+A₀ + λ₁A₁ + λ₂A₂ + … + λₖAₖ,
+```
+where ``A₀`` and ``A₁, …, Aₖ`` are real (or complex) matrices, and ``λ₁, …, λₖ``
+are intervals.
+
+### Fields
+
+- `A0` -- matrix
+- `A`  -- vector of matrices
+- `λ`  -- vector of intervals
+
+### Notes
+
+This type is the general case of the [`IntervalMatrixPencil`](@ref), which only
+contains one matrix proportional to an interval.
+
+### Examples
+
+The matrix pencil ``I + [1 1; -1 1] * (0 .. 1) + [0 1; 1 0] * (2 .. 3)`` is:
+
+```jldoctest
+julia> using LinearAlgebra
+
+julia> A0 = Matrix(1.0I, 2, 2);
+
+julia> A1 = [1 1; -1 1.]; A2 = [0 1; 1 0];
+
+julia> λ1 = 0 .. 1; λ2 = 2 .. 3;
+
+julia> P = AffineIntervalMatrix(A0, [A1, A2], [λ1, λ1]);
+
+julia> P[1, 1]
+[1, 3]
+
+julia> P[1, 2]
+[0, 2]
+```
+"""
+struct AffineIntervalMatrix{T, IT, MT0<:AbstractMatrix{T}, MT<:AbstractMatrix{T}, MTA<:AbstractVector{MT}, VIT<:AbstractVector{IT}} <: AbstractIntervalMatrix{IT}
+    A0::MT0
+    A::MTA
+    λ::VIT
+end
+
+IndexStyle(::Type{<:AffineIntervalMatrix}) = IndexLinear()
+size(M::AffineIntervalMatrix) = size(M.A0)
+getindex(M::AffineIntervalMatrix, i::Int) = getindex(M.A0, i) + sum(M.λ[k] * getindex(M.A[k], i) for k in eachindex(M.λ))
+function setindex!(M::AffineIntervalMatrix{T}, X::T, inds...) where {T}
+    setindex!(M.A0, X, inds...)
+    @inbounds for k in 1:length(M.A)
+        setindex!(M.A[k], zero(T), inds...)
+    end
+end
+copy(M::AffineIntervalMatrix) = AffineIntervalMatrix(copy(M.A0), copy(M.A), M.λ)
diff --git a/test/pencils.jl b/test/pencils.jl
new file mode 100644
index 00000000..a6a09f4c
--- /dev/null
+++ b/test/pencils.jl
@@ -0,0 +1,25 @@
+@testset "Interval matrix pencil" begin
+    A0 = [1 0; 0 1.]
+    A1 = [0 1; 1 0.]
+    λ = 0 .. 1
+    P = IntervalMatrixPencil(A0, A1, λ)
+
+    @test size(P) == (2, 2)
+    @test P[1, 1] == interval(1)
+    @test P[1, 2] == interval(0, 1)
+    P[1, 1] = 2.0
+    @test P[1, 1] == interval(2)
+end
+
+@testset "Affine interval matrix" begin
+    A0 = [1 0; 0 1.]
+    A1 = [0 1; 1 0.]; A2 = copy(A1)
+    λ1 = 0 .. 1; λ2 = copy(λ1)
+    P = AffineIntervalMatrix(A0, [A1, A2], [λ1, λ2])
+
+    @test size(P) == (2, 2)
+    @test P[1, 1] == interval(1)
+    @test P[1, 2] == interval(0, 2)
+    P[1, 1] = 5.0
+    @test P[1, 1] == interval(5)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 78a497ff..54de463c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,3 +11,4 @@ include("constructors.jl")
 include("arithmetic.jl")
 include("setops.jl")
 include("exponential.jl")
+include("pencils.jl")

From c5951da986d0f1e4994440305c427daa0814db97 Mon Sep 17 00:00:00 2001
From: Marcelo Forets <mforets@gmail.com>
Date: Sun, 10 Oct 2021 16:26:07 -0400
Subject: [PATCH 2/5] Update src/pencils.jl

Co-authored-by: Christian Schilling <git@christianschilling.net>
---
 src/pencils.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pencils.jl b/src/pencils.jl
index 4c47d904..18f3b503 100644
--- a/src/pencils.jl
+++ b/src/pencils.jl
@@ -76,7 +76,7 @@ contains one matrix proportional to an interval.
 
 ### Examples
 
-The matrix pencil ``I + [1 1; -1 1] * (0 .. 1) + [0 1; 1 0] * (2 .. 3)`` is:
+The affine matrix ``I + [1 1; -1 1] * (0 .. 1) + [0 1; 1 0] * (2 .. 3)`` is:
 
 ```jldoctest
 julia> using LinearAlgebra

From b7412c173cb3a47f67ec9b81d4b019736249ada2 Mon Sep 17 00:00:00 2001
From: Marcelo Forets <mforets@gmail.com>
Date: Sun, 10 Oct 2021 17:47:25 -0400
Subject: [PATCH 3/5] fix test

---
 src/pencils.jl | 17 +++++++++++++----
 1 file changed, 13 insertions(+), 4 deletions(-)

diff --git a/src/pencils.jl b/src/pencils.jl
index 18f3b503..7a5339e0 100644
--- a/src/pencils.jl
+++ b/src/pencils.jl
@@ -46,9 +46,7 @@ size(M::IntervalMatrixPencil) = size(M.A0)
 getindex(M::IntervalMatrixPencil, i::Int) = getindex(M.A0, i) + M.λ * getindex(M.A1, i)
 function setindex!(M::IntervalMatrixPencil{T}, X::T, inds...) where {T}
     setindex!(M.A0, X, inds...)
-    if !iszero(getindex(M.A1, inds...)) || !iszero(M.λ)
-        setindex!(M.A1, zero(T), inds...)
-    end
+    setindex!(M.A1, zero(T), inds...)
 end
 copy(M::IntervalMatrixPencil) = IntervalMatrixPencil(copy(M.A0), copy(M.A1), M.λ)
 
@@ -90,7 +88,7 @@ julia> λ1 = 0 .. 1; λ2 = 2 .. 3;
 julia> P = AffineIntervalMatrix(A0, [A1, A2], [λ1, λ1]);
 
 julia> P[1, 1]
-[1, 3]
+[1, 2]
 
 julia> P[1, 2]
 [0, 2]
@@ -100,6 +98,17 @@ struct AffineIntervalMatrix{T, IT, MT0<:AbstractMatrix{T}, MT<:AbstractMatrix{T}
     A0::MT0
     A::MTA
     λ::VIT
+
+    # inner constructor with dimension check
+    function AffineIntervalMatrix(A0::MT0, A::MTA, λ::VIT) where {T, IT, MT0<:AbstractMatrix{T}, MT<:AbstractMatrix{T}, MTA<:AbstractVector{MT}, VIT<:AbstractVector{IT}}
+        k = length(A)
+        @assert k == length(λ) "expected `A` and `λ` to have the same length, got $(length(A)) and $(length(λ)) respectively"
+        n = checksquare(A0)
+        @inbounds for i in 1:k
+            @assert n == checksquare(A[i]) "each matrix should have the same size $n × $n"
+        end
+        return new{T, IT, MT0, MT, MTA, VIT}(A0, A, λ)
+    end
 end
 
 IndexStyle(::Type{<:AffineIntervalMatrix}) = IndexLinear()

From f41bb876d0649869effbec5719fef73d479d5297 Mon Sep 17 00:00:00 2001
From: Marcelo Forets <mforets@gmail.com>
Date: Thu, 14 Oct 2021 14:47:50 -0400
Subject: [PATCH 4/5] update names

---
 docs/src/lib/types.md          |  7 +------
 src/IntervalMatrices.jl        |  5 +++--
 src/{pencils.jl => affine.jl}  | 24 ++++++++++++------------
 test/{pencils.jl => affine.jl} |  6 +++---
 test/runtests.jl               |  2 +-
 5 files changed, 20 insertions(+), 24 deletions(-)
 rename src/{pencils.jl => affine.jl} (74%)
 rename test/{pencils.jl => affine.jl} (76%)

diff --git a/docs/src/lib/types.md b/docs/src/lib/types.md
index 8b69a9d3..4354e3f7 100644
--- a/docs/src/lib/types.md
+++ b/docs/src/lib/types.md
@@ -29,14 +29,9 @@ IntervalMatrix
 IntervalMatrixPower
 ```
 
-## Matrix pencil
-
-```@docs
-IntervalMatrixPencil
-```
-
 ## Affine interval matrix
 
 ```@docs
+AffineIntervalMatrix1
 AffineIntervalMatrix
 ```
diff --git a/src/IntervalMatrices.jl b/src/IntervalMatrices.jl
index 00430fc2..0b256264 100644
--- a/src/IntervalMatrices.jl
+++ b/src/IntervalMatrices.jl
@@ -13,7 +13,7 @@ using Reexport
 # Type defining an interval matrix
 # ================================
 include("matrix.jl")
-include("pencils.jl")
+include("affine.jl")
 
 # =================================
 # Operations for interval matrices
@@ -62,7 +62,8 @@ export IntervalMatrixPower,
        base,
        index
 
-export IntervalMatrixPencil, AffineIntervalMatrix
+export AffineIntervalMatrix,
+       AffineIntervalMatrix1
 
 set_multiplication_mode(config[:multiplication])
 
diff --git a/src/pencils.jl b/src/affine.jl
similarity index 74%
rename from src/pencils.jl
rename to src/affine.jl
index 7a5339e0..eb1ec792 100644
--- a/src/pencils.jl
+++ b/src/affine.jl
@@ -1,5 +1,5 @@
 """
-    IntervalMatrixPencil{T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}} <: AbstractIntervalMatrix{IT}
+    AffineIntervalMatrix1{T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}} <: AbstractIntervalMatrix{IT}
 
 Interval matrix representing the matrix
 
@@ -16,39 +16,39 @@ where ``A₀`` and ``A₁`` are real (or complex) matrices, and ``λ`` is an int
 
 ### Examples
 
-The matrix pencil ``I + [1 1; -1 1] * (0 .. 1)`` is:
+The matrix ``I + [1 1; -1 1] * (0 .. 1)`` is:
 
 ```jldoctest
 julia> using LinearAlgebra
 
-julia> P = IntervalMatrixPencil(Matrix(1.0I, 2, 2), [1 1; -1 1.], 0 .. 1);
+julia> P = AffineIntervalMatrix1(Matrix(1.0I, 2, 2), [1 1; -1 1.], 0 .. 1);
 
 julia> P
-2×2 IntervalMatrixPencil{Float64, Interval{Float64}, Matrix{Float64}, Matrix{Float64}}:
+2×2 AffineIntervalMatrix1{Float64, Interval{Float64}, Matrix{Float64}, Matrix{Float64}}:
   [1, 2]  [0, 1]
  [-1, 0]  [1, 2]
 ```
 """
-struct IntervalMatrixPencil{T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}} <: AbstractIntervalMatrix{IT}
+struct AffineIntervalMatrix1{T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}} <: AbstractIntervalMatrix{IT}
     A0::MT0
     A1::MT1
     λ::IT
 
     # inner constructor with dimension check
-    function IntervalMatrixPencil(A0::MT0, A1::MT1, λ::IT) where {T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}}
+    function AffineIntervalMatrix1(A0::MT0, A1::MT1, λ::IT) where {T, IT, MT0<:AbstractMatrix{T}, MT1<:AbstractMatrix{T}}
         @assert checksquare(A0) == checksquare(A1) "the size of `A0` and `A1` should match, got $(size(A0)) and $(size(A1)) respectively"
         return new{T, IT, MT0, MT1}(A0, A1, λ)
     end
 end
 
-IndexStyle(::Type{<:IntervalMatrixPencil}) = IndexLinear()
-size(M::IntervalMatrixPencil) = size(M.A0)
-getindex(M::IntervalMatrixPencil, i::Int) = getindex(M.A0, i) + M.λ * getindex(M.A1, i)
-function setindex!(M::IntervalMatrixPencil{T}, X::T, inds...) where {T}
+IndexStyle(::Type{<:AffineIntervalMatrix1}) = IndexLinear()
+size(M::AffineIntervalMatrix1) = size(M.A0)
+getindex(M::AffineIntervalMatrix1, i::Int) = getindex(M.A0, i) + M.λ * getindex(M.A1, i)
+function setindex!(M::AffineIntervalMatrix1{T}, X::T, inds...) where {T}
     setindex!(M.A0, X, inds...)
     setindex!(M.A1, zero(T), inds...)
 end
-copy(M::IntervalMatrixPencil) = IntervalMatrixPencil(copy(M.A0), copy(M.A1), M.λ)
+copy(M::AffineIntervalMatrix1) = AffineIntervalMatrix1(copy(M.A0), copy(M.A1), M.λ)
 
 """
     AffineIntervalMatrix{T, IT, MT0<:AbstractMatrix{T}, MT<:AbstractMatrix{T}, MTA<:AbstractVector{MT}} <: AbstractIntervalMatrix{IT}
@@ -69,7 +69,7 @@ are intervals.
 
 ### Notes
 
-This type is the general case of the [`IntervalMatrixPencil`](@ref), which only
+This type is the general case of the [`AffineIntervalMatrix1`](@ref), which only
 contains one matrix proportional to an interval.
 
 ### Examples
diff --git a/test/pencils.jl b/test/affine.jl
similarity index 76%
rename from test/pencils.jl
rename to test/affine.jl
index a6a09f4c..e649c5d1 100644
--- a/test/pencils.jl
+++ b/test/affine.jl
@@ -1,8 +1,8 @@
-@testset "Interval matrix pencil" begin
+@testset "Affine interval matrix in one interval" begin
     A0 = [1 0; 0 1.]
     A1 = [0 1; 1 0.]
     λ = 0 .. 1
-    P = IntervalMatrixPencil(A0, A1, λ)
+    P = AffineIntervalMatrix1(A0, A1, λ)
 
     @test size(P) == (2, 2)
     @test P[1, 1] == interval(1)
@@ -11,7 +11,7 @@
     @test P[1, 1] == interval(2)
 end
 
-@testset "Affine interval matrix" begin
+@testset "Affine interval matrix in several intervals" begin
     A0 = [1 0; 0 1.]
     A1 = [0 1; 1 0.]; A2 = copy(A1)
     λ1 = 0 .. 1; λ2 = copy(λ1)
diff --git a/test/runtests.jl b/test/runtests.jl
index 54de463c..9910562a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,4 +11,4 @@ include("constructors.jl")
 include("arithmetic.jl")
 include("setops.jl")
 include("exponential.jl")
-include("pencils.jl")
+include("affine.jl")

From c6fb550c4b8eacf3a7ae6ab9d0d94626b91b0a13 Mon Sep 17 00:00:00 2001
From: Marcelo Forets <mforets@gmail.com>
Date: Sat, 16 Oct 2021 08:03:17 -0400
Subject: [PATCH 5/5] add test with complex matrix

---
 test/affine.jl | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/test/affine.jl b/test/affine.jl
index e649c5d1..c9f45fff 100644
--- a/test/affine.jl
+++ b/test/affine.jl
@@ -9,6 +9,11 @@
     @test P[1, 2] == interval(0, 1)
     P[1, 1] = 2.0
     @test P[1, 1] == interval(2)
+
+    # complex interval
+    λc = λ + im*λ/2
+    Pc = AffineIntervalMatrix1(A0, A1, λc)
+    @test eltype(Pc) == Complex{Interval{Float64}}
 end
 
 @testset "Affine interval matrix in several intervals" begin
@@ -22,4 +27,9 @@ end
     @test P[1, 2] == interval(0, 2)
     P[1, 1] = 5.0
     @test P[1, 1] == interval(5)
+
+    # complex interval
+    λ1c = λ1 + im*λ1/2; λ2c = λ2 - im*λ2/2
+    Pc = AffineIntervalMatrix(A0, [A1, A2], [λ1c, λ2c])
+    @test eltype(Pc) == Complex{Interval{Float64}}
 end