JuliaLinearAlgebra / HierarchicalMatrices.jl

Julia package for hierarchical matrices
Other
24 stars 5 forks source link

2D with reshaped array #12

Open skeletonyk opened 6 years ago

skeletonyk commented 6 years ago

Really nice work! I am wondering if this is also suitable for the convolution of a 2D field? Do you require the indices of x be related to the physical position it represents? Thanks!

Ke

MikaelSlevinsky commented 6 years ago

Could you add a MWE that you would like to speed up with hierarchical matrices?

skeletonyk commented 6 years ago

Sorry a dumb question.. What is MWE? I am only solving a Poisson equation, using something called Lattice Green's Function which is a regularized discrete kernel which doesn't blow up at x_i = y_i. Thanks.

MikaelSlevinsky commented 6 years ago

MWE = minimal working example

skeletonyk commented 6 years ago

A naive version would be

  1. padding the matrix A. 2. FFT both the kernel and the padded matrix. 3. Dot product the two matrix. 4. IFFT Just a simple convolution using FFT. Does it answer the question regarding the MWE? Hopefully I am understanding it correctly~

Ke

MikaelSlevinsky commented 6 years ago

What does this have to do with hierarchical matrices? (How can you do convolution better than with the FFT?)

Cheers,

Mikael

On Mar 16, 2018, at 8:27 PM, Ke YU notifications@github.com<mailto:notifications@github.com> wrote:

A naive version would be

  1. padding the matrix A. 2. FFT both the kernel and the padded matrix. 3. Dot product the two matrix. 4. IFFT Just a simple convolution using FFT. Does it answer the question regarding the MWE? Hopefully I am understanding it correctly~

Ke

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/JuliaMatrices/HierarchicalMatrices.jl/issues/12#issuecomment-373883813, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AHzBpT5bbJ2ynT9Ds0veCmp0JfkMtM-Gks5tfGaNgaJpZM4StAK6.

skeletonyk commented 6 years ago

Sorry for the late reply. Maybe I understood it wrong. I thought this package is similar to a multipole algorithm where the hierarchical matrices are referring the quadtree or octree structure.. Sorry for the interruption if I misunderstood the package.

MikaelSlevinsky commented 6 years ago

It implements precisely that!

But there are many different applications/uses which is why if you show me a MWE (in Julia code) then I might be able to identify the precise hierarchy you are referring to/need.

Cheers,

Mikael

On Mar 20, 2018, at 5:56 PM, Ke YU notifications@github.com<mailto:notifications@github.com> wrote:

Sorry for the late reply. Maybe I understood it wrong. I thought this package is similar to a multipole algorithm where the hierarchical matrices are referring the quadtree or octree structure.. Sorry for the interruption if I misunderstood the package.

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/JuliaMatrices/HierarchicalMatrices.jl/issues/12#issuecomment-374785298, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AHzBpcHH2AbH26EAVxo5fp4WV28DDQleks5tgYkigaJpZM4StAK6.

skeletonyk commented 6 years ago

I don't have a Julia code for that. Our lab has a Fortran version of the following algorithm and I was thinking of implement a lightweight version of it in Julia. The method is described in the paper, 'Liska, Sebastian, and Tim Colonius. "A parallel fast multipole method for elliptic difference equations." Journal of Computational Physics 278 (2014): 76-91.'

Maybe it's too detailed to require this package to do. Basically It solves the Poisson equation on a Cartesian grid using polynomial interpolation and anterpolation to bridge the different hierarchies. Do you have any advices on how this can be most easily incorporated into the framework of this package maybe?

Ke

skeletonyk commented 6 years ago

If I could briefly introduce this method. The domain is divided into many 'blocks'. Each block contains n^3 cells. n is typically a small integer like 10. When coarsened, 8 blocks (3d) are merged into a bigger block which also have n^3 cells but of twice the size. The values of the cells of the bigger blocks are interpolated using polynomials from the small cells.

I am also open to other approaches for solving Poisson equations on regular Cartesian grids. Do you have any suggestions?

maltezfaria commented 5 years ago

Hi,

I was wondering if it is possible/easy to modify the example of the kernel matrix in 1d to take a point say in 2 or 3d instead. Something like generating a hierarchical representation of the matrix A below:

using LinearAlgebra
N = 20
u = range(-1,1,length=N)
v = range(-1,1,length=N)
pts = Vector{NTuple{3,Float64}}()
for u in u, v in v
    x = sinpi(u/2)*cospi(v)
    y = sinpi(u/2)*cospi(v)
    z = cospi(u/2)
    push!(pts,(x,y,z))
end
f(x,y) = inv(norm(x.-y))
f(pts[1],pts[2])
A = [f(pt1,pt2) for pt1 in pts, pt2 in pts]

I think it is somewhat related to this open ticket, so I decided to post here. Thanks!

MikaelSlevinsky commented 5 years ago

Hi @maltezfaria, it appears your points are all on the unit 2-sphere, neat! I think that currently the data structures can handle the particular matrix A, but there are no constructors yet. Such constructors would be available if the algorithms described in Introduction to hierarchical matrices with applications were implemented.