JuliaApproximation / PiecewiseOrthogonalPolynomials.jl

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

QL might need strong form for optimal complexity #48

Open dlfivefifty opened 11 months ago

dlfivefifty commented 11 months ago

@ioannisPApapadopoulos @karsknook

I've been playing with QL instead of reverse Cholesky for non-SPD operators. It suffers from fill in in the top row of blocks so the complexity would be $O(np+n^3)$.

BUT if we did strong form (so non-symmetric matrices) I believe everything apart from the first block column will be diagonal. Then QL will be optimal complexity.

Note what we need is ContinuousPolynomial{-1} which would be P^(1,1) instead of Legendre but also with delta functions between the elements.

KarsKnook commented 11 months ago

I am assuming you mean symmetric indefinite matrices when you say non-SPD. John and I agreed a couple weeks ago that in that case you can still do reverse LU and it will give sparse factors. This was relevant for his cylindrical case where he was using ContinuousPolynomial{1} only in the z-direction. Obviously when you do tensor products you get overlapping spectra so Fortunato & Townsend's ADI result fails.

dlfivefifty commented 11 months ago

You can but it will be unstable without pivoting. You want QL for reliable solves

ioannisPApapadopoulos commented 11 months ago

In general I agree but I've had good success without pivoting with the disk/annulus stuff. It worked for the complex valued matrices too :)

KarsKnook commented 11 months ago

What would a pivoted reverse LU look like? If I remember correctly, the entries are largest on the diagonal and increase with $p$, so I don't think there will be any pivoting until you reach the lowest order blocks

KarsKnook commented 11 months ago

Partial pivoting that is (I am not sure what kind of pivoting is usually preferred)

dlfivefifty commented 11 months ago

In general I agree but I've had good success without pivoting with the disk/annulus stuff. It worked for the complex valued matrices too :)

Works for general potentials ? (Yeah right! There's a reason people invented pivoting.)

If I remember correctly, the entries are largest on the diagonal and increase with , so I don't think there will be any pivoting until you reach the lowest order blocks

Depends on the potential. The model problem is high frequency Helmholtz where diagonal dominance is lost. (Its possible its still stable without pivoting for constant coefficients but with other potentials this won't be the case.)

What would a pivoted reverse LU look like?

Pivoting will likely have the exact same problems with fill-in as QL. (Note "reverse LU" = "UL")