Merge pull request #401 from avik-pal/ap/banded_qr

Banded QR for non square matrices

Author: ChrisRackauckas

.github/workflows/CI.yml

@@ -28,6 +28,8 @@ jobs:
             group: 'LinearSolveHYPRE'
           - version: '1'
             group: 'LinearSolvePardiso'
+          - version: '1'
+            group: 'LinearSolveBandedMatrices'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1

Project.toml

@@ -107,4 +107,4 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "IterativeSolvers", "InteractiveUtils", "JET", "KrylovKit", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff", "HYPRE", "MPI", "BlockDiagonals", "Enzyme", "FiniteDiff"]
+test = ["Test", "IterativeSolvers", "InteractiveUtils", "JET", "KrylovKit", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff", "HYPRE", "MPI", "BlockDiagonals", "Enzyme", "FiniteDiff", "BandedMatrices"]

ext/LinearSolveBandedMatricesExt.jl

@@ -5,8 +5,14 @@ import LinearSolve: defaultalg,
     do_factorization, init_cacheval, DefaultLinearSolver, DefaultAlgorithmChoice
 
 # Defaults for BandedMatrices
-function defaultalg(A::BandedMatrix, b, ::OperatorAssumptions)
-    return DefaultLinearSolver(DefaultAlgorithmChoice.DirectLdiv!)
+function defaultalg(A::BandedMatrix, b, oa::OperatorAssumptions)
+    if oa.issq
+        return DefaultLinearSolver(DefaultAlgorithmChoice.DirectLdiv!)
+    elseif LinearSolve.is_underdetermined(A)
+        error("No solver for underdetermined `A::BandedMatrix` is currently implemented!")
+    else
+        return DefaultLinearSolver(DefaultAlgorithmChoice.QRFactorization)
+    end
 end
 
 function defaultalg(A::Symmetric{<:Number, <:BandedMatrix}, b, ::OperatorAssumptions)

test/banded.jl

@@ -0,0 +1,38 @@
+using BandedMatrices, LinearAlgebra, LinearSolve, Test
+
+# Square Case
+n = 8
+A = BandedMatrix(Matrix(I, n, n), (2, 2))
+b = ones(n)
+A1 = A / 1;
+b1 = rand(n);
+x1 = zero(b);
+A2 = A / 2;
+b2 = rand(n);
+x2 = zero(b);
+
+sol1 = solve(LinearProblem(A1, b1; u0 = x1))
+@test sol1.u ≈ A1 \ b1
+sol2 = solve(LinearProblem(A2, b2; u0 = x2))
+@test sol2.u ≈ A2 \ b2
+
+# Square Symmetric
+A1s = Symmetric(A1)
+A2s = Symmetric(A2)
+
+sol1s = solve(LinearProblem(A1s, b1; u0 = x1))
+@test sol1s.u ≈ A1s \ b1
+sol2s = solve(LinearProblem(A2s, b2; u0 = x2))
+@test sol2s.u ≈ A2s \ b2
+
+# Underdetermined
+A = BandedMatrix(rand(8, 10), (2, 2))
+b = rand(8)
+
+@test_throws ErrorException solve(LinearProblem(A, b)).u
+
+# Overdetermined
+A = BandedMatrix(ones(10, 8), (2, 0))
+b = rand(10)
+
+@test_nowarn solve(LinearProblem(A, b))

test/runtests.jl

@@ -15,6 +15,7 @@ if GROUP == "All" || GROUP == "Core"
     @time @safetestset "Default Alg Tests" include("default_algs.jl")
     VERSION >= v"1.9" && @time @safetestset "Enzyme Derivative Rules" include("enzyme.jl")
     @time @safetestset "Traits" include("traits.jl")
+    VERSION >= v"1.9" && @time @safetestset "BandedMatrices" include("banded.jl")
 end
 
 if GROUP == "LinearSolveCUDA"