Skip to content

Commit d60e8f8

Browse files
authored
reduce allocations and improve performance in sminreal (#652)
* reduce allocations and improve performance in sminreal * add pointers between minreal and sminreal docstrings
1 parent 4012c10 commit d60e8f8

File tree

2 files changed

+15
-6
lines changed

2 files changed

+15
-6
lines changed

src/simplification.jl

Lines changed: 13 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -13,10 +13,11 @@ trunc_zero!(A) = A[abs.(A) .< 10eps(maximum(abs, A))] .= 0
1313
trunc_zero!(sys.A); trunc_zero!(sys.B); trunc_zero!(sys.C)
1414
sminreal(sys)
1515
```
16+
See also [`minreal`](@ref)
1617
"""
1718
function sminreal(sys::StateSpace)
1819
A, B, C, inds = struct_ctrb_obsv(sys)
19-
return StateSpace(A, B, C, sys.D, sys.timeevol)
20+
return basetype(sys)(A, B, C, sys.D, sys.timeevol)
2021
end
2122

2223
# Determine the structurally controllable and observable realization for the system
@@ -39,10 +40,16 @@ Compute a bit-vector, expressing whether a state of the pair (A, B) is
3940
structurally controllable based on the location of zeros in the matrices.
4041
"""
4142
function struct_ctrb_states(A::AbstractVecOrMat, B::AbstractVecOrMat)
42-
bitA = .!iszero.(A)
43+
size(A,1) > typemax(UInt16) && error("Maximum size of A exceeded. If you encounter this error, please open an issue. This limit is not fundamental and excists for performance reasons only.")
44+
# UInt16 can only store up to 65535, so if A is completely dense and of size larger than 65535, the computations below might overflow. This is exceedingly unlikely though.
45+
bitA = UInt16.(.!iszero.(A)) # Convert to Int because mutiplying with a bit matrix is slow
4346
x = vec(any(B .!= 0, dims=2)) # index vector indicating states that have been affected by input
44-
for i = 1:size(A, 1) # apply A nx times, similar to controllability matrix
45-
x = x .| .!iszero.(bitA * x)
47+
xi = bitA * x
48+
xi2 = similar(xi)
49+
@. xi = (xi != false) | !iszero(x)
50+
for i = 2:size(A, 1) # apply A nx times, similar to controllability matrix
51+
mul!(xi2, bitA, xi)
52+
@. xi = (xi2 != false) | !iszero(xi)
4653
end
47-
x
48-
end
54+
xi .!= false # Convert back to BitVector
55+
end

src/types/StateSpace.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -389,6 +389,8 @@ end
389389
Minimal realisation algorithm from P. Van Dooreen, The generalized eigenstructure problem in linear system theory, IEEE Transactions on Automatic Control
390390
391391
For information about the options, see `?ControlSystems.MatrixPencils.lsminreal`
392+
393+
See also [`sminreal`](@ref), which is both numerically exact and substantially faster than `minreal`, but with a much more limited potential in removing non-minimal dynamics.
392394
"""
393395
function minreal(sys::T, tol=nothing; fast=false, atol=0.0, kwargs...) where T <: AbstractStateSpace
394396
A,B,C,D = ssdata(sys)

0 commit comments

Comments
 (0)