diff --git a/Project.toml b/Project.toml index e7be8b0..c565f47 100644 --- a/Project.toml +++ b/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"] diff --git a/docs/Project.toml b/docs/Project.toml index 9da1eee..fcaf036 100644 --- a/docs/Project.toml +++ b/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" diff --git a/src/SpecialMatrices.jl b/src/SpecialMatrices.jl index e4ce29c..b10c630 100644 --- a/src/SpecialMatrices.jl +++ b/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 diff --git a/src/companion.jl b/src/companion.jl index 646b006..3e81822 100644 --- a/src/companion.jl +++ b/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 diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 0000000..adbac9d --- /dev/null +++ b/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] diff --git a/test/companion.jl b/test/companion.jl index bee2a83..c6c6f1a 100644 --- a/test/companion.jl +++ b/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