研究简报与论文选登

P1Implementation of a Parallel Multi-dimensional Fast Fourier TransformationP20 COUPL+：并行PDE求解函数库P37 PETSc求解二维拉普拉斯方程P50 初值问题的微分-代数方程系统求解器P64 SMP集群系统上可扩展特征问题并行求解器的研究P73 基础并行软件平台建设与应用综述（续）P81 加卸载响应比理论（LURR）预测中国大陆地区2005年地震趋势取得新进展P82 2005年中科院超级计算用户成果交流会圆满结束征 稿 示启《超级计算通讯》是中国科学院计算机网络信息中心超级计算中心主办的综合性刊物，每季一期。

本刊兼具学术性、科普、要闻报道和趣味性，不定期开辟有“研究简报、论文选登”、“学习园地”、“网格论坛”、“超级计算机”、“专家论坛”、“信息港”、“交流之家”、“他山之石——文摘”、“本期人物”、“企业之星”、“会议信息”等栏目。

希望可以为广大高性能计算研究人员、应用开发人员、各级管理人员和高等院校师生提供一个相互学习交流的园地，一个获取信息的窗口，一个展示高性能计算及其应用成果的平台。

热烈欢迎所有致力于高性能计算研发的工作人员，以及社会各界人士为本刊撰稿。

来稿请注明作者姓名、工作单位、通讯地址、邮政编码、email地址。

稿件一经录用，酌付稿酬，并赠送本期《超级计算通讯》。

