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

About

Download Source Code

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;

   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:

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.