[Propose] Add new method for LU decomposition

[yoshoku]

I would like to add new method for LU decomposition like scipy.linalg.lu. The already existing lu_fact method returns the factor and permutation matrices in a form different from the form of LU decomposition. In contrast, the lu method in this PR returns three matrices according to the form of LU decomposition that is A = P L U:

> a = Numo::DFloat.new(4, 3).rand
=> Numo::DFloat#shape=[4,3]
[[0.0617545, 0.373067, 0.794815],
 [0.201042, 0.116041, 0.344032],
 [0.539948, 0.737815, 0.165089],
 [0.0508827, 0.108065, 0.0687079]]

> p, l, u = Numo::Linalg.lu(a)
=> [Numo::DFloat#shape=[4,4]
[[0, 1, 0, 0],
 [0, 0, 1, 0],
 [1, 0, 0, 0],
 [0, 0, 0, 1]], Numo::DFloat(view)#shape=[4,3]
[[1, 0, 0],
 [0.114371, 1, 0],
 [0.372336, -0.549651, 1],
 [0.0942363, 0.133488, -0.0711185]], Numo::DFloat(view)#shape=[3,3]
[[0.539948, 0.737815, 0.165089],
 [0, 0.288682, 0.775934],
 [0, 0, 0.709057]]]

> p.dot(l.dot(u))
=> Numo::DFloat#shape=[4,3]
[[0.0617545, 0.373067, 0.794815],
 [0.201042, 0.116041, 0.344032],
 [0.539948, 0.737815, 0.165089],
 [0.0508827, 0.108065, 0.0687079]]

> (a - p.dot(l.dot(u))).abs.max
=> 5.551115123125783e-17