Implementation of the Tri-Diagonal Matrix Algorithm (TDMA) with NumPy

Question

I am using TDMA in Python using NumPy. The tridiagonal matrix is ​​stored in three arrays:

a = array([...]) b = array([...]) c = array([...]) 

I would like to efficiently calculate alpha coefficients. The algorithm is as follows:

 # n = size of the given matrix - 1 alpha = zeros(n) alpha[0] = b[0] / c[0] for i in range(n-1): alpha[i+1] = b[i] / (c[i] - a[i] * alpha[i]) 

However, this is inefficient due to the Python for loop. I want me to want something like this approach:

 # n = size of the given matrix - 1 alpha = zeros(n) alpha[0] = b[0] / c[0] alpha[1:] = b[1:] / (c[1:] - a[1:] * alpha[:-1]) 

In this latter case, the result is incorrect, because NumPy stores the right side of the last expression in a temporary array and then assigns links to its elements on alpha[1:] . Therefore, a[1:] * alpha[:-1] is just an array of zeros.

Can I tell NumPy to use the alpha values ​​calculated in the previous steps inside my inner loop?

Thanks.

Answer 1

There seems to be no way in Python to do this without using C or its python variations.

Answer 2

If its tridiagonal systems that you want to solve is solve_banded() in numpy.linalg . Not sure if this is what you are looking for.