Skip to content

Add banded paths to generic solve #1372

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Conversation

jishnub
Copy link
Member

@jishnub jishnub commented Jun 2, 2025

Currently, the generic solve only checks for diagonal and triangular matrices when dispatching to the faster methods. We also have solve methods for Bidiagonal and SymTridiagonal matrices, so we may add these branches as well.

With these,

julia> n = 5000; dv, ev = rand(n), rand(n-1);

julia> A = diagm(0=>dv, 1=>ev, -1=>ev); # symmetric tridiagonal by construction

julia> @btime $A \ $dv;
  656.843 ms (9 allocations: 190.81 MiB) # master
  11.549 ms (22 allocations: 235.09 KiB) # PR

After this change, almost all the time is spent on checking if the matrix is tridiagonal.

@stevengj
Copy link
Member

stevengj commented Jun 2, 2025

(This all seems a bit weird to me; I feel like if you want a structured solve algorithm, you should use a structured-matrix type. But I guess since we're already checking for triangular matrices there is some argument for this…)

@jishnub
Copy link
Member Author

jishnub commented Jun 3, 2025

I see this as a generic fallback, where the structure may not be obvious from the type, but may be obtained by looking at the values. Suppose we have a fast solver for a Diagonal, but we call the method with a Transpose of a Diagonal, we may still use the banded structure even when we aren't directly dispatching to the method for the Diagonal. One way around this is to define \ for a Transpose{<:Any, <:Diagonal}, which would also work. The alternate, as in this PR, is to check for the banded structure from the values in the fallback method, and dispatch to the fast paths. This makes sense when the alternate approach --- factorizing the matrix --- is significantly more expensive than checking for a banded structure. In addition, this method would now work with arbitrary banded matrices that have 3 or fewer populated bands, even when the package might not be implementing a fast solver for these.

Copy link

codecov bot commented Jun 3, 2025

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.81%. Comparing base (c64c328) to head (f348c95).

Additional details and impacted files
@@           Coverage Diff           @@
##           master    #1372   +/-   ##
=======================================
  Coverage   93.81%   93.81%           
=======================================
  Files          34       34           
  Lines       15717    15723    +6     
=======================================
+ Hits        14745    14751    +6     
  Misses        972      972           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

return Bidiagonal(A, :U) \ B
end
return UpperTriangular(A) \ B
end
if isbanded(A, -1, 1)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't you get this "for free" with istrium1 && istril1?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, but there were some triangular lu-related test failures for unusual types, so I've removed the triangular branch for now. We may add this back once those issues are resolved.

jishnub added a commit that referenced this pull request Jun 3, 2025
Working on #1372, I came across some unusual tests with non-Number
`eltype`s, for which this `fill!` fails unless the appropriate zero is
used. Probably missing convert methods, but we might as well use the
zero of the eltype here.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants