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}];
Developer`ToPackedArray@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
];
Developer`ToPackedArray@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.
Fourieralready 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.) $\endgroup$Doby your two suggestions and finally by aparallelTableon theTranspose[]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. $\endgroup$MemberQand the syncronization actually is the bottleneck but i can't adapt the hint you linked, because my Arrays are still unpacked. $\endgroup$