JuliaApproximation / PiecewiseOrthogonalPolynomials.jl

A Julia package for piecewise spectral methods such as p-FEM
MIT License
6 stars 2 forks source link

Introduce ArrowheadMatrix #32

Closed dlfivefifty closed 1 year ago

codecov[bot] commented 1 year ago

Codecov Report

Patch coverage: 98.75% and project coverage change: +1.30 :tada:

Comparison is base (b8486b9) 97.40% compared to head (53d307a) 98.70%.

Additional details and impacted files ```diff @@ Coverage Diff @@ ## main #32 +/- ## ========================================== + Coverage 97.40% 98.70% +1.30% ========================================== Files 1 2 +1 Lines 154 386 +232 ========================================== + Hits 150 381 +231 - Misses 4 5 +1 ``` | [Impacted Files](https://app.codecov.io/gh/JuliaApproximation/PiecewiseOrthogonalPolynomials.jl/pull/32?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaApproximation) | Coverage Δ | | |---|---|---| | [src/PiecewiseOrthogonalPolynomials.jl](https://app.codecov.io/gh/JuliaApproximation/PiecewiseOrthogonalPolynomials.jl/pull/32?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaApproximation#diff-c3JjL1BpZWNld2lzZU9ydGhvZ29uYWxQb2x5bm9taWFscy5qbA==) | `98.41% <97.72%> (+1.01%)` | :arrow_up: | | [src/arrowhead.jl](https://app.codecov.io/gh/JuliaApproximation/PiecewiseOrthogonalPolynomials.jl/pull/32?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaApproximation#diff-c3JjL2Fycm93aGVhZC5qbA==) | `98.98% <98.98%> (ø)` | |

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.

dlfivefifty commented 1 year ago

@ioannisPApapadopoulos @KarsKnook Got an efficient reverse Cholesky working with arrowhead matrices! Here we can see that with two elements and p = 100k we can expect it to cost around 0.007662s:

julia> n = 3; p = 100_000;

julia> A = ArrowheadMatrix(BandedMatrix(0 => randn(n) .+ 10, 1 => randn(n-1), -1 => randn(n-1)), 
                               [BandedMatrix((0 => randn(n-1), -1 => randn(n-1)), (n,n-1)) for _=1:2],
                               BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}[],
                            fill(BandedMatrix((0 => randn(p) .+ 10, 2 => randn(p-2)), (p, p)), n-1));

julia> @time reversecholesky!(Symmetric(copy(A)));
  0.007662 seconds (18 allocations: 4.579 MiB)

It will be significantly faster if we exploit the fact that the operators are the same in each element (this might even happen automatically if we replace fill with Fill but if not its an easy change).

So just need the correct formulae for weak Laplacian / mass matrix