Iterative solver vs. FFT
Joaquin Drut
The laplacian is diagonal in Fourier space. You can Fourier transform
the rhs, multiply by the inverse laplacian and FT back. This is
typically faster than an iterative solver, at least for quantum
mechanics applications where one uses a multipoint laplacian.
Has anyone tried this approach on a GPU? Is there any reason not to
try this? (e.g. FFT is a highly non-local operation).
Cheers,
Joaquin.
Liwen Chang
Joaquin Drut
Liwen Chang
Joaquin Drut
wyang
> Also, the matrix shown in the slides actually has no well-defined
> inverse (it has a zero mode, because the laplacian is translation
> invariant).
you correctly pointed out, if p1(x) is a solution, p1(x)+p0 is also a
solution where p0 is a constant. In case of Fourier transform, p0 will
be added to the k=0 term (the constant term).
For this problem, I would anchor the value at a grid point. For
example, I would enforce P[0] to be 0. People enforce a boundary
condition in many ways. One way is to explicitly eliminate P[0] (and
P[N-1]) from the matrix equation, which means one needs to reformulate
to solve an (N-2)x(N-2) matrix equation instead of NxN (since
P[N-1]=P[0]=0). Another method is a "penalty method" (which I learned
from a FEM course).
One can solve a pde in various ways including using FFT. One reason
that people use grid-based methods (FDM, FEM, FVM, etc.), despite some
strong benefits with the Fourier transform method, is that operators
are local, so we don't do lots of global reductions. Another point
people often quote is that grid-based approaches are generally
superior where there are some sort of discontinuities because it would
be impossible to represent such a discontinuity in terms of waves.
On Aug 4, 11:36 am, Joaquin Drut <joaquind...@gmail.com> wrote:
> I'm sorry, I failed to understand your message. My question was about
> how to regularize the laplacian shown in the slides. How would you use
> iterative methods for that?
>
> Indeed, that discretization of the laplacian shown is a special case,
> but it appears in many problems. That's why I was wondering if there
> is a particularly good way to regularize the operator when solving the
> problem on a GPU. Other multipoint discretizations have the same
> problem (again b/c they are all translation invariant).
>
> Cheers,
> Joaquin.
>
> On Aug 4, 2010, at 2:25 PM, Liwen Chang wrote:
>
> > In that case,
> > I probably would use iterative methods.
> > and BTW, Laplacian is just a case of tridiagonal matrices.
> > Mr.Cohen's talk can be applied for any tridiagonal matrix
>
> > Liwen
>
Joaquin Drut
On Aug 6, 2010, at 12:56 AM, wyang wrote:
> I am speaking from my old memory, so details may be incorrect.
>
>> Also, the matrix shown in the slides actually has no well-defined
>> inverse (it has a zero mode, because the laplacian is translation
>> invariant).
>
> I don't think FFT is free from the non-uniqueness problem, either. As
> you correctly pointed out, if p1(x) is a solution, p1(x)+p0 is also a
> solution where p0 is a constant. In case of Fourier transform, p0 will
> be added to the k=0 term (the constant term).
This is correct. I didn't mean to imply that FFT will save you from
the fact that the "bare" laplacian is a singular operator.
> For this problem, I would anchor the value at a grid point. For
> example, I would enforce P[0] to be 0. People enforce a boundary
> condition in many ways. One way is to explicitly eliminate P[0] (and
> P[N-1]) from the matrix equation, which means one needs to reformulate
> to solve an (N-2)x(N-2) matrix equation instead of NxN (since
> P[N-1]=P[0]=0). Another method is a "penalty method" (which I learned
> from a FEM course).
Ok, now we are talking. This is a basic implementation of something
more general that in physics we call "gauge fixing". It makes sense.
The boundary conditions of the laplacian shown in the slides were
periodic, which was the source of the problem. If you choose
antiperiodic, or if you are interested in some other kind of boundary
condition, then you might be safe without specifying the constant or
anchoring point you mentioned above. Otherwise you would have an
infinite number of solutions and the iterative solver would run forever.
> One can solve a pde in various ways including using FFT. One reason
> that people use grid-based methods (FDM, FEM, FVM, etc.), despite some
> strong benefits with the Fourier transform method, is that operators
> are local, so we don't do lots of global reductions. Another point
> people often quote is that grid-based approaches are generally
> superior where there are some sort of discontinuities because it would
> be impossible to represent such a discontinuity in terms of waves.
Ok, I can see your points. I'm interested in these issues because we
use FFT a lot in Quantum Mechanics simulations. Part of the reason is
that we use a full representation of the laplacian, i.e. we use ALL
the points in the grid to represent it... I guess that could be called
a "dense" stencil operator. In that case you have extreme non-locality
anyway, so FFT is much better than matrix-vector multiplication.
In any case, thanks for the comments!