Speed up the legendre root finding (used in gauss-legendre grid) #629
Conversation
… harmonics and integration
src/utils/legendre_roots.hpp
Outdated
|
|
||
| // Parameters for Newton iteration | ||
| double tol = 1e-14; | ||
| int max_iter = 100000; |
There was a problem hiding this comment.
I think 100,000 iterations is way too many. You should be gaining 1 bit of precision at least per iteration and there are only 52 bits in a double precision number.
There was a problem hiding this comment.
Done! I left it the same as the N for bisection and forgot to update. The typical number of iteration to reach machine precision is 3, but I have set it to 20 for safety.
src/utils/legendre_roots.hpp
Outdated
| // For a large n, you might not need extremely large N | ||
| // if the function is well-behaved. | ||
| // But here we keep your original 60000 intervals. | ||
| roots_weights[1][i] = Ci(i, n, xi, -1.0, 1.0, 60000); |
There was a problem hiding this comment.
Use analytic expression for the roots weights, e.g., Canuto, Hussaini, Quarteroni, Zang, Eq. (2.3.10)
There was a problem hiding this comment.
Done. Thanks for the reference!
src/utils/legendre_roots.hpp
Outdated
| std::vector<double> xi(n); | ||
|
|
||
| // Parameters for Newton iteration | ||
| double tol = 1e-14; |
There was a problem hiding this comment.
I would iterate until floating point precision here. In my experience when
There was a problem hiding this comment.
Now set to std::numeric_limits::epsilon()
|
@jmstone Ready to be merged! |
updating the bisection solver by a newton raphson finder as @dradice suggested. Now speed up the code by quite a bit. I also included a unit test for doing cross integration of spherical harmonics on Gauss-Legendre grid. The integral is calculated robustly up to l_max=20, loss of precision seen at l_max>=25, generating an error on the order of 1e-8 for higher l's. This is quite enough for CCE but should be kept in mind for future usage. Perhaps using precomputed values for the legendre roots and weights, or using the gsl implementation would be more robust.