Direct solver for Poisson's equation in C++

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

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