-
-
Save TheoChristiaanse/d168b7e57dd30342a81aa1dc4eb3e469 to your computer and use it in GitHub Desktop.
Tridiagonal Matrix Algorithm solver in Python. I've modified the code from cbellei so, it works with python 3.0+ and implemented the use of jit to increase the speed.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import numpy as np | |
from numba import jit, f8 | |
## Tri Diagonal Matrix Algorithm(a.k.a Thomas algorithm) solver | |
@jit(f8[:] (f8[:],f8[:],f8[:],f8[:] )) | |
def TDMAsolver(a, b, c, d): | |
''' | |
TDMA solver, a b c d can be NumPy array type or Python list type. | |
refer to http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm | |
and to http://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm) | |
''' | |
nf = len(d) # number of equations | |
ac, bc, cc, dc = map(np.array, (a, b, c, d)) # copy arrays | |
for it in range(1, nf): | |
mc = ac[it-1]/bc[it-1] | |
bc[it] = bc[it] - mc*cc[it-1] | |
dc[it] = dc[it] - mc*dc[it-1] | |
xc = bc | |
xc[-1] = dc[-1]/bc[-1] | |
for il in range(nf-2, -1, -1): | |
xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il] | |
return xc |
I run your code, and got different speed results.
Test results:
[ 0.14877589 0.75612053 -1.00188324 2.25141243]
[ 0.14877589 0.75612053 -1.00188324 2.25141243]
Speed results:
jit_new 6.220152854919434
control 0.7551865577697754I am not familiar with the numba or jit. So I don't know if it is caused by this. Anyway, thank you for your improvement.
Indeed the code gives wrong answers. I tried to compare the results also.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
@ahwillia The file is forked from two other works that did not have any license file. I think you need to go to the original work ask him to license it and then I can add a license based on his license. Currently, under Github policy I can fork code and make it my own thing but since it is a derivative of some original code I think that is way it needs to go. But I'm no lawyer so, if you can advise me on this would be awesome. If I can I would release it under the most open license I can find.