FFTW 3.3 MPI Cheat Sheet

Because I believe they provide a nice way to perform P3DFFT- and 2DECOMP&FFT-like MPI-parallel pencil decompositions, I have have spent some time staring at FFTW's 3.3 alpha and beta releases. In particular, poring over their Distributed-memory FFTW with MPI capabilities.

Using 3.3-alpha1 I wrote a small library (called underling) which mimicked P3DFFT 2.3's data movement capabilities. I wholly isolated the data movement from computing the FFTs. This, unlike P3DFFT and 2DECOMP, provides well-defined memory layout at any stage during the pencil transposes. Functionally my first approach was quite sound. However, it performs suboptimal on-node memory reshuffling (a design mistake on my part) and can stand to be improved.

I am revisiting the assumption of separating the parallel FFTs from the parallel data movement. Doing so allows using the higher-level 2D r2c/c2r planning APIs freshly documented for 3.3-beta.  It should require less needless memory reshuffling when combined with appropriate FFTW_MPI_TRANPOSED_IN and FFTW_MPI_TRANSPOSED_OUT flag choices.

Because, though wonderfully written, the relevant sections of the FFTW MPI documentation do not provide a cheat sheet, I have cooked my own for the various piece parts I may use in a second attempt. Perhaps someone will find it useful. Be sure to check the FFTW MPI reference for full details, especially the local_size calls necessary to obtain data distribution information.

Please tell me if you catch any mistakes. Many thanks to Steven G. Johnson for correcting my earlier misunderstanding of the real-to-complex and complex-to-real 2D DFT semantics.

							transposed in					transposed out					transposed in|out
transpose	in	n0/P	×	n1	×	nc	n1	×	n0/P	×	nc	n0/P	×	n1	×	nc	n1	×	n0/P	×	nc
	out	n1/P	×	n0	×	nc	n1/P	×	n0	×	nc	n0	×	n1/P	×	nc	n0	×	n1/P	×	nc
c2c 2D	in	ñ0/P	×	ñ1	×	ñc	ñ1/P	×	ñ0	×	ñc	ñ0/P	×	ñ1	×	ñc	ñ1/P	×	ñ0	×	ñc
	out	ñ0/P	×	ñ1	×	ñc	ñ0/P	×	ñ1	×	ñc	ñ1/P	×	ñ0	×	ñc	ñ1/P	×	ñ0	×	ñc
r2c 2D	in	n0/P	×	2(n1/2+1)	×	nc	2(n1/2+1)/P	×	n0	×	nc	n0/P	×	2(n1/2+1)	×	nc	2(n1/2+1)/P	×	n0	×	nc
	out	ñ0/P	×	(ñ1/2+1)	×	ñc	ñ0/P	×	(ñ1/2+1)	×	ñc	(ñ1/2+1)/P	×	ñ0	×	ñc	(ñ1/2+1)/P	×	ñ0	×	ñc
c2r 2D	in	ñ0/P	×	ñ1/2+1	×	ñc	(ñ1/2+1)/P	×	ñ0	×	ñc	ñ0/P	×	ñ1/2+1	×	ñc	(ñ1/2+1)/P	×	ñ0	×	ñc
	out	n0/P	×	2(n1/2+1)	×	nc	n0/P	×	2(n1/2+1)	×	nc	2(n1/2+1)/P	×	n0	×	nc	2(n1/2+1)/P	×	n0	×	nc

Notation: Real-valued directions are denoted n0 and n1 while complex-valued directions are ñ0 and ñ1. Half-complex storage is denoted ñ/2+1 and its padded, real-valued counterpart is 2(n/2+1). Directions decomposed along a communicator with P processes are denoted n/P. nc stands for "number of components" and corresponds to the advanced planning API's howmany arguments.

My divergence divergence.

I've spent three full work days chasing down a verification failure within my spectral, compressible Navier—Stokes code.  I can't remember the last time I was so happy to average 1/3rd of a source line per day.

Using Boost Spirit 2.1+ to evaluate constant arithmetic expressions

I used Boost Program_options for my research code's input parsing requirements. As expected, it works well. For floating point quantities, I've found having to specify 4.18879020478639 for 4π/3 is a bit cumbersome. When specifying arguments on the command line, I could use any one of a number of hacks to perform floating point arithmetic in bash. But that didn't help me when taking input from response files. I wanted an arithmetic expression utility that could process a C-like grammar. Long double accuracy was a requirement. I looked around a bit. GNU libmatheval and muParser appeared to be good fits but looked to be double precision-only. Though straightforward, I decided against combining the Shunting-yard algorithm with some in-place RPN logic because I wasn't enthusiastic about either finding or writing an appropriate (and admittedly toy simple) lexer.

