solve_banded is very slow as a tridiagonal solver
Hello, I noticed that using solve_banded to solve a tridiagonal linear system is at least an order of magnitude slower than it should be based on a straightforward operation count of the Thomas algorithm. In the python source code, solve_banded is pretty much a straightforward call into gbsv from lapack, but I did not dare to follow up on what is going on in there. Maybe it's doing pivoting, or there is some inefficiency due to being coded for more general band-diagonal problems. In any event, in a lot of textbook applications of tridiagonal solvers, pivoting is not necessary (see, e.g., Higham, Accuracy and Stability of Numerical Algorithms, SIAM 1996, Section 9.5) and the Thomas algorithm does not vectorize efficiently. So it would be very good if either solve_banded could be fixed if it does something stupid, or otherwise if a plain Thomas algorithm tridiagonal solver could be added as this would cover a lot of run-off-the-mill problems while solve_banded could handle the rest. In any event, just to demonstrate that the performance gap is really there, I coded up a Thomas algorithm in cython and here is the output for a tridiagonal system of size 10^6: CPU time Thomas algorithm in python: 3.47022 CPU time solve_banded: 0.428928 CPU time Thomas algorithm in cython: 0.032048 Speedup python/cython: 108.281952072 Speedup solve_banded/cython: 13.3839241138 The code is attached, but it's quickly put together and I am not claiming that it's well. In fact, it's the first time I try cython which was easy enough, except that @cython.boundscheck(False) makes the code run slower, not faster... Regards, Marcel
First I would try updating to version 0.16. Solve_banded was updated to use dgtsv when given a tridiagonal system. Dgtsv does do partial pivoting, but it should be much faster than dgbsv for this case. For a system of 1e6, it is about 6 times faster for me ( 39.1 ms vs 270 ms). Eric On Fri, Aug 14, 2015 at 12:45 PM, Marcel Oliver < m.oliver@jacobs-university.de> wrote:
Hello,
I noticed that using solve_banded to solve a tridiagonal linear system is at least an order of magnitude slower than it should be based on a straightforward operation count of the Thomas algorithm.
In the python source code, solve_banded is pretty much a straightforward call into gbsv from lapack, but I did not dare to follow up on what is going on in there. Maybe it's doing pivoting, or there is some inefficiency due to being coded for more general band-diagonal problems.
In any event, in a lot of textbook applications of tridiagonal solvers, pivoting is not necessary (see, e.g., Higham, Accuracy and Stability of Numerical Algorithms, SIAM 1996, Section 9.5) and the Thomas algorithm does not vectorize efficiently.
So it would be very good if either solve_banded could be fixed if it does something stupid, or otherwise if a plain Thomas algorithm tridiagonal solver could be added as this would cover a lot of run-off-the-mill problems while solve_banded could handle the rest.
In any event, just to demonstrate that the performance gap is really there, I coded up a Thomas algorithm in cython and here is the output for a tridiagonal system of size 10^6:
CPU time Thomas algorithm in python: 3.47022 CPU time solve_banded: 0.428928 CPU time Thomas algorithm in cython: 0.032048
Speedup python/cython: 108.281952072 Speedup solve_banded/cython: 13.3839241138
The code is attached, but it's quickly put together and I am not claiming that it's well. In fact, it's the first time I try cython which was easy enough, except that @cython.boundscheck(False) makes the code run slower, not faster...
Regards, Marcel
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
Eric Moore writes:
First I would try updating to version 0.16. Solve_banded was updated to use dgtsv when given a tridiagonal system. Dgtsv does do partial pivoting, but it should be much faster than dgbsv for this case. For a system of 1e6, it is about 6 times faster for me ( 39.1 ms vs 270 ms).
Eric
Dear Eric, thank you very much for the quick reply. I'll follow this issue and, if necessary, report back. For now, we'll probably just go with my cython hack as it's unreasonable to expect my collaborators to hand-update scipy. But the tridiagonal issue comes up often enough that it's worth some thought and optimization. Even if it's only a factor two in the end (I'll check once I get to update), a plain Thomas tridiagonal solver would probably be worth it. With best regards, Marcel
participants (2)
-
Eric Moore -
Marcel Oliver