来稿请寄：100080北京市中关村南四街四号 北京349信箱中国科学院计算机网络信息中心超级计算中心《超级计算通讯》编辑部 徐 薇（收）或发送电子邮件至 wxu@研究简报与论文选登Implementation of a Parallel Multi-dimensional FastFourier TransformationX.B. Chi E. van den Berg Q. Cheng J. ChenComputer Network Information CenterChinese Academy of Sciences, Beijing, 100080, P.R. Chinachi@ ewout@ walls@Abstract: In this article we consider the parallel implementation of themultidimensional fast Fourier transformation and describe in detail a number ofaspects that need addressing; data flow planning, communication and data layout. Weintroduce three different data layouts and discuss their implications on the parallelFFT algorithm. We then consider the performance of the algorithms on a clustercomputer and evaluate the resulting speedup curves.1 IntroductionThe Fourier transformation is a widely known mathematical tool and has applications in a wide range of scientific areas such as signal processing, speech recognition and synthesis, geophysical research, and image processing, but also in less obvious fields such as the solution of certain partial differential equations, due to a number of useful properties of the transformation [1]. The 1-dimensional Fourier transformation and its inverse are defined as22()() (forward transformation)()() (inverse trasformation)kx kx F k f x e dx f x F k e dk ππ∞−−∞∞−−∞==∫∫ (1)When considering problems that require computational, rather than analytical mathematics, either because no explicit solution to the above functions can be found for the given f(x) or F(k), or because the functions are altogether unknown and we only know their evaluation at a number of regularly spaced values for x or k, or simply because it is computationally more efficient, we can turn to the discrete Fourier transformation (DFT) which is defined by210(),, with 0n i rk n n r k k y F x y x e r n π−===≤∑<(2)for 1-dimensional problems, with x m , the known sample values in the spatial domain—1—《超级计算通讯》 Vol.3,No.3,2005,9and y k , 0 ≤ k ≤ n - 1, the dual values in the frequency domain. By interchanging x k and y r , and k and r k at all other positions, and negating the exponent of e , we get the inverse transformation x = F n −1 (y ). The DFT defined above can easily be extended to an arbitrarily high d -dimensional DFT as follows;121122*********11,,...,,,...,000...... (with 0,1,...,-1, for 1)d d d d d d d n n n r k r k r k r r r n n n k k k i i k k k y x r n ωωω−−−=====∑∑∑i d ≤≤ (3)where we defined2i n n e πω= (4)for conciseness of notation. For the 1-dimensional DFT, we need to evaluate n summations, each consisting of n terms, so that the overall complexity of this operations seems to be O (n 2). However, in 1965 Cooley and Tukey [2] showed that by exploiting the inherent symmetry in the transformation, this complexity can be reduced to O (n log n ) using an algorithms commonly known at the fast Fourier transformation (FFT). For two and higher dimensional transformations similar symmetry exists and the concept underlying the FFT can be extended for the fast computation of these transforms. Perhaps more convenient is the decoupling of the dimensions reducing the problem to a series of 1-dimensional DFTs. This idea is best illustrated on the 2-dimensional DFT;121211221122121212121210201201111,,0(n n n n r k r k r k r k r r n n k k n n k k k k k k y x ωωωω===−−−−===∑∑∑∑,).x (5)For a fixed k 1, the term between brackets reduces to22212212201,,n r k r r n k k k y x ω=−=∑ (6)which is simply a 1-dimensional DFT, similar to the one defined in equation 2. Since we have n1 different values for k 1, and hence n 1 independent transformations, thecomputation of has a complexity of O (n 1 · n 2 log n 2). Using these intermediate values, we can complete the original DFT by computingy11112112101,,n r k r r n k k k y y ω=−=∑ (7)for all values of k 2. This amounts to a time complexity of O (n 2n 1 log n 1), leading to an overall complexity of O (n 1n 2 log n 1n 2). In terms of matrices, this means that we can first apply Fourier transforms on all rows, followed by transformations over all columns on the results of the first step. Decomposition of higher-dimensional transforms proceeds completely analogously, successively computing 1-dimensional—2—研究简报与论文选登 discrete Fourier transformations along each dimension.2 Sequential algorithmThe sequential algorithm of a multi-dimensional fast Fourier transformation could be as simple as subsequently performing one-dimensional transformations along each dimension. Direct implementation of this approach, however, will result in an algorithm with rather poor performance due to the hierarchical nature of computer memories in which high memory locality is preferred. To achieve this desired level of memory locality an additional matrix transposition step for two-dimensional transformations and dimension reordering steps for higher-order dimensions is needed, so that the input data for the Fourier transformation are stored in adjacent locations in memory. Such an operation, even at its cost of O(nm) for an n × m matrix, contributes significantly to the performance of the FFT as a whole. As such, the resulting sequence of operations for transformations becomes as is outlined in figure 1(a) for a two-dimensional, and as figures 1(c,d) for multi-dimensional transforms.The letters in the figures represent arrays that serve as the input and output of the different operations, with A denoting the array holding the initial data to be transformed, and B and C serving as intermediate arrays. Array B, in addition is designated to store the final transformation result. At any time only two different arrays are used which means that either B or C can coincide with A to save the amount of memory used.(a)(b)—3—《超级计算通讯》 Vol.3,No.3,2005,9(c)(d)Figure 1: Possible data-flows for (a) two-dimensional, and (c) even and (d) odd number higher- dimensional FFT transformations. In parallel execution a distributed transpose (b) replaces the common transpose operation.The arrows indicate operations on input data and point to the array to store the result of the operation, and hence represent the data flow of the algorithm. For a convenient discussion of this data flow, we introduce a parity value that is associated with each operation, with the exception of the first FFT. For outof-place operations, ie operations for which the input and output arrays differ, the parity value is set to 1, whereas for in-place operations, ie those where the input and output arrays coincide, the parity is set to 0.Since the initial and final arrays as well as all operations are known in advance, the only thing that remains to be determined is the exact data flow path. Although the given graphs show all possible connections, some paths will result in better performance than others and in certain cases some transitions may not exist and therefore cannot be used in the path.Fourier transformations can generally be done both in-place and out-of-place, depending on the implementation of the FFT routines. The FFT routines provided by FFTW 3.0.1 support both possibilities, and it was found that transformations with a length 2n were best done in-place, whereas all other transformations were faster out-of-place.—4—研究简报与论文选登The transposition operations are more demanding in that in-place transposition can only be done efficiently when the matrix is square, which also proved to be the best choice. For all other matrices we are required to use out-of-place transposition, reducing the number of possible paths.If all preferred transitions were taken the situation might well arise that the final output end up in array C, instead of B. This problem can be overcome in two ways; either add an addition step in which the data are simply copied to B or change the destination array of one of the oper- ations. A number of time measurement clearly showed that the second solution was superior to the first one, giving rise to the question which operation to change the parity of. Since transpositions can only be changed when the matrix is square, it is best to select FFT operations for this purpose.The use of transposition for reorganizing the data fails to work for the third dimension and above and a complete reordering is needed as illustrated in figures 1(c,d). Reordering is done by sequentially reading through the data and writing it with strides and offsets such that the desired dimension order is obtained. By grouping dimensions in pairs we can limit the use of reorder operations to roughly half of the higher-order dimensions. As an example, consider an n1 × n2 × n3 × . . . × n d array. After the Fourier transformation has been applied on all n1 × n2 sub-arrays, we reorder the array such that all dimensions along which no Fourier transformation has been applied are moved to the front, ie to locations with smaller strides between elements. In the first step, for an array with d = 5 dimensions, this would leave us an array of size n3×n4×n5×n1×n2. Then, we can compute the FFTs for n3, transpose the first two dimensions and compute the FFTs over the fourth dimension, resulting in an array of size n4×n3×n5×n1×n2. We continue with another reordering step, shifting the frontmost two dimensions towards the back of the dimensions that have not yet been transformed and shift these latter dimensions to the front; n5×n4×n3×n1×n2. In this case only the fifth dimension remains so we simply compute the FFT along this dimension. As a final step the dimension order of the array is restored by special reordering operation that reverses all but the original first two dimensions. This is done to allow the transformation to return data in transposed state (with the first two dimensions interchanged) which saves time, especially when an inverse Fourier transformation is performed at a later stage. Transposed data results from omitting the rightmost transposition step in figure 1(a). By deciding this final order at this stage of the algorithm allows for a strict separation of the transformation of the first two—5—《超级计算通讯》 Vol.3,No.3,2005,9dimensions and the remaining higher-order dimensions, as shown in parts (c) and (d) of the given figure.When considering higher-order dimensions for data planning, we only need to take into account the reorder operations which are always done out-of-place. The overall planning of the data flow for the sequential algorithm is done in reverse, starting at output array B and going back towards the output array of the first FFT, granting the preferred parity values for all operations. This means that in some cases the parity desired by the first FFT operation may not be granted. However, we had already determined that the best operation to act as a spill would be the FFT, such that we obtain the desired result.A final, optional operation not shown in figure 1 is a scaling operation that multiplies all data elements by a fixed constant. This operation is always done in-place and therefore does not affect the planning phase.3 Parallel algorithmIn the parallel algorithm, the workload needs to be distributed among the available processors. This is done by partitioning the data along the first dimension ( the number of rows ) into blocks of roughly equivalent size. Each partition is then assigned to a processor and the Fourier transformation can commence. In the first step a one-dimensional FFT is applied on the vectors covering all columns, for each rows that is part of the local data. The second step again comprises a matrix transposition in which the rows and columns are interchanged and array is partitioned and distributed along the now-transposed columns. In order to transpose a column, we first need to gather its elements that are distributed across the processors and because each processor does this, an all-to-all communication operation is required. As a result, the transposition operations as depicted in figure 1(a) can no longer be used and needs to be replaced by the steps given in (b) which consists of a communication step followed by a local transposition step. An example of this process can be seen in figure 2 showing the initial data distribution on the left, the intermediate distribution, after data transfer in the middle and the final result of the global transposition operation at the right side. In the third step of the parallel algorithm Fourier transformation is performed on the transposed columns with a vector length equal to the original number of rows. The fourth, optional step is similar to the second step—6—研究简报与论文选登transposing the data back to the initial form and can be omitted when the results are to be given in transposed form.Communication in the global transposition phase can be done in blocking and non-blocking mode. In blocking mode, communication primitive functions do not return until after all communication is completed, whereas in nonblocking mode, the functions return much faster requiring additional function calls at a later stage to check the result of these operations. The advantage of non-blocking communication is that some computations can be done during data transfer thus allowing partial overlap of the two and reducing the total execution time of the algorithm.In the Fourier transformation phase a potentially large number of independent transformations are performed. After one or more FFTs are completed, transfer of the results can be initiated and another set of transformations computed. After that the transfer can be finalized and the newly available results transmitted, until no more transformations need to be computed. This way nonblocking communication can be effectively used to reduce the communication overhead that puts an upper bound to the possible speed-up.Once processing of the first two dimensions has been completed, the higherorder dimensions can be transformed. At this point the only difference between the sequential and parallel algorithms is the extend of the first or second dimension, depending on whether the results are given in normal or transposed form. Therefore the transformation along these dimensions is done as shown in figures 1(c,d) without the need for any modifications or additional steps at the end.Planning of the data flow on the other hand needs to be changed in order to make the parallel algorithm perform well and correctly. When using the planning method given for the sequential algorithm in conjunction with nonblocking communication primitives problems may arise. This is due to the lack of arrays available during the transformation. Suppose what happens when an out- of-place Fourier transformation is combined with the non-blocking communication step. First, an FFT operation is applied to a number of vectors in, say, array B and the results are stored in C. These results are then transferred to other processors and, because communication is out-of-place, incoming data must be stored in B. This means that the data that have just been transformed are overwritten by partially transposed data. As long as the amounts of data send and received are equal for all processors this will work, but as soon as these amounts differ there will be at least one processor for which the data—7—《超级计算通讯》 Vol.3,No.3,2005,9that have yet to be transformed are overwritten by other data thus corruption the whole computation.There are two ways to avoid this problem; the first one is to finish all transformations prior to the communication step at the cost of the advantage that could be gained by overlapping the two steps. The second solution is to compute the Fourier transformations that are followed by commu- nication in-place whenever possible. Since communication is a more time consuming operation than transformation, this second solution is preferred over the first one which should be used only in case an in-place transformation is absolutely impossible.The above scheme can be realized by planning according to a number of rules. First, communication that is followed by non-blocking communication is strongly preferred in-place if the dimensions of the data require this. Second, the preferred parity of the hindmost, rather than that the first, FFT operation is discarded to allow non-blocking communication as much as possible. The assignment of arrays to operations is still done from the end back to the start of the total transformation, which means that a pre-planning phase is needed to detect possible problems and set a flag to modify the parity of the first FFT operation encountered. If non-blocking commu- nication occurs when the flag is still set this must also be changed to retain correctness of the parallel transform. In case blocking communication is used we merely need to check if all prefer- red parities can be granted, otherwise change that of the last FFT using the above method.4 Data structureIn order for the parallel algorithm to perform well, the communication overhead should be as small as possible. Communication time in part depends on the data layout, ie how the rows are distributed over the processors. In the following we will consider three different data layouts and discuss their advantages, disadvantages and the effects on communication.4.1 Blocked layoutThe first data layout, illustrated in figure 2(a), partitions the data into regularly shaped blocks. The partitioning along the rows is a physical one in that data in—8—(a)(b)(c)Figure 2: Data layout for the regular block-padded structure.different partitions reside in the memory of different processors, the column-wise partitioning is logical and reflects the assign- ment of rows after a global transposition operation. Higher dimensional data can be imagined as replications of the matrix shown in the above figure, such that a three-dimensional block of data (the third dimension holding the aggregated higher-order dimensions) is formed. The partitioning of this data remains as indicated by the dashed lines. The size of the blocks formed by the lines is rows np ⎡⎤⎢⎥ by columns np ⎡⎤⎢⎥, where np is the number ofprocessors. Padding is required if the number of rows or columns does not divide by np and is appended after the row or column with the highest index, ie right and bottom in the figure. The computational load distribution arising from this partition scheme is slightly skewed towards the lower numbered processors that hold full blocks of data. This, however, does not significantly affect the performance of the algorithm as the maximum number of rows processed by any processor, a value that dictates potential parallel performance, is minimal.(a)(b)(c)Figure 3: Data layout for the irregular block structure.Communication For the communication phase we interpret the data as a two-dimensional array, where the number of rows, including padding, is multiplied by the product of all higher-order dimensions, as illustrated in figure 4(a), which agrees with a row-major memory layout. Because the dimensions of the data operated on by each processor are the same, no problems arise from using non-blocking communication primitives allowing maximum overlap of computation and communication. For the best performance, communication is divided into a number of smaller operations each with a fixed upper bound on the number of rows transferred at one time. The given number of rows is determined from the optimal message size for a given machine, the total number of aggregated rows, and the number of columns in each partition.4.2 Strided layoutThe blocked layout described above has the problem that padding is explicitly incorporated in the input and output data. This means that either the user needs to take this artefact into account, or that an additional reordering step needs to be performed prior to and following Fourier transformation possibly requiring an additional block of memory. Since neither of these two solutions is entirely satisfactory, a new data layout is needed. The strided layout, depicted in figure 3(a), partitions the data such that no padding is required in the input and output arrays. The rows and columns are partitioned similarly which is a necessary requirement for the transposition phase or an inverse Fourier transformation operating on transposed data. The only difference is that the division along the rows is physical, ie these partitions are distributed amongthe processors, while the columns partitions exist only logically and do not come intoexistence until the communication phase. Given np processors, a set of x rows (or columns) is partitioned by successively assigning a number of them to a processor 0 ≤ p < np . This number is determined as follows; each processor is initially assigned x np ⎢⎥⎣⎦ rows while one additional row is added to all processors p < x mod np , so that the remainder of the above division is accounted for. Higher dimensions are dealt with in the same way as in the blocked layout.Communication Unlike the blocked layout where data communication proceeds fairly naturally with fixed block sizes, the strided layout introduces a number of issues that complicate the communication phase. These all arise from the variations in block sizes sent and received. The most obvious complication is the all-to-all communication itself needing to deal with different block sizes and offsets. This problem is amplified when dealing with higher- dimensional transformations, as shown in figure 4(b) where the blocks received from processor p = 3 have fewer rows than the others. These gaps need to be padded before higher- dimensional instances of the two-dimensional array can be stored. Consequentially, due to the lack of padding in the data sent, the higher-order dimensions can no longer be merged with the first dimension (number of rows), unless this dimension divides by the number of processors. This too implies that for each instance of the two-dimensional array, a separate communication phase is needed which may seriously affect parallel performance. A final consequence of the uneven block size is that in-place transpositions ( with source and destination arrays equal ) are possible only if the number of rows and columns are equal and divide by the number of processors used.(a)(b)(c)Figure 4: Data layout multi-dimensional arrays (a) block padded; (b) strided; (c) continuous. 4.3 Contiguous layoutAlthough the strided layout eliminated the need for explicit padding it did introduce intermediate padding in the communication phase and as a result, problems with communication of higher-order dimensional data. The contiguous data layout was conceived to overcome these problems while maintaining a layout that does not require padding at the input and output level. The contiguous and strided data layouts are similar and differ only in the way received data are stored prior to the transposition phase. Thus, the assignment of rows and columns to processors remains unaltered from the strided layout. In figure 4(c) the data layout after receiving multi- dimensional data is given. Unlike the strided layout where data are stored in abutted blocks, the contiguous layout groups data according to the processor it originates from. Due of the lack of gaps in data, higher-order dimensions can once collated with the first matrix dimension allowing for more efficient communication. Another potentialadvantage of the contiguous layout is that the data need not be reordered into blocks that are subsequently transposed to yet other blocks which may have a negative effecton cache efficiency. Instead, the data is left in its sequential form and read likewisefor transposition.5 Performance evaluationEach of the three different data layouts were incorporated in the Parallel Multi-Dimensional Fast Fourier Transformation (pmdfft) library for testing. Writtenin C, the library benefits from the routines defined by the MPI standard and uses the one-dimensional Fouriertransformation routines provided by the FFTW library [3]. The library is freely available on-line and can be downloaded at . Similar to FFTW, pmdfft has a planning and an execution phase in which the transformations is prepared respectively computed. The planning phase includes the preparation of all one-dimensional Fourier transformations, the definition of all required MPI data types and the designation of arrays to operations. The execution phase takes the information that is stored in the plan and performs the actual transformation, communication and transposition operations.For the parallel performance evaluation, only the execution times are taken into consideration despite the fact that, depending on the dimension extends of the data,the planning phase can occasionally take a significant amount of time. However, when a great number of transformations are done, this overhead becomes negligiblefor individual transformations and can be safely omitted. As a reference, the planning and execution times of the sequential FFT are given in table 1 for a number of different array sizes.Table 1: Sequential planning and execution times for different matrix sizes.Normal TransposedExecution Planning Execution Size Planning1024×1024 0.010728 0.267979 0.009331 0.223688 2571×438 0.187368 2.511004 0.147429 2.418615 31×72×32×20 0.014203 1.070320 0.013491 0.987648 240×490×1000 0.008226 58.500558 0.008387 54.087522 The platform used for testing is the Deepcomp 6800 installed at theSupercomputing Center of the Computer Network Information Center – Chinese Academy of Sciences (SCCAS) in Beijing. The machine is an SMP cluster with a total of 256 nodes, equipped with four Itanium II processors running at 1.3 GHz and interconnected by a QsNet network, good for a peak performance of 5.324 GFlops.(a) 1024 × 1024, contiguous(b) 2571 × 438, contiguousFigure 5: Speed-up plotted against number of processors for parallel Fourier transformation of arrays of size 1024 × 1024 (a), and 2571 × 438 (b) for the contiguous data layout.Figures 5 and 6 give the speed-up curves for the arrays given in table 1. In figure 5 the results for the two-dimensional arrays 1024×1024 and 2571×438 with contiguous layout are given. The other two layouts are omitted here because of the similarity of their results. From the table we can see that although the amounts of dataof the two arrays differ only by approximately 7%, the executionand planning times differ considerable. This difference is caused by the favourable dimensions of the former matrix (powers of two) which are particularly suited for FFT, whereas the dimensions of the latter matrix (containing large prime factors) require different FFT algorithms that are much slower. This also results in a better parallel performance for the 2571×438 matrix, compared to the 1024×1024 matrix, because the relative communication overhead for the former is less and therefore has little impact on the speed-up of the algorithm; computation is the dominant factor. Note also that the non-blocking versions of the algorithm are generally faster than the blocking versions and as a result have better speed-up values. Leaving the output in tansposed form is not only faster in the sequential program but also leads to higher speed-up values because of the reduced communication overhead.In figure 6 the differences in performance between the three data layouts are more apparent. Most striking is the poor performance when using the strided layout. This is caused mainly because of its general inability to group together the higher-order dimensions thus prohibiting more efficient communication. The blocked and contiguous layouts show much better results both for the small 4-dimensional array, as well as for the large 3-dimensional array that requires much more communication. Also apparent is the difference in blocking and nonblocking communication modes, the latter typically being faster. In accordance with expectations, leaving output in transposed order gives better speed-up figures as the result of the reduction in communication.(a)31 × 72 × 32 × 20, blocked。