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

• 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.

The source code is released under the GPLv2 or above:

#### 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;

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);
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();

return EXIT_SUCCESS;
}


Compile:
 $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 time

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