luSolve fails for matrices with zeros on the main diagonal

The implementation of `luSolve` in `Linear.Matrix` is sensitive to the ordering of equations. In particular, it fails in cases where the coefficient matrix has at least one zero on the main diagonal, even if the rows could be permuted so that there are no zeros on the main diagonal.

Consider the following code:

```haskell
{-# LANGUAGE DataKinds #-}

import Data.Vector qualified as V
import Linear.Matrix (luSolve)
import Linear.V

vector = V . V.fromList

matrix = V . V.fromList . map (V . V.fromList)

ex1 =
let a :: V 2 (V 2 Double)
a = matrix [[1, 2], [0, 4]]
b = vector [5, 6]
in luSolve a b

ex2 =
let a :: V 2 (V 2 Double)
a = matrix [[0, 4], [1, 2]]
b = vector [6, 5]
in luSolve a b
```

`ex1` and `ex2` represent the same system of equations, so I'd expect them to be equivalent. But

```haskell
>>> ex1
V {toVector = [2.0,1.5]}
>>> ex2
V {toVector = [NaN,NaN]}
```

Implementations of `luSolve` in other packages (e.g. hmatrix, scipy) are robust to ordering; for example, in hmatrix:

```haskell
import Numeric.LinearAlgebra

ex1' =
let a :: Matrix Double
a = (2 >< 2) [1, 2, 0, 4]
b = (2 >< 1) [5, 6]
lu = luPacked a
in luSolve lu b

ex2' =
let a :: Matrix Double
a = (2 >< 2) [0, 4, 1, 2]
b = (2 >< 1) [6, 5]
lu = luPacked a
in luSolve lu b
```

```haskell
>>> ex1'
(2><1)
[ 2.0
, 1.5 ]
>>> ex2'
(2><1)
[ 2.0
, 1.5 ]
```