From 4ad6fc7d4be3e72a83428ed53665b854c4093414 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 2 Jun 2025 22:01:31 +0530 Subject: [PATCH 1/5] Add banded paths to generic solve --- src/generic.jl | 12 ++++++++++++ src/tridiag.jl | 5 +++++ test/generic.jl | 12 ++++++++++++ 3 files changed, 29 insertions(+) diff --git a/src/generic.jl b/src/generic.jl index c8336f87..f5b7e695 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -1239,13 +1239,25 @@ function (\)(A::AbstractMatrix, B::AbstractVecOrMat) if istril(A) if istriu(A) return Diagonal(A) \ B + elseif istriu(A, -1) + return Bidiagonal(A, :L) \ B else return LowerTriangular(A) \ B end end if istriu(A) + if istril(A, 1) + return Bidiagonal(A, :U) \ B + end return UpperTriangular(A) \ B end + if istriu(A, -1) && istril(A, 1) + T = Tridiagonal(A) + if issymmetric(T) + return SymTridiagonal(T) \ B + end + return T \ B + end return lu(A) \ B end return qr(A, ColumnNorm()) \ B diff --git a/src/tridiag.jl b/src/tridiag.jl index 0e8a5119..80078aeb 100644 --- a/src/tridiag.jl +++ b/src/tridiag.jl @@ -1119,6 +1119,11 @@ function ldiv!(A::Tridiagonal, B::AbstractVecOrMat) return B end +function Base.:(\)(A::Tridiagonal, B::AbstractVecOrMat) + T = Base.promote_op(\, eltype(A), eltype(B)) + ldiv!(T.(A), copy(B)) +end + # combinations of Tridiagonal and Symtridiagonal # copyto! for matching axes function _copyto_banded!(A::Tridiagonal, B::SymTridiagonal) diff --git a/test/generic.jl b/test/generic.jl index 521beaf4..59760db0 100644 --- a/test/generic.jl +++ b/test/generic.jl @@ -937,4 +937,16 @@ end @test B == A2 end +@testset "linear solve for dense banded matrices" begin + b = randn(5) + for A in (diagm(0=>1:5, -1=>1:4, 1=>1:4), # SymTridiagonal + diagm(0=>1:5, -1=>1:4, 1=>5:8), # Tridiagonal + diagm(0=>1:5, 1=>1:4), # Upper Bidiagonal + diagm(0=>1:5, -1=>1:4), # Lower Bidiagonal + ) + x = A \ b + @test A * x ≈ b + end +end + end # module TestGeneric From 4e3eabfb34e7affc4632ec9fff178d1a15910402 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 2 Jun 2025 22:28:50 +0530 Subject: [PATCH 2/5] Remove tridiag method --- src/generic.jl | 5 ++--- src/tridiag.jl | 5 ----- 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/src/generic.jl b/src/generic.jl index f5b7e695..91b412fd 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -1251,12 +1251,11 @@ function (\)(A::AbstractMatrix, B::AbstractVecOrMat) end return UpperTriangular(A) \ B end - if istriu(A, -1) && istril(A, 1) - T = Tridiagonal(A) + if isbanded(A, -1, 1) + T = Tridiagonal(diagview(A,-1), diagview(A, 0), diagview(A, 1)) if issymmetric(T) return SymTridiagonal(T) \ B end - return T \ B end return lu(A) \ B end diff --git a/src/tridiag.jl b/src/tridiag.jl index 80078aeb..0e8a5119 100644 --- a/src/tridiag.jl +++ b/src/tridiag.jl @@ -1119,11 +1119,6 @@ function ldiv!(A::Tridiagonal, B::AbstractVecOrMat) return B end -function Base.:(\)(A::Tridiagonal, B::AbstractVecOrMat) - T = Base.promote_op(\, eltype(A), eltype(B)) - ldiv!(T.(A), copy(B)) -end - # combinations of Tridiagonal and Symtridiagonal # copyto! for matching axes function _copyto_banded!(A::Tridiagonal, B::SymTridiagonal) From 0c995d776e6245a9f19c1dbaea47cc9a522b8fdf Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 2 Jun 2025 22:34:58 +0530 Subject: [PATCH 3/5] Use `lu(::Tridiagonal)` in generic solve --- src/generic.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/generic.jl b/src/generic.jl index 91b412fd..a739d803 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -1256,6 +1256,7 @@ function (\)(A::AbstractMatrix, B::AbstractVecOrMat) if issymmetric(T) return SymTridiagonal(T) \ B end + return lu(T) \ B end return lu(A) \ B end From c19ebaa721e03918096e251ca9d77e4c4f6856ed Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 2 Jun 2025 23:08:49 +0530 Subject: [PATCH 4/5] Materialize `Triangular` matrix --- src/generic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generic.jl b/src/generic.jl index a739d803..f6b42c9c 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -1252,7 +1252,7 @@ function (\)(A::AbstractMatrix, B::AbstractVecOrMat) return UpperTriangular(A) \ B end if isbanded(A, -1, 1) - T = Tridiagonal(diagview(A,-1), diagview(A, 0), diagview(A, 1)) + T = Tridiagonal(A) if issymmetric(T) return SymTridiagonal(T) \ B end From f348c955fbc2da52b20c755dc9f9e4ab8458ba6a Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 3 Jun 2025 23:52:30 +0530 Subject: [PATCH 5/5] Restrict to `Bidiagonal` --- src/generic.jl | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/src/generic.jl b/src/generic.jl index f6b42c9c..83c2c1c8 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -1236,28 +1236,24 @@ function (\)(A::AbstractMatrix, B::AbstractVecOrMat) require_one_based_indexing(A, B) m, n = size(A) if m == n - if istril(A) - if istriu(A) - return Diagonal(A) \ B - elseif istriu(A, -1) + istrium1 = istriu(A, -1) + istril1 = istril(A, 1) + if istril1 && iszero(diagview(A,1)) # istril(A) + if istrium1 + if iszero(diagview(A, -1)) + return Diagonal(A) \ B + end return Bidiagonal(A, :L) \ B else return LowerTriangular(A) \ B end end - if istriu(A) - if istril(A, 1) + if istrium1 && iszero(diagview(A, -1)) # istriu(A) + if istril1 return Bidiagonal(A, :U) \ B end return UpperTriangular(A) \ B end - if isbanded(A, -1, 1) - T = Tridiagonal(A) - if issymmetric(T) - return SymTridiagonal(T) \ B - end - return lu(T) \ B - end return lu(A) \ B end return qr(A, ColumnNorm()) \ B