Default non-square operators to LSMR

[ChrisRackauckas]

Is there any better choice with no other information a priori?

[codecov]

Codecov Report

Merging #188 (ef07182) into main (79eea00) will decrease coverage by 1.25%.
The diff coverage is 71.42%.

@@            Coverage Diff             @@
##             main     #188      +/-   ##
==========================================
- Coverage   64.25%   62.99%   -1.26%     
==========================================
  Files           9        9              
  Lines         621      627       +6     
==========================================
- Hits          399      395       -4     
- Misses        222      232      +10     

Impacted Files	Coverage Δ
src/LinearSolve.jl	0.00% <ø> (-75.00%)	⬇️
src/default.jl	40.90% <0.00%> (-7.53%)	⬇️
src/factorization.jl	83.25% <100.00%> (+1.16%)	⬆️
src/iterative_wrappers.jl	70.81% <100.00%> (+1.41%)	⬆️
src/common.jl	78.94% <0.00%> (-8.78%)	⬇️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[KlausC]

Obtained a list of potential Krylov- solvers, which can be called with non-square matrices,
there is no dependency on sign(m - n):

solver cgls: ok
solver cgne: ok
solver craig: ok
solver craigmr: ok
solver crls: ok
solver crmr: ok
solver lnlq: ok
solver lslq: ok
solver lsmr: ok
solver lsqr: ok

Code:

julia> using Krylov

julia> function checkcallable(s::Symbol, args...)
           st = string(s)
           islowercase(first(st)) && last(st) == '!' || return
           s = Symbol(st[1:end-1])
           isdefined(Krylov, s) || return
           f = getproperty(Krylov, s)
           try
               res = f(args...)
               println("solver $f: $(res isa Tuple && length(res) >= 2 ? "ok" : "nok")")
           catch ex
               if ex isa MethodError
                   println("solver $f(A, b) not callable")
               else
                   println("solver $f: failed with $ex")
               end
           end
       end
checkcallable (generic function with 1 method)

julia> m = 2; n = 12; A = rand(m, n); b = rand(m);

julia> checkcallable.(names(Krylov), Ref(A), Ref(b));
...

[KlausC]

Still see following issue:

julia> m, n = 13, 3
(13, 3)

julia> A = sprand(m, n, 0.5); b = rand(m);

julia> prob = LinearProblem(A, b)
LinearProblem. In-place: true
...

julia> solve(prob, KrylovJL_LSMR())
u: 3-element Vector{Float64}:
 1.7737328916933053
 0.3025612503302547
 0.4379756520127739

julia> solve(prob)
ERROR: BoundsError
Stacktrace:
 ...

That is:

1. default algorithm (A sparse, m > n) is still QRFactorization()
2. with alg = QRFactorization() there is still the BoundsError in _ldiv
3. Test cases for (sparse, full) x ( m < n, m = n , m > n) missing