From 8d6494c623fa163410968278b31d6a72530e61c6 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Thu, 18 Apr 2024 11:58:00 -0400
Subject: [PATCH] Change the default QR to ColumnNorm

---
 Project.toml         |  2 +-
 src/LinearSolve.jl   |  1 +
 src/default.jl       | 14 ++------------
 test/default_algs.jl | 34 +++++++++++++++++++++++++++++++++-
 4 files changed, 37 insertions(+), 14 deletions(-)

diff --git a/Project.toml b/Project.toml
index d4295eecd..89602a87d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "LinearSolve"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 authors = ["SciML"]
-version = "2.28.0"
+version = "2.29.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl
index e274521c9..a89d1ff50 100644
--- a/src/LinearSolve.jl
+++ b/src/LinearSolve.jl
@@ -70,6 +70,7 @@ needs_concrete_A(alg::AbstractKrylovSubspaceMethod) = false
 needs_concrete_A(alg::AbstractSolveFunction) = false
 
 # Util
+## This is a check exclusively based on the size and not on the actual rank of the matrix
 is_underdetermined(x) = false
 is_underdetermined(A::AbstractMatrix) = size(A, 1) < size(A, 2)
 is_underdetermined(A::AbstractSciMLOperator) = size(A, 1) < size(A, 2)
diff --git a/src/default.jl b/src/default.jl
index 4621edcff..47fe280d7 100644
--- a/src/default.jl
+++ b/src/default.jl
@@ -216,19 +216,9 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool})
     elseif assump.condition === OperatorCondition.WellConditioned
         DefaultAlgorithmChoice.NormalCholeskyFactorization
     elseif assump.condition === OperatorCondition.IllConditioned
-        if is_underdetermined(A)
-            # Underdetermined
-            DefaultAlgorithmChoice.QRFactorizationPivoted
-        else
-            DefaultAlgorithmChoice.QRFactorization
-        end
+        DefaultAlgorithmChoice.QRFactorizationPivoted
     elseif assump.condition === OperatorCondition.VeryIllConditioned
-        if is_underdetermined(A)
-            # Underdetermined
-            DefaultAlgorithmChoice.QRFactorizationPivoted
-        else
-            DefaultAlgorithmChoice.QRFactorization
-        end
+        DefaultAlgorithmChoice.QRFactorizationPivoted
     elseif assump.condition === OperatorCondition.SuperIllConditioned
         DefaultAlgorithmChoice.SVDFactorization
     else
diff --git a/test/default_algs.jl b/test/default_algs.jl
index 8d3ac261e..d9fb79141 100644
--- a/test/default_algs.jl
+++ b/test/default_algs.jl
@@ -33,7 +33,7 @@ solve(prob)
 
 @test LinearSolve.defaultalg(nothing, zeros(5),
     LinearSolve.OperatorAssumptions(false)).alg ===
-      LinearSolve.DefaultAlgorithmChoice.QRFactorization
+      LinearSolve.DefaultAlgorithmChoice.QRFactorizationPivoted
 
 A = spzeros(2, 2)
 # test that solving a singular problem doesn't error
@@ -112,3 +112,35 @@ prob = LinearProblem(funcop, b)
 sol1 = solve(prob)
 sol2 = solve(prob, LinearSolve.KrylovJL_CRAIGMR())
 @test sol1.u == sol2.u
+
+# Default for Underdetermined problem but the size is a long rectangle
+A = [2.0 1.0
+     0.0 0.0
+     0.0 0.0]
+b = [1.0, 0.0, 0.0]
+prob = LinearProblem(A, b)
+sol = solve(prob)
+
+@test sol.u ≈ [0.4, 0.2]
+
+## Show that we cannot select a default alg once by checking the rank, since it might change
+## later in the cache
+## Common occurrence for iterative nonlinear solvers using linear solve
+A = [2.0 1.0
+     1.0 1.0
+     0.0 0.0]
+b = [1.0, 1.0, 0.0]
+prob = LinearProblem(A, b)
+
+cache = init(prob)
+sol = solve!(cache)
+
+@test sol.u ≈ [0.0, 1.0]
+
+cache.A = [2.0 1.0
+           0.0 0.0
+           0.0 0.0]
+
+sol = solve!(cache)
+
+@test sol.u ≈ [0.4, 0.2]