Linear system solution (LU)
I did some digging around and it seems LU decomposition was recently added to Arraymancer (#382).

*Somewhat sidetracking: any reason why the LAPACK LU result is duplicated in two tensors in `lu_permuted`? It seems like an easy memory optimization to avoid this duplication  (each matrix is half zeros). Scipy's wrapper to `DGETRF`, `scipy.linalg.lu_factor`, returns a single array lu as the wrapped function does.*

Ok, back to the main focus: solving a linear system of equations `AX = B` is a very common thing to do, and LU factorization is a classic, widely applicable and decently efficient way to do it. We could wrap `DGETRS`, except:

1. `DGETRS` wants to do the permutation which is already done in arraymancer's `lu_permuted` wrapper.

2. `DGETRS` expects a single LU matrix as returned by `DGETRF`.

Proposals:

1. Merge `L` and `U` back into a single matrix (`L + U`) and pass a no-op array `IPIV` to `DGETRS`. Easiest but least efficient option.

2. Directly call `DTRSM` from LAPACK twice, which is the routine that actually does the forward-backward subs.

3. Write a pure-nim FBS proc(s) equivalent to `DTRSM`. A naive implementation seems rather straightforward, but I'm not sure if `DTRSM` has additional optimization that would involve more complex implementation and have significant benefits. Some investigation needed here.

4. Don't use `lu_permuted` at all, and wrap `DGESV` directly.

Comments?