Skip to content

Commit dbe1787

Browse files
stevengjstevengj
committed
fix load-order issues by moving evalpoly methods to a separate file
1 parent d5be83b commit dbe1787

File tree

3 files changed

+83
-80
lines changed

3 files changed

+83
-80
lines changed

‎src/LinearAlgebra.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -561,6 +561,8 @@ include("schur.jl")
561561
include("structuredbroadcast.jl")
562562
include("deprecated.jl")
563563

564+
include("evalpoly.jl")
565+
564566
const ⋅ = dot
565567
const × = cross
566568
export ⋅, ×

‎src/evalpoly.jl‎

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
# matrix methods for evalpoly(X, p) = ∑ₖ Xᵏ⁻¹ p[k]
2+
3+
# non-inplace fallback for evalpoly(X, p)
4+
function _evalpoly(X::AbstractMatrix, p)
5+
Base.require_one_based_indexing(p)
6+
p0 = isempty(p) ? Base.reduce_empty_iter(+, p) : p[end]
7+
Xone = one(X)
8+
S = Base.promote_op(*, typeof(Xone), typeof(Xone))(Xone) * p0
9+
for i = length(p)-1:-1:1
10+
S = X * S + @inbounds(p[i] isa AbstractMatrix ? p[i] : p[i] * I)
11+
end
12+
return S
13+
end
14+
15+
_scalarval(x::Number) = x
16+
_scalarval(x::UniformScaling) = x.λ
17+
18+
"""
19+
evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p)
20+
21+
Evaluate the matrix polynomial ``Y = \\sum_k X^{k-1} p[k]``, storing the result
22+
in-place in `Y`, for the coefficients `p[k]` (a vector or tuple). The coefficients
23+
can be scalars, matrices, or [`UniformScaling`](@ref).
24+
25+
Similar to `evalpoly`, but may be more efficient by working more in-place. (Some
26+
allocations may still be required, however.)
27+
"""
28+
function evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p::Union{AbstractVector,Tuple})
29+
@boundscheck axes(Y,1) == axes(Y,2) == axes(X,1) == axes(X,2)
30+
Base.require_one_based_indexing(p)
31+
32+
N = length(p)
33+
pN = iszero(N) ? Base.reduce_empty_iter(+, p) : p[N]
34+
if pN isa AbstractMatrix
35+
Y .= pN
36+
elseif N > 1 && p[N-1] isa Union{Number,UniformScaling}
37+
# initialize Y to p[N-1] I + X p[N], in-place
38+
Y .= X .* _scalarval(pN)
39+
for i in axes(Y,1)
40+
@inbounds Y[i,i] += p[N-1] * I
41+
end
42+
N -= 1
43+
else
44+
# initialize Y to one(Y) * pN in-place
45+
for i in axes(Y,1)
46+
for j in axes(Y,2)
47+
@inbounds Y[i,j] = zero(Y[i,j])
48+
end
49+
@inbounds Y[i,i] += one(Y[i,i]) * pN
50+
end
51+
end
52+
if N > 1
53+
Z = similar(Y) # workspace for mul!
54+
for i = N-1:-1:1
55+
mul!(Z, X, Y)
56+
if p[i] isa AbstractMatrix
57+
Y .= p[i] .+ Z
58+
else
59+
# Y = p[i] * I + Z, in-place
60+
Y .= Z
61+
for j in axes(Y,1)
62+
@inbounds Y[j,j] += p[i] * I
63+
end
64+
end
65+
end
66+
end
67+
return Y
68+
end
69+
70+
Base.evalpoly(X::AbstractMatrix, p::Tuple) = _evalpoly(X, p)
71+
Base.evalpoly(X::AbstractMatrix, ::Tuple{}) = zero(one(X)) # dimensionless zero, i.e. 0 * X^0
72+
Base.evalpoly(X::AbstractMatrix, p::AbstractVector) = _evalpoly(X, p)
73+
74+
Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{Union{Number, UniformScaling}, Vararg{Union{Number, UniformScaling}}}) =
75+
evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p)
76+
Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{AbstractMatrix{<:Number}, Vararg{AbstractMatrix{<:Number}}}) =
77+
evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p)
78+
Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:Union{Number, UniformScaling}}) =
79+
isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p)
80+
Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:AbstractMatrix{<:Number}}) =
81+
isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p)

‎src/generic.jl‎

