The title is probably misleading, but this is a lesson I needed to talk about.
I wrote out some simple code for the quadrature over the reference triangle last time, which involves a double loop. To my chagrin, my immediate reaction to speeding up the code was to put it into Cython, and give it some type declaration.
This did speed up my integrals, but not as much as vectorization. By simply condensing one of the loops into a dot product, and using vector-function evaluation, I sped up my code a substantial amount, especially with higher order integration of “hard” functions.
def quadtriangle_vector(f, w00, w10, x00, x10):
total = 0
for i in range(len(w00)):
total += w00[i] * np.dot(w10 / 2, f([(1 + x00[i]) * (1 - x10) / 2 - 1, x10]))
To see what I mean, consider the following function
from scipy.special import eval_jacobi as jac
return jac(2, 1, 1, np.sin(x - x))
p = 20
x00, w00 = rj(p + 1, 0, 0)
x10, w10 = rj(p + 1, 1, 0)
The speedup I get is staggering.
%timeit quadtriangle(f, w00, w10, x00, x10)
3.11 ms ± 24.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit quadtriangle_vector(f, w00, w10, x00, x10)
348 µs ± 1.67 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Also, I tried to fully vectorize by removing the outer-loop. This actually slowed down the code a bit. Maybe I did it wrong? But for now, I’m decently happy with the speed.