Clean up linear algebra (related to modes)

I initially planned to implement a computation of modes with the sparse unit cell here. This, however, prompted code review and subsequent cleanup.

In summary, this MR:

• Eliminates our custom wrappers for LAPACK lu and gen_eig. The former I reimplemented using scipy.linalg.lapack, and the latter wasn't used anywhere.
• Shortens linalg tests by using pytest parameterization.
• Uses a @ b instead of np.dot(a, b) or a.dot(b). This helps both readability and customization, because @ is properly dispatched by sparse matrices (unlike .dot).
• Requires that setup_linsys takes square matrices as input. For some time already, we only passed square matrices to it, see this comment.