I’m developing a package for physics simulations, where the core routine solves a dense linear system of complex equations,
\mathbf{A} \mathbf{X} = \mathbf{B}
The matrix \mathbf{A} can be small in some applications, for example 6\times6, but other use cases at the other extreme would be more like 10,000\times10,000. I’m wondering if it makes sense to have a switch in the high-level (user-facing) routine depending on the size of the problem,
function simulation(N)
# using random numbers as dummy calculation
if N > 100
A = rand(Complex{Float64}, (N,N))
Y = rand(Complex{Float64}, (N, 5))
# A = Matrix{Complex{Float64}}(undef, N, N)
# Y = Matrix{Complex{Float64}}(undef, N, 5)
@info "using regular Array"
elseif N <= 100
A = rand(SMatrix{N,N,Complex{Float64},N^2})
Y = rand(SMatrix{N,5,Complex{Float64},5N})
@info "using static Array"
end
X = A \ Y
# [...] more calculations with X and Y in the real-world situation
return abs2.(X)
end
simulation(3)
simulation(300)
Does this make sense? Is there something already in StaticArrays to take care of this sort of situation?
this naive implementation will have type instability, so be sure to do it in a separate function. but theoretically it makes sense, if you use the small matrix enough times it can be worth while
Can you please clarify what part(s) I should do in a separate function? At some point there will have to be a if(N<100) returning either Array or SArray, depending on N. The matrix is typically used a hundred times, so it seems worth contemplating, especially when size ~ 3x3 etc. where solve() has special cases.
This doesn’t seem worth optimizing for small matrices — you can solve 100 6x6 systems in a few microseconds, so why do you care how fast it is? If you are solving only 100 times, only the bigger matrix sizes will matter.
The matrix is used a few hundred times in this function, but the function itself might be called 1000s of times by the user. Granted it’s still not that long to run, but I’d like to see if there’s a good idiomatic way to deal with this kind of situation, if just for my own learning.
Basically, there are two cases: (1) you’re calling it enough times that the costs of compiling a specialized version for a specific matrix size (a StaticArray) and dispatching to it dynamically are worth it, or (2) you aren’t. Since even calling it 100000 times probably takes less than a second for a 6x6 matrix, you may be in case (2).