Already having a Boost pre-requisite in my code, I decided to try out Boost Spirit since it could provide a header-only, precision-agnostic solution. Spirit has gone through a lot of significant revisions and many articles seemed to be written for the older Spirit Classic. The latest official reference manual gave several excellent examples that almost fit. I adapted Spirit's calc3.cpp and roman.cpp sample grammars to accomplish what I wanted. I added ** as an operator for exponentiation and pi as a built-in constant. I also added a smattering of special functions (e.g. exp). The task required about 150 lines of Spirit-based code. Further extensions should be fairly straightforward.
My sample code does the following:

1. Declares X-macros for the unary and binary functions built into the grammar.
2. Declares Boost Phoenix-ready functors for the unary and binary functions.
3. Declares expression::grammar consisting of the rules expression, term, factor, and primary.
4. Declares a helper called parse to hide needing boost::spirit::ascii::space when performing whitespace-agnostic parsing.
5. Declares a small test macro called EXPRTEST to reduce test case boilerplate.
6. Runs a collection of tests to ensure the parser is behaving as expected.

Compared to the samples on which I built, the biggest changes to the grammar were adding primary and reworking factor. I could have shortened up some details with a few using boost::spirit::qi; declarations, but I wanted to keep track of exactly which parts of Spirit I used. Feedback and suggestions are most welcome.

Updated 9 Jan 2012: During a second pass motivated by reducing compilation times, I removed the X-macro nastiness and replaced it with some qi::symbols-based logic. The code in the anonymous namespace could likely be removed if I could figure out the right boost::phoenix::bind hoodoo to use to replace lazy_ufunc_ and lazy_bfunc_ (any suggestions?). I have not yet measured the difference in either memory footprint or runtime performance. Compilation time and memory is greatly reduced for both GNU and Intel compilers.

Updated 31 Mar 2013: I have moved the code into a gist to facilitate both maintenance and gathering feedback. I have added the constant e and slightly reordered the primary grammar per the comments.

Updated 17-19 Sept 2013: I have added an MPL2 header to the gist. In the comments is a suggestion for how to use standalone MPL2 logic building atop this grammar from my project Suzerain to incorporate such parsing into your own, possibly non-MPL2, source tree.

Adding enthought.traits.api.HasTraits as an ancestor using Python metaclasses

Today, while banging my head against a wall trying to get a SWIG-generated Python module to play nice with Python Traits, I cooked a way to use Python metaclasses to modify a class to add enthought.traits.has_traits.HasTraits as a base class:

from enthought.traits.api import HasTraits
from enthought.traits.has_traits import MetaHasTraits

# See http://onlamp.com/lpt/a/3388 for an intro to Python metaclasses and
# http://code.activestate.com/recipes/204197-solving-the-metaclass-conflict/
# for what can go wrong here
class AddHasTraitsAncestor(MetaHasTraits):
    def __new__(cls, name, bases, dct):
        return MetaHasTraits.__new__(cls, name, cls.modify(bases), dct)

    def __init__(cls, name, bases, dct):
        super(AddHasTraitsAncestor, cls).__init__(name, cls.modify(bases), dct)

    @staticmethod
    def modify(bases):
        if (bases == (object,)):
            return (HasTraits,)
        else:
            new_bases = list(bases)
            new_bases.insert(0, HasTraits)
            return tuple(new_bases)

Adding the line __metaclass__ = AddHasTraitsAncestor to a class definition is enough, in some cases, for the class to now have HasTraits as a base class. Most of the frustrating bits were before I realized HasTraits has its own metaclass MetaHasTraits which AddHasTraitsAncestor had to subclass before the approach worked.

"Uhhh...?" you say. That puzzled look asking "Why didn't you just add the base class yourself?". I used this approach because SWIG was generating the Python class definitions and I was abusing SWIG's pythoncode feature to tweak the SWIG wrappers to be Traits-ready a la

%extend someCPlusPlusClass {

    %pythoncode %{
        __metaclass__ = AddHasTraitsAncestor
    %}

}

Turns out it works and my SWIG wrappers are now Traits-friendly. I'm still working on getting Traits to see the SWIG-based members, however. Suggestions much appreciated.