From eb1da3d41d29d7f985bf85d7d100c659fcc8f582 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 7 Nov 2023 15:34:07 -0500 Subject: [PATCH 1/3] Add support for almost banded matrix --- Project.toml | 5 ++- ext/LinearSolveFastAlmostBandedMatricesExt.jl | 33 +++++++++++++++++++ 2 files changed, 37 insertions(+), 1 deletion(-) create mode 100644 ext/LinearSolveFastAlmostBandedMatricesExt.jl diff --git a/Project.toml b/Project.toml index 65e3ff9b5..5283d1361 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.19.0" +version = "2.20.0" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" @@ -34,6 +34,7 @@ BandedMatrices = "aae01518-5342-5314-be14-df237901396f" BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" +FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e" HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" @@ -54,6 +55,7 @@ LinearSolveKrylovKitExt = "KrylovKit" LinearSolveMetalExt = "Metal" LinearSolvePardisoExt = "Pardiso" LinearSolveRecursiveArrayToolsExt = "RecursiveArrayTools" +LinearSolveFastAlmostBandedMatricesExt = ["FastAlmostBandedMatrices"] [compat] ArrayInterface = "7.4.11" @@ -92,6 +94,7 @@ julia = "1.9" BandedMatrices = "aae01518-5342-5314-be14-df237901396f" BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0" Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" +FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" diff --git a/ext/LinearSolveFastAlmostBandedMatricesExt.jl b/ext/LinearSolveFastAlmostBandedMatricesExt.jl new file mode 100644 index 000000000..187f2b454 --- /dev/null +++ b/ext/LinearSolveFastAlmostBandedMatricesExt.jl @@ -0,0 +1,33 @@ +module LinearSolveFastAlmostBandedMatricesExt + +using FastAlmostBandedMatrices, LinearAlgebra, LinearSolve +import LinearSolve: defaultalg, + do_factorization, init_cacheval, DefaultLinearSolver, DefaultAlgorithmChoice + +function defaultalg(A::AlmostBandedMatrix, b, oa::OperatorAssumptions) + if oa.issq + return DefaultLinearSolver(DefaultAlgorithmChoice.DirectLdiv!) + else + return DefaultLinearSolver(DefaultAlgorithmChoice.QRFactorization) + end +end + +# For BandedMatrix +for alg in (:SVDFactorization, :MKLLUFactorization, :DiagonalFactorization, + :SparspakFactorization, :KLUFactorization, :UMFPACKFactorization, + :GenericLUFactorization, :RFLUFactorization, :BunchKaufmanFactorization, + :CHOLMODFactorization, :NormalCholeskyFactorization, :LDLtFactorization, + :AppleAccelerateLUFactorization, :CholeskyFactorization, :LUFactorization) + @eval begin + function init_cacheval(::$(alg), ::AlmostBandedMatrix, b, u, Pl, Pr, maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + return nothing + end + end +end + +function do_factorization(alg::QRFactorization, A::AlmostBandedMatrix, b, u) + return alg.inplace ? qr!(A) : qr(A) +end + +end From 24e618b924690b66b18b98afe6c614aaa8683a03 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Wed, 8 Nov 2023 15:02:44 -0500 Subject: [PATCH 2/3] Fallback for non-square A --- src/factorization.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/factorization.jl b/src/factorization.jl index a0c7abf71..786b202c3 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -12,6 +12,8 @@ _ldiv!(x, A, b) = ldiv!(x, A, b) function _ldiv!(x::Vector, A::Factorization, b::Vector) # workaround https://github.com/JuliaLang/julia/issues/43507 + # Fallback if working with non-square matrices + length(x) != length(b) && return ldiv!(x, A, b) copyto!(x, b) ldiv!(A, x) end From e7b1f7d43458b6253db48610f36c23e1e4cd595a Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Sun, 12 Nov 2023 10:25:07 -0500 Subject: [PATCH 3/3] Add tests --- Project.toml | 2 +- test/banded.jl | 18 +++++++++++++++++- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 5283d1361..71f0ca822 100644 --- a/Project.toml +++ b/Project.toml @@ -113,4 +113,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", "BandedMatrices"] +test = ["Test", "IterativeSolvers", "InteractiveUtils", "JET", "KrylovKit", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff", "HYPRE", "MPI", "BlockDiagonals", "Enzyme", "FiniteDiff", "BandedMatrices", "FastAlmostBandedMatrices"] diff --git a/test/banded.jl b/test/banded.jl index 0bbf5e65b..89c104966 100644 --- a/test/banded.jl +++ b/test/banded.jl @@ -1,4 +1,4 @@ -using BandedMatrices, LinearAlgebra, LinearSolve, Test +using FastAlmostBandedMatrices, BandedMatrices, LinearAlgebra, LinearSolve, Test # Square Case n = 8 @@ -16,6 +16,12 @@ sol1 = solve(LinearProblem(A1, b1; u0 = x1)) sol2 = solve(LinearProblem(A2, b2; u0 = x2)) @test sol2.u ≈ A2 \ b2 +A = AlmostBandedMatrix(BandedMatrix(fill(2.0, n, n), (1, 1)), fill(3.0, 2, n)) +A[band(0)] .+= 1:n + +sol1ab = solve(LinearProblem(A, b; u0 = x1)) +@test sol1ab.u ≈ Matrix(A) \ b + # Square Symmetric A1s = Symmetric(A1) A2s = Symmetric(A2) @@ -31,8 +37,18 @@ b = rand(8) @test_throws ErrorException solve(LinearProblem(A, b)).u +A = AlmostBandedMatrix(BandedMatrix(fill(2.0, n - 2, n), (1, 1)), fill(3.0, 2, n)) +A[band(0)] .+= 1:(n - 2) + +@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)) + +A = AlmostBandedMatrix(BandedMatrix(fill(2.0, n + 2, n), (1, 1)), fill(3.0, 2, n)) +A[band(0)] .+= 1:n + +@test_nowarn solve(LinearProblem(A, b))