Slow when BDF is used

[Nicholaswogan]

Great package, thanks! I've noticed that when the BDF method is used, the code is slow compared to the original fortran. This might have to do with the use of std::vector<std::vector<double>> to hold the jacobian? It should be faster to use a old style multidimensional C array, or a std::vector<double> with the proper indexing.

EDIT:
I compared to https://github.com/sdwfrost/liblsoda , and the speed was the same. So there isn't a problem with std::vector .

For context, I made a python wrapper to the c++: https://github.com/Nicholaswogan/NumbaLSODA . I was benchmarking it against the scipy wrapper of LSODA, and noticed that when the BDF method is used, my python wrapper was much slower than scipy. This should not be the case because my wrapper should never invoke the python interpreter, while scipy's does.

[Nicholaswogan]

The reason scipy's LSODA implementation is faster is because it uses faster BLAS (openblas or intel mkl). The speed differences between Scipy's implementation and your C++ implementation only become apparent for big problems (>500 odes) with big matrix operations, where it actually maters that you are using fast BLAS.

[dilawar]

Thanks for notes @Nicholaswogan . This fell through the cracks at my end.

I ported it to C++ for the MOOSE simulator which rarely uses this method: the results differ wildly vis-a-vis default RK5 (chemical ODE systems). I am still not sure if the issue is with LSODA or due to some bug in the port. At least the number computed by this port and the scipy gives roughly the same numbers.

Feel free to send a PR with files that you used to benchmark this code.