# Parallelization of distinct array write access from subkernels

I'm working on an implementation of a multivariate FFT, which is (or at least should be) highly parallelizable due to the row-column-algorithm. However, i can't figure out how to implement that. The serialized idea looks like the following

MFFTs[dims_,data_] := Block[{e,f,locdat=data},
If[Dimensions[dims] == {1}, Return[Fourier[locdat]];
, (*else: row column*)
Do[locdat[[e,All]] = Fourier[locdat[[e,All]]], {e,1,dims[[1]]}];
Do[locdat[[All,f]] = Fourier[locdat[[All,f]]];,{f,1,dims[[2]]}];
Return[locdat];
]
];


which would take on a sample of e.g. $2^{14}$ data points on a grid of $4\times2^{12}$ samples (quite anisotrope) about 0.089 seconds. Note that the operations in the first loop work on complete distinct sets of data. After that the second Do-Loop does also never write on an array entry twice. So that should be parallelizable quite far. However even for that short set of data the code

MFFTp[dims_, data_] := Block[{e, f, locdat = data},
If[Dimensions[dims] == {1},Return[Fourier[locdat]];
,
locdat = ParallelTable[Fourier[a],{a,locdat}];
Do[locdat[[All, f]] = Fourier[locdat[[All, f]]];, {f, 1, dims[[2]]}];
Return[locdat];
]
];


would need about 0.0549 seconds (running on 4 subkernels).

Do you know a way to do something similar with the second one (or in more dimensions that would be encapsulated Dos) ?

Any Usage of ParallelDo would require to SetSharedVariable the locdat, which kills any timing due to synchronization (about 4 seconds in this example using 4 subkernels). Due to the distinct write access it would be nice to be able to write in parallel onto an array in the Master kernel or something like that.

### Edit: Update

One little improvement - at least limited to the 2D-FFT is using twice the command locdat = Transpose[ParallelTable[Fourier[a],{a,locdat}]]; to perform the FFT on each row and column. That won't work for a 3DFFT though. And it doesn't improve the timing that much.

As @Ajasja pointed out the bottleneck is the collection of the results, so some kind of Parallelwrite would be really great.

