How can I increase the solution speed of system that looks like this:

import numpy as np
from scipy.sparse import diags
import scipy.sparse.linalg as spla

n = 1_000_000  
diagonals = [-0.01*np.ones(n-1), np.ones(n), -0.01*np.ones(n-1)]
offsets = [-1, 0, 1]  

mat = diags(diagonals, offsets, shape=(n, n), format='csc')

LU_mat = spla.splu(mat)  

b = np.random.randn(n)

for i in range(654):
    b = LU_mat.solve(b) 

On my machine this example appears to be single threaded.

How can I efficiently use the resources of my machine to solve this system?
Is there a threaded sparse LU solver for python?