Add a `linalg.lu` function for the LU decomposition

[rgommers]

Only the default (partial pivoting) algorithm that is implemented in all libraries and for all devices is added here. gh-627 has details on the no-pivoting case, but it's not universally supported and the only reason to add it would be that it's more performant in some cases where users know it will be numerically stable. Such an addition can be done in the future, but it seems like a potentially large amount of work for implementers for limited gain.

Closes gh-627

EDIT: I'd like to look a bit more closely at how potentially to support this in NumPy, so converted the PR to draft for now.