-
Are you aware that Fourier already is multivariate? E.g., for $2^{24}$ values, n = 2^4; data = RandomReal[{0, 1}, {n, n, n, n, n, n}]; Timing[ Fourier[data];] takes 0.6 seconds on a single kernel; most of that is to claim RAM for the calculation. (It doesn't parallelize; if you try, it will take ten times longer.) –  whuber Mar 23 '12 at 19:40
Could the bottleneck be the transfer of data as here mathematica.stackexchange.com/questions/2886/… Also, in your example code it is probably possible to replace Do[locdat[[All, f]] = Fourier[locdat[[All, f]]];, {f, 1, dims[[2]]}]; by a ParallelMap or ParallelDo. –  Ajasja Mar 23 '12 at 21:36
@whuber i wasn't aware of that, that migth at least speed up the sequential algorithm, though my emphasize is on parallelization –  Ronny Mar 23 '12 at 23:33
@Ajasja yes, you're right, the bottleneck is the transfer of data. I replaced the second Do by your two suggestions and finally by a parallelTable on the Transpose[] data, which would not solve the problem for more than 2 dimensions. I haven't had an idea yet how to apply your advice to my problem; will try that tomorrow, though. –  Ronny Mar 23 '12 at 23:36
@Ajasja - well the number of MemberQ and the syncronization actually is the bottleneck but i can't adapt the hint you linked, because my Arrays are still unpacked. –  Ronny Mar 24 '12 at 10:10

Firstly, I want to point out not only that Fourier works for arbitrary n-dimensional arrays (as already mentioned by whuber), but also that it's already very efficiently parallelized using threading (via the Intel MKL). Therefore, attempts to parallelize it further are futile unless you intend to distribute the workload over a cluster--and then, even using algorithms that require no synchronization, the cost of communication must be considered and carefully minimized. As Amdahl's law makes clear, the performance one may gain through parallelization is very strongly constrained by the remaining serial portion of the workload, and communication is inherently difficult to parallelize within the "scatter/gather" paradigm offered by the Parallel package. It is possible to implement MPI-style message passing on top of MathLink, but even then, Mathematica and MathLink are not ideal for the purpose, and you will definitely need a cluster with a high-performance interconnect such as InfiniBand in order to make this reasonably scalable.

All this being said, there's no reason not to at least try to implement a row-column decomposition to see what can be done with it. Taking your MFFTs to start with, and assuming that the input array will always have quite a high aspect ratio as in your question, we can get much better performance by doing the column-wise transforms "by hand" with a tensor dot product rather than calling Fourier once for each column. This is because the overhead of the function call is significant, and for small $N$ there's not that much to be gained from an $O(N M \text{log} N)$ algorithm versus one that works in $O(N^2 M)$ time (although this still implies $\approx$ 3 times as many operations for $N = 4$, so it helps that the dot product is efficient and well parallelized). Let's write:

DFTMatrix[n_] := 1/Sqrt[n] Table[Exp[(2 Pi I (r - 1) (s - 1))/n], {r, 1, n}, {s, 1, n}];

MFFTs2[data_?ArrayQ] /; ArrayDepth[data] > 1 :=
Block[{local = data, len = Length[data], i},
Do[local[[i, All]] = Fourier[local[[i, All]]], {i, 1, len}];
DeveloperToPackedArray@DFTMatrix@N[len].local
];

MFFTs2[data_] := Fourier[data]; (* fallthrough for 1-d arrays *)


which is a lot faster (and still works for arrays of arbitrary dimensionality). For an array with $(N, M) = (4, 2^{20})$ , the AbsoluteTiming of MFFTs2 is 0.52 seconds versus 12.16 seconds for MFFTs, and (this surprised me!) MFFTs2 is even faster than Fourier for arrays of certain non-power-of-2 dimensions such as $4 \times (2^{20}+1)$ , where AbsoluteTimings are 0.78 seconds for MFFTs2 and 1.45 seconds for Fourier. Here I've kept your Do construct because it doesn't unpack, unlike the simpler Fourier /@ local, which unpacks down to level 1 (although isn't significantly slower as a result). It's worth remembering, though, that both approaches still copy their input (i.e., the FFT is not performed in-place).

Now, thinking about parallelizing it, the main consideration is avoiding unpacking since there's really no way to avoid the communication cost of sending each of the sub-arrays at level 1 to the subkernels. As a result, performance will always be poor, even on a cluster; parallel FFTs are ordinarily used when the data do not need to be sent in their entirety as they have already been distributed, and it is then desired to form the distributed FFT. In any case, doing the best we can, we need both withModifiedMemberQ from this answer, and to avoid the distribution of the whole of the array when only the parts at level 1 are actually required by each subkernel. Writing that down, we get:

MFFTp2[data_?ArrayQ] /; ArrayDepth[data] > 1 :=
Block[{local = data, i},
local = withModifiedMemberQ@ParallelTable[
Fourier[i], {i, local},
DistributedContexts -> None
];
DeveloperToPackedArray@DFTMatrix@N@Length[data].local
];

MFFTp2[data_] := Fourier[data];


which, for an array of $(N, M) = (4, 2^{21} + 1)$ , gives an AbsoluteTiming of 4.59 seconds (of which 0.81 seconds is spent on communication), versus Fourier's AbsoluteTiming of 3.04 seconds. So, at least communication costs aren't totally overwhelming, and perhaps on a cluster with a fast interconnect it might outperform Fourier for certain (large) inputs. However, performance is still mediocre at best, which conforms to our generally low expectations for this kind of approach.

-
Thanks for this comprehensive answer. I know about Ahmdahl Law and don't think, the gain in 2D might be that much, but it might be in higher dimensions. Thanks for the hint, that Fourier@locdat unpacks, but do you have an Idea how to perform that in-place? That's what I can't figure out in Mathematica. Your parallel version looks great, I think it's surprising that the dot product is that much optimized. I already read about the withModifiedMemberQ but couldn't adapt it. I'm quite eager to try your idea on higher dimensional data. –  Ronny Mar 27 '12 at 7:46
So i tried your approach and it works great for the posted question. Though, if the $N$ (from your remarks) gets too large (say it's bigger that $M/2$ but less than $M$), the DFT-Matrix is quite slow. For those cases, using the ParallelTable on the columns of the data and transposing the result, is much faster. With that distinction, the 2D case might reach its limits, when $N<M$ are arbitrary. How to perform the second idea on 3D or 4D data, will be one of my next things to think about. –  Ronny Mar 27 '12 at 9:10
As I mentioned in the answer, I doubt if you can beat Fourier in general, but you might be able to find something that works well for you by experimentation. In regard to how to perform the FFT in-place: this is not possible using any legitimate approach, although I have an idea how one might do it (in version 8 only) by using highly questionable undocumented techniques. If I have time to get this working, I will update this answer and leave another comment here. –  Oleksandr R. Mar 27 '12 at 10:32
Okay, I'm also not sure whether performing in-place gains that much speed, I'm not that far away from Fourier, that's enough, I think, combining the two ideas mentioned, where $N>\log M$ seems to be a good barrier to change form the dot product to a second parallelTable (where the second parallelTable does not give good times for the original setting). –  Ronny Mar 27 '12 at 10:43
Well, if you find something you're happy with, don't forget to either add an answer or update your question! –  Oleksandr R. Mar 27 '12 at 10:52

The answer of Oleksandr is quite nice and solves the stated initial question, if the dimensions of the array (let's assume it is transposed such that the Dimensions are in increasing order) are quite different, e.g. as stated $2\times 2^{14}$. Though, when increasing the first, there is a point, when the explicit construction of the Fourier Matrix is too much time comsuming. Let $N$ and $M$, $N<M$ denote these dimensions (see answer of Oleksandr), for me that was the case, when $N \geq M/4$, but I'd only check powers of $2$ and not that many different values of $M$. For these cases the following performs better, where i only state the parallel version, the serial should be obvious and I'm also using the withModifiedMemberQ mentioned here .

MFFTp3[data_?ArrayQ]/;(ArrayDepth[data]>1):=
Block[{local=data,i},
local = withModifiedMemberQ@ParallelTable[
Fourier[i], {i, local},
DistributedContexts -> None
];
local = Transpose[local];
local = withModifiedMemberQ@ParallelTable[
Fourier[i], {i, local},
DistributedContexts -> None
];
Transpose[local]
];
`

This is of course (just) the Row-Column Algorithm in a parallel manner inspired by the first step given by Oleksandr.

-