eigsh smallest eigenvalues: sigma=0 shift-invert is 15-27x faster than which='SM'
SciPy 1.18.1 / NumPy 2.5.4, 2-D Dirichlet Laplacian (5-point), k=6 smallest, checked against the analytic eigenvalues. which='SM' and which='SA' both converged (no ArpackNoConvergence) but slowly: 0.36 s / 6.9 s / 51 s at n=1e4 / 4e4 / 9e4, with rel. error ~1e-13. eigsh(A, k=6, sigma=0) (shift-invert, one sparse LU) took 0.08 / 0.49 / 1.9 s with rel. error ~1e-14. Times are from a shared 8-core box under load, so treat the ratios as rough. The usual worry that sigma=0 fails on a singular matrix did not hold here: for a Neumann Laplacian (exact zero eigenvalue, n=1e4) sigma=0 returned all 6 eigenvalues to 5e-17 with no warning, because rounding gives the LU a tiny nonzero pivot. That is luck, not a guarantee, so a slightly negative shift (sigma=-1e-3) is the safer default and was just as accurate. The cost of shift-invert is the LU fill-in, so for 3-D grids check memory first.