## Fast direct solver for the 2D-Poisson's equation in C++

#### About

- Direct solver for the 2D-Poisson's equation, $u_{xx} + u_{yy} = f$, based on the Fast Fourier Transform, using the FFTW library.
- Dirichlet and Neumann boundary conditions are supported.
- Dirichlet boundary conditions lead to DST-I and Neumann lead to DCT-I, but only if the boundary derivatives are approximated by 2-non-adjacent grid points, i.e. $g'(x_0) \approx \frac{g(x_2)-g(x_0)}{2 h}$. Although unusual, this has been necessary in order to be able to use the FFT for solving the resulting equation system.
- Only uniform grids are supported, and grid values are stored in a 2D array using boost::multi_array.

#### Download Source Code

The source code is released under the GPLv2 or above:- source: laplace.h, laplace.cpp
- other: Makefile, example.cpp, reg_test.cpp

#### Dependencies

FFTW library and boost::multi_array.#### Usage example

#include <cstdio> #include <boost/multi_array.hpp> #include "laplace.h" int main(int argc, char** argv) { size_t n1=5472, n2=3648; double h1=1.0, h2=1.0, a1=1.0, a2=1.0; pde::types::boundary bdtype=pde::types::Neumann; double bdvalue=0.0; pde::fftw_threads(4); // number of threads for internal fft routine boost::multi_array<double,2> F(boost::extents[n1][n2]); boost::multi_array<double,2> U; double trunc; // return value for poisolve(), indicates truncation error arr::runif(F);Compile:if(bdtype==pde::types::Neumann) { // find compatible boundary value, otherwise truncation error>0 bdvalue=pde::neumann_compat(F,a1,a2,h1,h2); } // Poisson equation solver trunc=pde::poisolve(U,F,a1,a2,h1,h2,bdvalue,bdtype,false); // ... // minor cleanup in case fftw-threads were used pde::fftw_clean();returnEXIT_SUCCESS; }

$ g++ -Wall -O2 -c laplace.cpp $ g++ -Wall -O2 -c example.cpp $ g++ example.o laplace.o -lm -lfftw3 -lfftw3_threads -o example

#### Software integration

One application is in the area of High Dynamic Range (HDR) photography. The paper Gradient Domain High Dynamic Range Compression, 2002, by R. Fattal, D. Lischinski, M. Werman, describes a method which requires the solution of Poisson's equation. For this purpose, a slightly modified version of the source code on this page has been integrated into:- pfstmo
since version
1.5,
within the command line tool
`pfstmo_fattal02`

- Luminance HDR since version 2.3.0, within the tonemap operator Fattal

#### Benchmark

When compiled with the flag`-DTIME_REPORT`

the Poisson solver prints timing information. For example for
an array size of 5472x3648 the solver needs 16s using 1 thread
and 5s using 4 threads:
$ g++ -Wall -O2 -DTIME_REPORT laplace.cpp ...1 thread, 5472x3648, Intel Core i5 750, 2.67GHz, fftw library provided by linux distribution (not optimised):

poisolve(): 69 ms: rhs boundary conditions poisolve(): 8120 ms: rhs to EV space (inv fft) poisolve(): 324 ms: solve equation in EV space poisolve(): 7829 ms: solution to normal space (fft) poisolve(): 16341 ms: total time4 threads:

poisolve(): 68 ms: rhs boundary conditions poisolve(): 2421 ms: rhs to EV space (inv fft) poisolve(): 323 ms: solve equation in EV space poisolve(): 2070 ms: solution to normal space (fft) poisolve(): 4881 ms: total time

#### Other solvers

Fishpack which is a Fortran77 library. See also Fishpack at Netlib.#### Maths

To come ...Math symbols are rendered by MathJax which requires JavaScript.