Something Haskell has lacked for a long time is a good medium-duty
linear system solver based on the LU decomposition.  There are bindings
to the usual C/Fortran libraries, but not one in pure Haskell.  (An
example "LU factorization" routine that does not do partial pivoting has
been around for years, but lacking pivoting it can fail unexpectedly on
well-conditioned inputs.  Another Haskell LU decomposition using partial
pivoting is around, but it uses an inefficient representation of the
pivot matrix, so it's not suited to solving systems of more than 100 x
100, say.)

By medium duty I mean a linear system solver that can handle systems of
(1000s) x (1000s) and uses Crout's efficient in-place algorithm.  In
short, a program does everything short of exploiting SIMD vector
instructions for solving small subproblems.

Instead of complaining about this, I have written a little library that
implements this.  It contains an LU factorization function and an LU
system solver.  The LU factorization also returns the parity of the
pivots ( = (-1)^(number of row swaps) ) so it can be used to calculate
determinants.  I used Gustavson's recursive (imperative) version of
Crout's method.  The implementation is quite simple and deserves to be
better known by people using functional languages to do numeric work. 
The library can be downloaded from GitHub:
https://github.com/gwright83/luSolve