I need to solve many very large systems of first order ODEs, which describe some chemical reactions. The number of variables (in each system) is on the order of $n \sim 10^5$. I am using ALGLIB vector ODE solver, which under the hood implements Cash-Karp adaptive ODE solver.
Chemical equations are often considered as stiff and it is recommended to use implicit methods. Unfortunately, implicit methods cannot work with such number of variables at they will require inverting huge matrices. This, in turn, requires on the order of $\sim n^3$ operations just to start from and that's before rounding errors are taken into account. Last time when I "played" around with an order smaller matrices ($ n \sim 10^4$), I needed around 50 digits of precision in order to get reasonably correct results. Such precision is not natively supported by current processors and all together this effectively rules out implicit methods for such number of variables.
From another side, adaptive explicit methods, like the one used in ALGLIB, calculate the errors and decrease the step to a very small number. That often makes the solution time insanely large (like months).
My guess is that this happens when some of the variables start to approach zero. As all variables must be strictly non-negative numbers (they are concentrations of some substances), once some of them get close to zero, then the relative errors immediately become very large and the algorithm then decreases the step. In fact, the algorithm would routinely overshoot zero thus making some variables negative. And after that everything would blow up. So, I treat all variables less than zero as exact zeros, when calculating the derivative. That makes the algorithm stable but still does not delete the fact that it likely decreases the step too much.
The parallelization for a single ODE system does not work well, or, shall I say, does not work at all. This is due to the fact that the ODE solver repeatedly calls a fairly quick function, which calculates the derivative. Such functions are based on just pure math and for such functions (even with the input array of $n \sim 10^5$ size) the parallelization overhead just kills the purpose. Rather, the parallelization is achieved by spawning multiple models at the same time where each model runs under a single thread. This can be done spawning multiple threads or multiple external processes (with each thread / process running a single "model"). Both have some benefits and drawbacks and this is irrelevant to the question. The bottom line is that parallelizing a single model kills performance.
The system has an integral of motion: the total number or "atoms" (a sum over all variables multiplied by some weights) in the system must be constant. And since all variables are strictly non-negative, this results in the upper bound for a derivative.
I wonder what can be done to solve such systems of ODEs [much] faster.