Use 3-arg mul! for Companion and improve test coverage (#69)

* Use 3-arg mul! for Companion

* Test Companion

* Add test for AbstracVector

* Test isassigned, just to be sure

* Use `begin:end`

Co-authored-by: Mark Kittisopikul <[email protected]>

* Test size(C,3) == 1 just to be sure.

* Use begin:end-1

Co-authored-by: Mark Kittisopikul <[email protected]>

* Add jldoctest

Co-authored-by: Mark Kittisopikul <[email protected]>

* Use Base.require_one_based_indexing

* Remove unneeded T

* Use jldoctest

Co-authored-by: Mark Kittisopikul <[email protected]>

* Better docstring

* Add Polynomials for docs

* Fix jldoctest output

* Add OffsetArray test

* Use test/Project.toml

* Remove unneeded "isassigned" import from Base

* Add @test_throws for OffsetMatrix * Companion

* Use convert

Co-authored-by: Mark Kittisopikul <[email protected]>

* Revert to using collect

* Bump version to 3.0.0

* Update test/companion.jl to use no_offset_view

Co-authored-by: Mark Kittisopikul <[email protected]>

* Update src/companion.jl to use convert again

Co-authored-by: Mark Kittisopikul <[email protected]>

* Add import

Co-authored-by: Mark Kittisopikul <[email protected]>

Project.toml

@@ -1,6 +1,6 @@
 name = "SpecialMatrices"
 uuid = "928aab9d-ef52-54ac-8ca1-acd7ca42c160"
-version = "2.0.1"
+version = "3.0.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -9,9 +9,3 @@ Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 [compat]
 Polynomials = "1, 2, 3"
 julia = "1.6"
-
-[extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = ["Test"]

docs/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 
 [compat]
 Documenter = "0.27.3"

src/SpecialMatrices.jl

@@ -5,7 +5,7 @@ using Polynomials
 
 import LinearAlgebra: Matrix, inv, det, ishermitian, isposdef, mul!
 
-import Base: getindex, isassigned, size, *
+import Base: getindex, size, *
 import Base.\
 
 include("cauchy.jl") # Cauchy matrix

src/companion.jl

@@ -1,100 +1,108 @@
 export Companion
+using LinearAlgebra: dot
 
 """
-[`Companion` matrix](http://en.wikipedia.org/wiki/Companion_matrix)
+   Companion(v::Union{AbstractVector,Polynomial})::AbstractMatrix
 
-```julia
-julia> A=Companion([3,2,1])
-3x3 Companion{Int64}:
+Construct a lazy
+[`companion` matrix](http://en.wikipedia.org/wiki/Companion_matrix)
+from the coefficients of its characteristic (monic) polynomial.
+
+The matrix is `n × n` for a vector input of length `n`
+or an input `Polynomial` of degree `n`,
+but the storage here is only `O(n)`.
+This version puts the coefficients
+along the last column of the matrix.
+Some texts put the coefficients along the first row,
+the transpose of the convention used here.
+
+This type has efficient methods
+for `mul!` and `inv`.
+
+```jldoctest
+julia> A = Companion([3,2,1])
+3×3 Companion{Int64}:
  0  0  -3
  1  0  -2
  0  1  -1
 ```
 Also, directly from a polynomial:
 
-```julia
+```jldoctest
 julia> using Polynomials
 
-julia> P=Polynomial([2.0,3,4,5])
-Polynomial(2 + 3x + 4x^2 + 5x^3)
+julia> P = Polynomial(2:5)
+Polynomials.Polynomial(2 + 3*x + 4*x^2 + 5*x^3)
 
-julia> C=Companion(P)
+julia> C = Companion(P)
 3×3 Companion{Float64}:
  0.0  0.0  -0.4
  1.0  0.0  -0.6
  0.0  1.0  -0.8
 ```
 """
-struct Companion{T} <: AbstractArray{T, 2}
+struct Companion{T} <: AbstractMatrix{T}
     c :: Vector{T}
 end
 
-# Generate companion matrix from a polynomial
+Companion(v::AbstractVector{T}) where T = Companion{T}(v)
 
-function Companion(P::Polynomial{T}) where T
-   n = length(P)
-   c = Array{T}(undef, n-1)
-   d=P.coeffs[n]
-   for i=1:n-1
-       c[i]=P.coeffs[i]/d
-   end
-   Companion(c)
+# Construct companion matrix from a polynomial
+
+function Companion(P::Polynomial)
+    c = P.coeffs[begin:end-1] ./ P.coeffs[end]
+    return Companion(c)
 end
 
-#Basic property computations
-size(C::Companion, r::Int) = (r==1 || r==2) ? length(C.c) :
-    throw(ArgumentError("Companion is of rank 2"))
+# Basic properties
 
 function size(C::Companion)
     n = length(C.c)
-    n, n
+    return n, n
 end
 
-#XXX Inefficient but works
-# getindex(C::Companion, i, j) = getindex(Matrix(C), i, j)
-# isassigned(C::Companion, i, j) = isassigned(Matrix(C), i, j)
-getindex(C::Companion{T}, i::Int, j::Int) where T = (j==length(C.c)) ? -C.c[i] : (i==j+1 ? one(T) : zero(T) )
-isassigned(C::Companion, i::Int, j::Int) = (j==length(C.c)) ? isassigned(C.c,i) : true
+@inline Base.@propagate_inbounds function getindex(
+    C::Companion{T},
+    i::Int,
+    j::Int,
+) where T
+    @boundscheck checkbounds(C, i, j)
+    return (j == length(C.c)) ? (@inbounds -C.c[i]) :
+        i == j+1 ? one(T) : zero(T)
+end
 
 
-function Matrix(C::Companion{T}) where T
-    M = zeros(T, size(C)...)
-    M[:,end]=-C.c
-    for i=1:size(C,1)-1
-        M[i+1, i] = one(T)
-    end
-    M
-end
+# Linear algebra
 
-#Linear algebra stuff
-function mul!(C::Companion{T}, b::Vector{T}) where T
-    x = b[end]
-    y = -C.c[1]*x
-    b[2:end] = b[1:end-1]-C.c[2:end]*x
-    b[1] = y
-    b
+# 3-argument mul! mutates first argument: y <= C * x
+function mul!(y::Vector, C::Companion, x::AbstractVector)
+    @boundscheck length(y) == length(x) == size(C, 1) ||
+        throw(DimensionMismatch("mul! arguments incompatible sizes"))
+    z = x[end]
+    y[1] = -C.c[1] * z
+    y[2:end] = x[begin:end-1] - C.c[2:end] * z
+    return y
 end
-*(C::Companion{T}, b::Vector{T}) where T = mul!(C, copy(b))
 
-function mul!(A::Matrix{T}, C::Companion{T}) where T
-    v = Array{T}(undef, size(A,1))
-    for i=1:size(A,1)
-        v[i] =dot(vec(A[i,:]),-C.c)
+# A <= B * C
+function mul!(A::Matrix, B::AbstractMatrix, C::Companion)
+    @boundscheck (size(A) == (size(B,1), size(C,2)) && size(B, 2) == size(C,1)) ||
+        throw(DimensionMismatch("mul! arguments incompatible sizes"))
+    Base.require_one_based_indexing(B)
+    @views for j in 1:size(A,2)-1
+        @inbounds A[:,j] = B[:,j+1]
     end
-    for i=1:size(A,1), j=1:size(A,2)-1
-        A[i,j] = A[i,j+1]
-    end
-    A[:,end] = v
-    A
+    @inbounds mul!((@view A[:,end]), B, C.c, -1, 0) # A[:,end] <= - B * c
+    return A
 end
-*(A::Matrix{T}, C::Companion{T}) where T = mul!(copy(A), C)
+
 
 function inv(C::Companion{T}) where T
     M = zeros(T, size(C)...)
-    for i=1:size(C,1)-1
+    for i in 1:size(C,1)-1
         M[i, i+1] = one(T)
     end
-    d = M[end, 1] = -one(T)/C.c[1]
-    M[1:end-1, 1] = d*C.c[2:end]
-    M
+    d = M[end, 1] = -one(T) / C.c[1]
+    M[1:end-1, 1] = d * C.c[2:end]
+    return M
 end

test/Project.toml

@@ -0,0 +1,7 @@
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[compat]

test/companion.jl

@@ -1,22 +1,65 @@
-n = rand(1:10)
-Z = Companion(randn(n))
+# test/companion.jl
 
-#Special properties
-@test Matrix(inv(Z)) ≈ inv(Matrix(Z))
+import LinearAlgebra: mul!
+using OffsetArrays: OffsetArray
+import OffsetArrays # for no_offset_view
 
-#Matvec product
-b = randn(n)
-@test Z*b ≈ Matrix(Z)*b
+v = OffsetArray(1:3, -5) # test non-1-based indexing and AbstractVector
+
+Z = @inferred Companion(OffsetArrays.no_offset_view(v))
+
+n = 9
+c = rand(n)
+Z = @inferred Companion(c)
+
+@test Z isa Companion{Float64}
+@test (@inferred size(Z)) == (n,n)
+@test size(Z,1) == n
+@test size(Z,3) == 1
+
+@test Z[2,1] == 1
+@test (@inferred getindex(Z, 2, n)) == -c[2]
+
+@test isassigned(Z, 1)
+@test isassigned(Z, n^2)
+@test !isassigned(Z, 0)
+@test !isassigned(Z, n^2+1)
+
+# Special properties
+Zi = @inferred inv(Z)
+Zm = @inferred Matrix(Z)
+@test Matrix(Zi) ≈ inv(Zm)
+
+# Matvec product
+x = randn(n)
+@test Z * x ≈ Zm * x
+
+y = similar(x)
+@inferred mul!(y, Z, x)
+@test y ≈ Z * x
+
+
+# matrix * companion
+m = 8
+B = randn(m, n)
+@test B * Z ≈ B * Zm
+
+A = copy(B)
+@inferred mul!(A, B, Z)
+@test A ≈ B * Zm
+
+
+# OffsetMatrix * companion
+
+B = OffsetArray(randn(size(B)), 9, -5) # test non-1-based indexing and AbstractMatrix
+@test_throws ArgumentError mul!(A, B, Z)
 
-m = rand(1:10)
-A = randn(m, n)
-@test A*Z ≈ A*Matrix(Z)
 
 # Polynomial construction
 using Polynomials
 p = Polynomial([-1,0,1])
-p_c = Companion(p)
+C = @inferred Companion(p)
 v1 = [1,1]
 v2 = [-1,1]
-@test p_c*v1 == v1
-@test p_c*v2 == -v2
+@test C * v1 == v1
+@test C * v2 == -v2