Multi-column matrix solve

[tel]

Solve currently implements Ax = b for known NxN matrix A and known N-dimensional vector b, resulting in a N-dimensional solution vector x. A common generalization is Ax = B for the same A, but with B now being a NxK dimensional matrix and x being solved as an NxK dimensional matrix itself.

This generalization is on one hand somewhat trivial as it can be performed via iterated solving for each of the K columns of B independently and could be reasonably implemented in Rust via that method. However, the LAPACK dgesv method already handles this generalization and thus it may be a small generalization of existing solving code to allow it.

As a practical application, this generalized solve is useful for computing Kalman filters where the optimal Kalman gain is computed as PH'/S where PH' is a full matrix whenever the observation space is multidimensional. While S is often well-behaved, it would be ideal to compute the gain via a solve routine as opposed to an invert-and-multiply.