I need to get any solution of a bunch of underdetermined systems `A*x = b`

fastly and cheaply. My first thought is to compute a minimum norm solution using the LQ decomposition of `A`

(which can be used as workspace to minimize allocations) with something like this:

```
using LinearAlgebra
m = 5; n = 10
A = randn(m, n)
b = randn(m)
x = similar(A, n)
ldiv!(x, lq!(A), b)
```

but this results in a `DimensionMismatch`

error, which I think is due to the square matrix `Q`

of the LQ decomposition.

Since I want to make this operation as cheap as possible, I considered calling `LAPACK.gels!`

directly, but the documentation says that `b`

is updated with the solution. How does this work in the minimum norm case, when the size of `b`

is less than the size of the solution `x`

?

Right now my best clue is to iterate over the LQ decomposition of `A`

(allocating more memory than necessary) and perform the matrix products manually:

```
using LinearAlgebra
m = 5; n = 10
A = randn(m, n)
b = randn(m)
L, Q = lq!(A)
ldiv!(LowerTriangular(L), b)
Q'b
```

Note: using a preallocated `x`

and trying `mul!(x, Q', b)`

also produces an `DimensionMismatch`

error.

Is there any better way than this?