use LQ for "wide A" \ b, not QR

[stevengj]

As discussed in #34350, for a "wide" ("underdetermined") system of equations, A \ b seems to be more efficient if computed by "LQ" factorization, i.e. QR factorization of the transpose A'. (More than 2x faster on my machine for a 100×1000 matrix.)

• It still computes the minimum-norm solution as before, and should still be robust to ill-conditioned A (it doesn't square the condition number, and uses row-pivoting = column pivoting of the transpose). Using the LQ factorization to compute the minimum-norm solution is described in the LAPACK manual on LQ, is straightforward and is already implemented.
• The old solution uses pivoted QR of A (not its transpose) followed by "complete orthogonal factorization" of R. The latter step seems to pay a big performance price.
• Although we have lq(A), it does not currently support row-pivoting, so I use qr(A', ColumnNorm())' instead, which is logically equivalent.

Hopefully this is already covered by the existing tests?

[dkarrasch]

Although we have lq(A), it does not currently support row-pivoting, so I use qr(A', ColumnNorm())' instead

It turns out that this is (at least currently) the necessary way to go anyway: this works even for non-BLAS eltypes, whereas lq(A) fails due to

julia> lq(A) \ b
ERROR: MethodError: no method matching lq!(::Matrix{BigFloat})
The function `lq!` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  lq!(::StridedMatrix{var"#s4975"} where var"#s4975"<:Union{Float32, Float64, ComplexF64, ComplexF32})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/lq.jl:71

[dkarrasch]

The failing test is real, though the test itself is fragile. First, it would fail even without this PR if we had to rely on generic solves and could not use

julia/stdlib/LinearAlgebra/src/qr.jl

Line 546 in b6a2cc1

function ldiv!(A::QRPivoted{T,<:StridedMatrix}, B::AbstractMatrix{T}, rcond::Real) where {T<:BlasFloat}

:

julia> A = zeros(1, 2); B = zeros(1, 1); # from the failing testset

julia> A = big.(A); B = big.(B);

julia> qr(A, ColumnNorm()) \ B
ERROR: SingularException(1)

Second, we may need a similar implementation (which takes the computed rank into account) for AdjointFactorization{<:Any,<:QRPivoted}.