Lines changed: 0 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -2091,83 +2091,3 @@ end
20912091
function copytrito!(B::StridedMatrixStride1{T}, A::StridedMatrixStride1{T}, uplo::AbstractChar) where {T<:BlasFloat}
20922092
LAPACK.lacpy!(B, A, uplo)
20932093
end
2094-
2095-
# non-inplace fallback for evalpoly(X, p)
2096-
function _evalpoly(X::AbstractMatrix, p)
2097-
Base.require_one_based_indexing(p)
2098-
p0 = isempty(p) ? Base.reduce_empty_iter(+, p) : p[end]
2099-
Xone = one(X)
2100-
S = Base.promote_op(*, typeof(Xone), typeof(Xone))(Xone) * p0
2101-
for i = length(p)-1:-1:1
2102-
S = X * S + @inbounds(p[i] isa AbstractMatrix ? p[i] : p[i] * I)
2103-
end
2104-
return S
2105-
end
2106-
2107-
_scalarval(x::Number) = x
2108-
_scalarval(x::UniformScaling) = x.λ
2109-
2110-
"""
2111-
evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p)
2112-
2113-
Evaluate the matrix polynomial ``Y = \\sum_k X^{k-1} p[k]``, storing the result
2114-
in-place in `Y`, for the coefficients `p[k]` (a vector or tuple). The coefficients
2115-
can be scalars, matrices, or [`UniformScaling`](@ref).
2116-
2117-
Similar to `evalpoly`, but may be more efficient by working more in-place. (Some
2118-
allocations may still be required, however.)
2119-
"""
2120-
function evalpoly!(Y::AbstractMatrix, X::AbstractMatrix, p::Union{AbstractVector,Tuple})
2121-
@boundscheck axes(Y,1) == axes(Y,2) == axes(X,1) == axes(X,2)
2122-
Base.require_one_based_indexing(p)
2123-
2124-
N = length(p)
2125-
pN = iszero(N) ? Base.reduce_empty_iter(+, p) : p[N]
2126-
if pN isa AbstractMatrix
2127-
Y .= pN
2128-
elseif N > 1 && p[N-1] isa Union{Number,UniformScaling}
2129-
# initialize Y to p[N-1] I + X p[N], in-place
2130-
Y .= X .* _scalarval(pN)
2131-
for i in axes(Y,1)
2132-
@inbounds Y[i,i] += p[N-1] * I
2133-
end
2134-
N -= 1
2135-
else
2136-
# initialize Y to one(Y) * pN in-place
2137-
for i in axes(Y,1)
2138-
for j in axes(Y,2)
2139-
@inbounds Y[i,j] = zero(Y[i,j])
2140-
end
2141-
@inbounds Y[i,i] += one(Y[i,i]) * pN
2142-
end
2143-
end
2144-
if N > 1
2145-
Z = similar(Y) # workspace for mul!
2146-
for i = N-1:-1:1
2147-
mul!(Z, X, Y)
2148-
if p[i] isa AbstractMatrix
2149-
Y .= p[i] .+ Z
2150-
else
2151-
# Y = p[i] * I + Z, in-place
2152-
Y .= Z
2153-
for j in axes(Y,1)
2154-
@inbounds Y[j,j] += p[i] * I
2155-
end
2156-
end
2157-
end
2158-
end
2159-
return Y
2160-
end
2161-
2162-
Base.evalpoly(X::AbstractMatrix, p::Tuple) = _evalpoly(X, p)
2163-
Base.evalpoly(X::AbstractMatrix, ::Tuple{}) = zero(one(X)) # dimensionless zero, i.e. 0 * X^0
2164-
Base.evalpoly(X::AbstractMatrix, p::AbstractVector) = _evalpoly(X, p)
2165-
2166-
Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{Union{Number, UniformScaling}, Vararg{Union{Number, UniformScaling}}}) =
2167-
evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p)
2168-
Base.evalpoly(X::StridedMatrix{<:Number}, p::Tuple{AbstractMatrix{<:Number}, Vararg{AbstractMatrix{<:Number}}}) =
2169-
evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p)
2170-
Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:Union{Number, UniformScaling}}) =
2171-
isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), typeof(_scalarval(p[begin])))), X, p)
2172-
Base.evalpoly(X::StridedMatrix{<:Number}, p::AbstractVector{<:AbstractMatrix{<:Number}}) =
2173-
isempty(p) ? _evalpoly(X, p) : evalpoly!(similar(X, Base.promote_op(*, typeof(one(eltype(X))), eltype(p[begin]))), X, p)

0 commit comments

Comments
 (0)