I need to solve the following problem: $$ \begin{cases} \dot{\vec{x(t)}} = A\vec{x(t)} + u(t)D\vec{x(t)} + u(t)\vec{b}, & x \in (0, T), \\ \vec{x(0)} = \vec{0}, \end{cases}$$ where $u(t)$ is known function, $A$ and $D$ are known constant matrices and $b$ is known constant vector.

I solve this problem with explicit Euler and accuracy of the grid of $10^5$ nodes is quite enough for me. However, I need to solve such systems a lot of times, so my computations takes about a week. I am looking for a way of parallelizing numerical integration of this ODE. Can you advice some papers concerning this topic?

I have found the approach called Parareal, but I am not sure that it is the most efficient way.

Great thanks for any help, advice or papers!

  • 1
    $\begingroup$ You could take a look at "Bifurcation & chaos in nonlinear structural dynamics: Novel & highly efficient optimal-feedback accelerated Picard iteration algorithms" by Wang et al. sciencedirect.com/science/article/pii/S1007570418301515 . Not sure how easy it is to parallelize some of the operations suggested in that paper. $\endgroup$ Commented May 24, 2018 at 20:45
  • $\begingroup$ Do you want to parallelize your integrator? Or, do you want run different cases in parallel? $\endgroup$
    – nicoguaro
    Commented May 25, 2018 at 13:04
  • $\begingroup$ @nicoguaro, I want to parallelize my integrator. $\endgroup$ Commented May 27, 2018 at 20:02

1 Answer 1


First of all, don't even consider "optimizing" before you're using the right integration method. Chucking more computers at the problem may sound like the easiest way to solve it, but in reality it will be more difficult than you'd think (due to implicit BLAS multithreading explained later). If you're using explicit Euler, you can do a whole lot better just by changing the method. Proper adaptive time stepping (possibly implicit depending on the stiffness) algorithms can be orders of magnitude more efficient than explicit Euler. Don't even consider parallelizing and keeping explicit Euler unless you have a very good reason (i.e. for hyperbolic PDEs so you're concerned about SSP, but in that case just use an SSP-optimized integrator... you're not going to have a good reason).

Next, your derivative calculation function is by far the most expensive part of the calculation. First of all, you should make sure that it's optimized. Language won't matter much here since all of the time will be in the BLAS kernel for matrix multiplication, but make sure it's non-allocating via in-place operations, using BLAS or SIMD'd loops as much as possible, etc.

Once that's all good, consider whether you still need to parallelize. With the right ODE methods and with a decently performant derivative implementation, 10^5 can be reasonably speedy. If you ended up using an implicit method, you'll want to do things like swap out the linear solver to some iterative solver, and that will give you another nice speedup.

But since this is a post about parallelism, I'll keep going, assuming you want to do something extreme.

How to parallelize this: the details to know

Since $A$ and $D$ are matrices, if you're using a language with a standard BLAS setup then you're already parallelizing via multithreading. Just check your CPU usage (ex: htop on Linux) and you'll see that large enough matrix multiplications with multithread in Python, R, Julia, MATLAB, etc. So, getting more performance out isn't as simple as just parallelizing it since if you're on a single computer it's likely you're using all of your cores optimally already.

Parallel in time algorithms like parareal only make sense when the ODE is cheap since there's no other available parallelism. Parallel in time algorithms are not efficient, but they try to make up for it by allowing parallelism. The overhead cutoff can be quite high for these algorithms, I've read that 32 cores is a good low estimate. These are definitely not suitable to the way you've described your problem.

However, there's two other levels of parallelism. If your solving is expensive, then solving at multiple systems at the same time in parallel will be almost optimal speedup. i.e. do a parallel for loop or map that solves 16 systems at a time in a 16 core machine and it will be about 16x speedup since the data transfer times are almost 0 compared to the full integration. For example, if you're using Julia, you can write this up yourself quite easily, or use what's built into DifferentialEquations.jl, or find the right resource for whatever other language. Remember, BLAS is already multithreading (and most likely those matrix multiplications are the most expensive operation in the loop) so this only makes sense if you can farm out to other computers/nodes. Julia's native parallelism does that, Python's multiprocessing can do it, you can use MPI in C/Fortran, etc. If you have a big cluster and want lots of gains for little work, this is a way to solve your problem.

Finally, if you have an optimized derivative function with most of the time in matrix operations but don't have extra nodes hanging around (or if you do, you can still apply this), another level of parallelism you can exploit is within-method. You're already doing within-level parallelism via the multithreaded BLAS operations, but you can do better. GPUs are essentially parallel linear algebra machines, so if you have one hanging around then you can just perform the matrix multiplications on then it's a good idea. Setting it up isn't hard since things like CuBLAS are just BLAS on GPUs, so you just have to call it. However, I think you'll have memory issues if you do it naively, so this could be difficult.


Stop using explicit Euler: the first thing you should do is really optimize your code and your method. If that's not enough, solve the system on different nodes at the same time. If that's not enough, accelerate the linear algebra with GPUs/Xeon Phis/etc. since the limiting step is the matrix operations.

For reference, here's a notebook explaining the optimizing step and here's a blog post explaining within-method GPU parallelism.

  • $\begingroup$ Thank you very much for the detailed answer! Actually, I did not say my real problem. This dynamic system is just a model from which my problem derives. In my real problem I have a system of discrete equations which looks just like explicit Euler for dynamic system. Then I have one additional equation like $y_i = c^T x_i$. Given $10^5$ values $u_i$ and $\widehat y_i$ and I need to find the values of $A$, $D$, $b$ and $c$ which give the values of $y_i$ the most close to $\widehat y_i$. So I put unknown values to an optimizer and compute the values of $x_i$ a lot of times. $\endgroup$ Commented May 26, 2018 at 16:58
  • $\begingroup$ So the matrices in my experiments have the dimension about $10$ and parallel matrix multiplication does not help $\endgroup$ Commented May 26, 2018 at 16:59
  • $\begingroup$ @StanislavMorozov Wouldn't it be a lot faster to solve that using something like least squares. Namely you can always do a similarity transformation of $x$ such that it is a vector of the past $n$ outputs, such that the only unknowns in the difference equation are the constant matrices. $\endgroup$
    – fibonatic
    Commented Jun 2, 2018 at 11:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.