We begin with the following code snippet :

#include "workspace/workspace.hpp"
#include "transform/fftw/fftw_r2c.hpp"
#include "transform/fftw/fftw_c2r.hpp"
#include "transform/fftw/fftw_local_array.hpp"
#include "coefficients/coefficients.hpp
#include "arrays/cartesian/local_geometry.hpp"
#include "arrays/array.hpp"
#include "coefficients/algebra.hpp"
#include "functors/x_functors.hpp"
#include "functors/k_functors.hpp"
#include <cmath> // for M_PI
int main(int argc, char** argv)
{
	using namespace spectral;
	workspace::initialize(argc,argv);
	
	// ...
	// Insert upcoming code here.
	// ...
	
	workspace::finalize();
	return 0;
}		
	

For now, this looks a lot like our first "Hello world" example, except that we have added a bunch of #include directives. Now we are ready for a couple of typedefs:

typedef array<cartesian::local_geometry<1> , double > Array_t;
typedef fftw_r2c<Array_t> FFT_t;
typedef fftw_c2r<Array_t> IFFT_t;
typedef FFT_t::Source_t Coefs_t;
typedef FFT_t::Target_t FCoefs_t;
	

Array_t is the engine that will store the data in memory. Here, we want to use a local 1D array of double precision numbers, stored contiguously in memory, and that is accessed using normal indexing. Note that "double" is the default data type for arrays, so we could have ommited it. The next two lines define the Fourier transform classes associated with this engine type. All transforms have Source_t and Target_t as member typedefs, so we may as well take advantage of that to define the types of coefficients we will use, so that we are 100% sure that we can apply our Fourier transforms on them.

Coefs_t is the source type of the forward Fourier transform, so it must represent the values of the function at grid points. FCoefs_t is the target type of the forward Fourier transform, so it must represent Fourier coefficients.

Now we have to enter a very technical discussion that you should skip if you are not familiar with how FFTW stores its data. In fact, FFT_t::Target_t is equivalent to the following typedef:

typedef coefficients<array<cartesian::local_geometry<1> , double >, fourier_coefs_tag> FCoefs_t;
	

When using Fourier transforms, one must pay attention to the fact that they are actually done in-place, so that the source and target are stored in the same memory location. That shouldn't refrain us from referring to them via two different types, depending on whether they have been Fourier transformed already or not. The "array" engine that appears in the FCoefs_t typedef above means that it owns the data in memory. On the contrary, the type Coefs_t is only a view of the data : before applying Fourier transform, the last column is indeed meaningless, it serves only as a padding region. It is the user's responsibility to ensure that the data actually present in memory is consistent with the way it is referred to.

Let us now define our nodes and dig some paths between them:

Coefs_t f("f");
FCoefs_t f_hat("f_hat");
FFT_t().dig(f,f_hat);
IFFT_t().dig(f_hat,f);
	

We have declared two nodes of respective types Coefs_t and FCoefs_t, and we have dug a Fourier transform from f to f_hat, and an inverse Fourier transform from f_hat to f. The C-strings passed to the constructors of f and f_hat are user friendly names that we want to give them. At this point, no memory allocation has been done, and no Fourier transform has been initialized yet.

What do I mean by digging a transform between two nodes, and what actually happens when the 3rd line of code is executed ? First, a temporary object of type FFT_t is created, using the default constructor. This object represents the mathematical Fourier transform operator, and it can therefore exist happily by itself, without specifying anything about its source or target. Here, right after being created, the transform is used as a tool to dig a path between f and f_hat. It is as if we had simply drawn on a piece of paper two circles named "f" and "f_hat", and one arrow connecting f to f_hat. The temporary object representing the transform is then destroyed, but the digging action is remembered and will allow us to apply the Fourier transform later on. Another way to see it is that we have been digging tunnels between caves, and putting up little signs to indicate what they will be used for, but we have not sent any material yet for the machine to process.

Once we have declared our bunch of transforms, we can proceed to "late initialization". In this context, "late" means that the initialization can occur even after everything has been declared, instead of being intertwined with the declarations. Thanks to this many things can occur behind the scenes in a safe way, and a lot of hassle is avoided.

Let us initialize f_hat. For that, we retrieve the geometry_type, we choose the domain size, and we call the node::initialize member function :

typedef engine_traits<FCoefs_t>::geometry_type geometry_type;
geometry_type::global_domain_type domain;
domain[0] = 512;
f_hat.initialize(geometry_type(domain));
	

Once this has been done, you can think of f_hat as a smart pointer to an array of 512 double precision numbers that has been allocated somewhere on the heap, with the right memory alignment so that it can be used by FFTW without performance loss.

To initialize the f node and the two Fourier transforms, we only have to add the following call:

f_hat.propagate();	
	

Internally, nodes and paths are stored in a graph structure (based on the Boost Graph Library), keeping track of how they were connected to each other and in which order. In the context of Kicksey-Winsey, I refer to these graph structures as bunches. Calling the member function node::propagate on f will traverse the bunch to initialize all the nodes that are connected to f, and all the transforms between them.

Now that everything is ready, we can start pouring data into the bunch. We want to define things in physical space, so we will use the node f. To easily assign non-trivial coefficients, we will use an algebraic expression involving a functor:

functors::_x<Coefs_t,0u>::type x;		
f = sin(2*M_PI*x)/2.0;
	

The first statement declares a variable x whose type is determined by the compile-time type selector functors::_x. The first template parameter, here Coefs_t, is the type of coefficients for which we want to generate a functor. The second parameter is the direction of the coordinate, 0 meaning that we retrieve the first coordinate (there is only one anyway since we are in 1D).

The second statement assigns the values sin(2*pi*x) to the coefficients f, with x varying uniformly between 0 and 1 as one goes through the array. Thanks to all the "expression templates" machinery, and provided that suitable optimization flags are switched on, this statement will be simply translated by the compiler to a loop which could look like :

double* ptr = f->raw();
for(size_t i = 0; i != 512; ++i)
{
	ptr[i] = sin(2*M_PI*i/512);
}
	

To apply the whole chain of transforms that we have defined above in one go, we can make the following call:

f.inject();
	

The bunch of transforms will be traversed exactly in the same way as was done when calling "propagate", except that now the transforms will be applied instead of initialized. In our case, an FFT will be applied from f to f_hat, followed by an inverse FFT from f_hat to f.

If we want to derivate the function with respect to x, we must insert a multiplication by -ikx before the inverse Fourier transform:

f.push();
f_hat *= -functors::_ik<FCoefs_t,0>::type()/512;
f.pull();
f -= 2*M_PI*cos(2*M_PI*x);
trace_coefs(f);
	

When calling push() on a node, all the transforms having this node as their source are applied, in the order in which they were dug. Here, only the FFT will be applied since it is the only transform pointing out of f. The second instruction does the multiplication by ikx. The division by 512 takes care of the normalisation that is absent from the inverse Fourier transform. pull() is the counterpart of push() in that it calls all the transform having the node as their target. This will here result in an inverse Fourier transform, hence reconstructing the derivative of f.

Finally, we substract the expected form of the derivative from f and print its elements to standard output. Their values should be close to zero if everything works.