Abstract:

Abstract

1
Introduction
The Algebraic Path Problem (APP) is a general framework which unifies sev- eral solution procedures for a number of problems into a single algorithmic formulation. Specific instances of the APP include problems from linear al- gebra (e.g. computing the inverse of a real non-singular matrix), from graph theory (e.g.  the transitive and reflexive closure, the all-pairs shortest paths problem, the minimum-cost spanning tree problem, etc.) as well as from lan- guage processing (e.g. generation of regular languages). These problems are among the most important computational problems in computer science and engineering area. The applications of the APP can be found in bioinformatics, network routing, control theory, and many others. As a result, the APP has been widely studied [1, 2, 3, 4, 5, 6]. The difference between many instances of the APP is in the underlying algebraic structure expressed as a closed semiring. There are many algorithms for solving the different instances of the APP. Among APP algorithms, a simple algorithm is well-known because this uni- fied algorithm is generalization of Warshall's algorithm for transitive closure, Floyd's algorithm for all-pairs shortest-paths problem [7], and the Gauss- Jordan method for matrix inversion [8]. The APP algorithm consists of triple- nested loop as in the dense matrix-matrix multiply-add (MMA) algorithm, except it has more strict data dependencies such that the outermost loop can- not be interchanged. To effectively use memory hierarchy in state-of-art processors, blocking algorithm is necessary to relax the required memory bandwidth and to obtain the high performance.  Block Floyd-Warshall (FW) algorithm and its data dependency have been discussed in [9, 10, 11, 7]. The blocking approach can be applied to solve other APP instances. The most compute intensive part of the block FW algorithm is the MMA in different semirings. This fact enables us to leverage a highly tuned MMA implementation for efficient computing of the APP. In this paper, we extend the previously designed block FW algorithm [7] to other instances of APP and show a parallel implementation of the unified block APP algorithm on the Cell Broadband Engine processor (Cell/B.E.). The Cell/B.E. is a heterogeneous multicore processor consisting of a general purpose processor, the Power Processor Element (PPE), and eight short-vector SIMD processors, so called the Synergistic Processor Elements (SPEs). As other contribution, this paper discusses how different optimization tech- niques affect the performance on the Cell/B.E. The techniques include com- puting with a block data layout, overlapping the computation with the data communication, parallelizing the block algorithm, and reducing the number of required registers in matrix "multiply-add" kernel. The remainder of this paper is organized as follows.  In Section 2, the overview of the APP and its applications are given, and then the scalar and block algorithms for the APP are introduced. In Section 3, the architecture of the Cell/B.E. is described in details. In Section 4, the discussion how the block algorithm for the APP is implemented on the Cell/B.E. is conducted. In Section 5, the results of performance evaluation are presented.  Finally, Section 6 describes conclusions and future work.
2
Algebraic Path Problem
2.1
Notations and Definitions
Let us consider a weighted graph G = (V, E, w) with a set of n vertices V = {1, 2, . . . , n } , a set of edges E ⊆ V × V , and a weight function w : E → S which assigns to each edge a weight from a special algebraic structure referred to as closed semiring ( S, ⊕, ⊗, *, Ō, Ī) where S is a set of elements; the "addition" ⊕ and the "multiplication" ✆ are binary operations ( S × S → S ); * is a unary operation ( S → S ) called closure ; and the zero Ō and the unity Ī are constants in S . The further discussion including the properties of the closed semiring can be found in [12, 3, 13, 14]. Let a path p in G be an arbitrary sequence of vertices p = ( i, k 1, k 2, . . . , k m, j ) which begins with i and ends with j . We then define a weight w( p ) of the path p as the product of all edges of the path:
Let P(i,j) denote the set of all paths from i to j . Then the algebraic path problem (APP) is to compute the "sum" d i,j of P( i,j) for all ( i,j ) pairs:
The APP can also be formulated in a matrix form [12, 3] by introducing a matrix closed semiring ( S n×n , ⊕, ⊗, *, Ō, Ī ) over the scalar closed semiring ( S, ⊕, ⊗, *, Ō, Ī), where S n×n is a set of n × n matrices; the matrix "addition" ⊕ and the matrix "multiplication" & are binary operations ( S n x n x S n x n → Sm×n); * is a unary operation ( Sm×n → Sm×n ) called closure of a matrix ; the zero matrix Ō the elements of which are all Ōs, and a unity matrix Ī with Īs on the main diagonal and Ōs otherwise, are constant matrices. We associate an initial matrix A = ( α i,j ) in S n×n with the weighted graph G, where
The matrix operations of "addition" and "multiplication" are defined as in linear algebra: if A = ( ai.j ) and B = ( bi.j ) are n × n matrices, then the “addition” and the “multiplication” are
When we define an n × n matrix D = ( d i j ) of the elements specified in (1), we get the matrix formulation for the APP in terms of the matrix A , as follows:
The matrix A* , which is the closure of matrix A , can be rewritten as follows:
Then the APP in the matrix form is to solve (5) for A*. Note that for any idempotent semiring¹, the infinite sum (4) converges to a finite sum, because An+i = An for i > 0.
2.2
Scalar algorithm
We denote the initial elements of a matrix A by α ( ) = α ij . The following recurrence equation (6) is used to update the elements of a matrix A(k) = (a(k)) for k = 1 to n do αk,k ← (αk,k ′)*; % “Black” element update for all 1 ≤ i ≤ n do α i k ← α i k − 1) → α k k ; % “Red” elements update end for for all 1 ≤ i ≤ n (i ≠ k) do for all 1 ≤ j ≤ n (j ≠ k) do end for end for for all 1 ≤ j ≤ n do a(k), & a(k-1), % “Blue” elements update end for end for Figure 1: Scalar algorithm for solving the algebraic path problem (1 ≤ i,j ≤ n on each iteration k = 1,2, . . . ,n).
We can cast the recurrence (6) into a scalar algorithm shown in Fig. 1. As can be seen from (6) and Fig. 1, the algorithm consists of four kinds of updates on each k -th iteration ( k = 1,2, . . . , n ). The algorithm updates the "black" element (pivot) first. Then it updates the "red" elements (pivot column), the "white" elements (non-pivot), and the "blue" elements (pivot row) in order. For any idempotent semiring, the closure is a k k k = (a k − 1 )* = 1 and the pivot is α k − 1 ) = I for any k ∈ {1,2, . . . ,n} . This fact indicates that a a (k − 1) ▧ a k − 1) = aik −1) → (aik −1) * (a): k −1)* ≡ a (k − 1) for i = k, and a (k) for j = k . Thus, for the idempotent semiring, the scalar algorithm in Fig. 1 can be simplified as it is shown in Fig. 2. It is clear that the total number of combined (⊕, ⊗) operations in the simplified algorithm is
for k = 1 to n do for all 1 ≤ i,j ≤ n (i ≠ k and j ≠ k) do α ( k ) ← α i ,k 1) ⊗ α ( k − 1) ⊕ α i ,k − 1) , end for end for Figure 2: Scalar algorithm for solving the APP in an idempotent semiring
Instances of the Algebraic Path Problem
2.3
We can define different instances of the algebraic path problem by specializing closed semiring ( S, ⊕, ⊗, *, ō, ī). Several instances are detailed in following: ● Transitive and reflexive closure : the weights α i ,j are taken from S ⊆ {0,1}.  ⊕ = V, ⊗ = ∧, α* = 1 for all a in S, Ō = 0, and Ī = 1.  The closure of a matrix A = ( a i ,j ) gives the transitive and reflexive closure. ● All-pairs shortest paths (APSP): the weights aij are taken from S ⊆ R+U {∞}, where R+ is the set of positive real numbers. ⊕ = min, ☏ = +,α* = 0 for all a in S , Ō = ∞, and Ī = 0. The closure of a matrix A = (ai,j) gives the all-pairs shortest paths. ● Critical paths (CRP), also called maximum cost paths: the weights aij are taken from S ⊆ R + U {+∞, −∞}, where R + is the set of positive real numbers.  ⊕ = max, ⊗ = +, α* = 0 for all a in S, Ō = −∞, and I = 0. The closure of a matrix A = ( a i j ) gives the maximum cost paths, or +∞ if there are paths of unbounded cost [16]. ● Maximum capacity paths (MCP), also called tunnel problem or network capacity problem [1]:  the weights a j are taken from S ⊆ R + U {∞}, where R + is the set of positive real numbers. ⊕ = max, ☏ = min, a* = ∞ for all a in S , Ō = 0, and Ī = ∞. The closure of a matrix A = ( a i j ) gives the maximum capacity paths. ● Maximum reliability paths (MRP): the weights ai.j are taken from S ⊆ {0,1} and can be considered as reliabilities of the information transport between the two connected vertices. ⊕ = max, ⊗ = ×, α* = 1 for all a in S , 0 = 0, and 1 = 1. The closure of a matrix A = ( i,j ) gives the most reliability paths. ● Minimum cost spanning tree (MST): the weights aij = aj,i are taken from S ⊆ R+ U {∞}, where R+ is the set of positive real numbers. ⊕ = min, ⊗ = max, a* = 0 for all a in S, Ō = ∞, and Ī = 0.  Let D = (d i,j ) is the closure of a matrix A = ( a i,j ), then its entries d i,j = α,j are the edges of the minimum spanning tree[17]. ● Inverse of a real non-singular matrix : the weights ai.j are taken from S ⊆ R. ⊕ and ⊗ are the conventional arithmetic + and × operations on R respectively, a* = 1/(1 − a) for all a in R with a̸= 1 ( a* is undefined for a = 1), 0 = 0, and I = 1. The closure of a matrix A = (ai,j) gives (I − A)−¹ (see [3],[8]).
2.4
Block algorithm
Block algorithms are widely known to be more efficient in many computing en- vironments than the corresponding scalar algorithms because block algorithms enable to relax the required memory bandwidth. Figure 3 shows a block al- gorithm for solving the algebraic path problem. In the block algorithm, the initial n × n matrix A is divided into an N × N matrix of b × b blocks where N = n/b and b  (2 ≤ b ≤ n/2) is the blocking factor.  For simplicity, but without loss of generality, the matrix size is assumed to be in multiples of the blocking factor b. On each K -th iteration ( K = 1, 2, . . . , N ) of the outermost loop, the block algorithm works as follows: it first updates the "black" block A κ,κ (line 2) by applying the closure operation to the b × b subproblem A K,K , which is solved by the scalar algorithm of Fig. 2, then it updates N — 1 "red" blocks A i,k (lines 3-5) and N − 1 "blue" blocks A K,J (lines 6-8), and, finally, it updates ( N − 1)² "white" blocks A 1 , J (lines 9-13). On each iteration, the "black" block update requires b(b-1)" (= N(B,®)  (b) in (7)) "multiply-add" operations; the "red" blocks update and the "blue" blocks update require b 3 ( n/b − 1) "multiply-add" operations each; and the "white" blocks update requires b 3 (n/b − 1)² "multiply-add" operations. Thus the total number of “multiply-add” operations is
1: for K = 1 to N do 2: A(K) ← (A(K-1)*; % “Black” block update 3: for all 1 ≤ I ≤ N (I ≠ K) do 4: A(K) ← A(K–1) ➤ A(K) ⊕ Ō; % “Red” blocks update end for 5: 6: for all 1 ≤ J ≤ N (J ≠ K) do 7: A(K), ← A(K), ® A(K-1) ⊕ Ō; % “Blue” blocks update 8: end for 9: for all 1 ≤ I ≤ N (I ≠ K) do 10: for all 1 ≤ J ≤ N (J ≠ K) do 11: A(K) ← A(K) → A(K) ⊕ A(K) ;% “White” blocks update 12: end for 13: end for 14: end for Figure 3: Block algorithm for solving the algebraic path problem The total number of block operations is Ntotal(n,b) = n/b(1 + 2(n/b − 1) + (n/b − 1)²) = (n/b)³. On the other hand, the number of block operations for updating "red" (R), “blue” (B), and “white” (W) blocks is
Then, the ratio is
The ratio shows that, when ( n/b ) is relatively large, the block algorithm spends the most computing time to update "red", "blue", and "white" blocks; see Fig. 4, where the problem size n is measured in multiples of b . The "red", "blue" and "white" blocks are updated using matrix-matrix "multiply-add" (MMA) in different algebraic semirings. As it was shown, the MMA operation is the most compute-intensive part (kernel) in the block algo- rithm. Hence the computational performance of algorithm strongly depends on the performance of MMA operation. Figure 4: Ratio of the number of operations for updating "red", "blue", and “white” blocks compared to the total number Note that here the block algorithm is used for a solution of the APP with idempotent semiring, nevertheless, it can also be used with several modifica- tions as a solution for non-idempotent semirings [8]. For instance, to compute the inverse of a matrix, we have to make the following modifications to the algorithm in Fig. 3: ●  replacing the statement 4 with A 1 , K ← − A 1 , K ⑧ A K , K in the "red" blocks update; ● permuting the order of updates from "black → red → blue → white" to "black → red → white → blue", and changing the superscripts of statements accordingly; ● using the scalar algorithm (Fig. 1), instead of the simplified algorithm (Fig. 2), for the "black" block update.
3
Cell Broadband Engine
The Cell Broadband Engine (Cell/B.E.) is a heterogeneous multicore proces- The Cell/B.E. consists of one Power Processor Element (PPE), eight SOL. Synergistic Processor Elements (SPEs), a Memory Interface Controller (MIC), a Cell Broadband Engine Interface (BEI), and an Element Interconnect Bus (EIB) [18, 19]. PPE is a 64-bit Power Architecture based general purpose processor. PPE usually acts as the control center of the Cell/B.E.; it distributes computational workloads among SPEs and coordinates an entire operation. SPEs are 128-bit RISC processors with a dual-issue pipelined four-way SIMD unit.  Each SPE contains a Synergistic Processor Unit (SPU) and a Memory Flow Controller (MFC) ². SPU does not have data caches, but instead, has its own 256 kB local memory called local store (LS). The register file of SPU contains 128 128-bit registers. Each SPE can access the host memory and the LS of other SPEs with direct memory access (DMA) transfers. The single- precision floating point peak performance of each SPU at 3.2 GHz clock speed is 25.6 Gflop/s because each SPU can compute eight single-precision floating point operations per clock with the four-way SIMD fused multiply-add (FMA) instruction [20]. EIB is the communication path between all the components of Cell/B.E. and is composed of four data rings, where the theoretical peak data bandwidth at 3.2 GHz is 204.8 GByte/s. The peak data bandwidth between the XDR DRAM memory (host memory) and EIB is, however, 25.6 GByte/s.
4
Porting the APP Algorithm to the Cell/B.E.
Parallelization
4.1
The block algorithm (Fig. 3) has been implemented on the Cell/B.E. so that the PPE updates the "black" block and then the SPEs update the "red", "blue" and "white" blocks. It is essential for high-speed computation consider- ing how to parallelize the algorithm for the heterogeneous multi-core processor. The block algorithm sees an initial n × n matrix A = ( a ij ) as the n/b × n/b matrix A = ( A 1,J ) consisting of b × b blocks, where n is the problem size and b is the block factor. Hence we can parallelize the algorithm in terms of b × b independent block as a unit of computational workloads. To correctly parallelize the block algorithm, the synchronization between all the processing elements at regular time-points is needed. The algorithm has data dependencies between colored blocks updates: the "red" and "blue" blocks depend on the "black" block; and the "white" blocks depend on both the "red" and "blue" blocks. Thus, we need to make at least two synchronization points: the first point is after the "black" block update; and the second point is before the "white" blocks update. To implement this barrier synchronization, we use a mailbox mechanism, which offers a 32-bit data passing between PPE and SPEs [18].  At either of the synchronization point, each SPE sends a notification mail to the PPE, and after the PPE receives all the notification mails, the PPE sends an acknowledgment mail to each SPE. The "black" block update for K -th iteration ( K = 2, 3, . . . , N ) cannot be started till the A K,K in "white" blocks is updated for the previous ( K − 1)-th iteration. This means that PPE can immediately start updating the "black" block when an SPE finishes updating the A K,K block for the previous iteration. For that reason, we have implemented the scheduling that the A K,K block is updated in the beginning of all "white" blocks updates.
4.2
Block size and data layout
In the block algorithm, a given n × n matrix is partitioned with b × b square blocks.  This raises the question of the optimal block size for given architecture. We use 64 × 64 as the block size for the implementation. The size of block for data in single precision is 16 kB, which is the maximum data size to transfer with one DMA instruction; also, each the 16 kB block can be fitted in the 256 kB local store of SPE. Moreover, if block data layout [21] is used, which is discussed below, a whole block can be transferred with one DMA transfer Notice that for other existing block matrix algorithms on the instruction. Cell/B.E., it is a common practice to use blocks of 64 × 64 elements [22, 23, 24, 25]. A two-dimensional array of matrix data in C programming language is usually stored in row-major data layout (RDL), where all row elements are
Figure 5:  (a) Row-major Data Layout and (b) Block Data Layout; the number in each layout indicates the order of memory assignment. stored continuously (see Fig. 5a). The block algorithm, however, deals with blocks of a matrix data; hence the implementation requires to process the data as the matrix of blocks. For such data access pattern to block data, it is often effective to use a special data layout called as block data layout (BDL) [26, 21] (see Fig. 5b).
4.3
Matrix-matrix multiply-add in different semirings
Differences between MMA and APP
4.3.1
As it was discussed above, the most compute intensive part of the block APP algorithm is in the "red", "blue", and "white" blocks updates. We can uti- lize an existing fast matrix-matrix multiply-add implementation for updating the blocks. The differences of the MMA implementation and the APP imple- mentation are types and the number of instructions to obtain corresponding "multiply-add" result. The MMA implementation only requires a single fused multiply-add (FMA) instruction while other APP implementations need more number of instructions. The number and the type of instructions needed de- pend on an instance of the APP; each instance needs different kinds of algebraic operations, which include (+, ×) operation for matrix inversion, (min, +) for shortest paths, (max, +) for critical paths, (max, min) for maximum capacity paths, (max, ×) for maximum reliability paths, and (min, max) for spanning tree.  In the SPU of the Cell/B.E., the standard addition and the multiplication operations are implemented by a single instruction ( fa and fm assembly in- struction [27], respectively), but the minimum or maximum operation require two instructions (combination of fcgt and selb for each). For example, the all-pairs shortest paths (APSP) problem requires three instructions to imple- ment the addition operation and the minimum operation. The fact indicates that the APSP implementation on Cell/B.E. will take three times as many clock cycles as the MMA implementation; in other words, we can expect that the SPU's theoretical peak performance for APSP would be 8.533 (=25.6/3) Gflop/s.  Also, the maximum capacity paths (MCP) problem requires four instructions to implement a minimum operation and a maximum operation. The fact indicates that the MCP implementation will take four times as many clock cycles as the MMA implementation; in other words, we can expect that the SPU's theoretical peak performance for MCP would be only 6.4 (=25.6/4) Gflop/s.
MMA implementations on the Cell/B.E.
4.3.2
Currently, there are three available MMA implementations on the Cell/B.E. [28, 24, 23]. Each implementation uses the block size of 64 × 64 and is deeply tuned by using four-way short-vector SIMD instructions as well as loop un- rolling and software pipelining. The first MMA implementation was reported by Chen et al. [28]; it is written in C programming language and is distributed with IBM Cell SDK as a demo program. The other two implementations are written in an assembly language: Hackenberg reported a deeply unrolled im- plementation [24]; and Kurzak et al. [23] reported an ultimate MMA imple- mentation for the SPU. We have evaluated the performances of all the three implementations on the experimental environment explained in Section 5.1.  Table 1 shows the summary of the different MMA implementations. The implementation by Kurzak et al. [23] is the best one from both aspects of the performance and register usage; hence, we have selected this MMA implementation to be utilized for our APP algorithm. In fact, we have modified the selected MMA implementation so that the register usage is reduced. We need it because less register usage is desirable
for the APP implementation which requires more additional registers to store intermediate results. The modified MMA implementation uses only 49 regis- ters although the obtained performance of 25.44 Gflop/s is almost same as in [23] (see Table 1). When we utilize the modified MMA algorithm for our APP implementation, it uses 65 registers; 16 additional registers are consumed to keep the intermediate results.
5
Performance Evaluation
5.1
Evaluation environment
The examination results presented here were measured on the Cell/B.E. pro- cessor of a PlayStation 3 (PS3). In the PS3, six of eight SPEs are only available because the other two are disabled to improve chip yields, and to run some op- erations on the background of PS3 Linux. We have installed Yellow Dog Linux 6.2 as the operating system and used Cell SDK 3.1 as the software develop- ment environment. Cell SDK 3.1 contains gcc version 4.1.1 based compilers (ppu-gcc and spu-gcc). As the optimization options for compilation, we used -03 for the PPE program and -Os for the SPE program. The input data has been generated by Erdős-Rényi random graph generator in GTgraph [29]. The performance in Gflop/s reported in this section is calculated by using the following formula:
Each minimum or maximum operation is counted as a single floating-point operation; in other words, any (⊕, ⊗)-operation is counted as two flops, which is a common practice [30, 31, 32].
5.2
Performance of a 64 × 64 matrix-matrix "multiply-
add” in different semirings on one SPE
This section shows the performance of a 64 × 64 matrix-matrix "multiply-add" in different semirings on one SPE. Recall that the MMA in linear algebra with (+, ×)-algebra delivers 25.44 Gflop/s (99.38% of the peak). The performance of a 64×64 MMA for all-pairs shortest paths (APSP), critical paths (CRP), and maximum reliability paths (MRP) problem delivers 8.502 Gflop/s each (33.21% of the peak). As it was estimated in Section 4.3.1, the sustained performance is close to one-third of the peak.  On the other hand, the performance of a 64×64 MMA for maximum capacity paths (MCP) and minimum spanning tree (MST) problem delivers 6.378 Gflop/s (24.91% of the peak). The sustained performance of one-fourth of the peak is also coincident with the previous estimation. For comparison, we have also measured the performance of a 64×64 matrix- matrix "multiply-add" on PPE. PPE achieves only 1.801 Gflop/s and 1.625 Gflop/s for APSP, CRP and MRP, and for MCP and MST, respectively. PPE has a similar SIMD instruction set to SPE, and the potential peak performance of PPE is the same as the SPE; the achieved performance is, however, several times less than on SPE. The low performance is considered to be due to the less number of registers; the available 32 128-bit vector registers of PPE are not enough to keep all block data elements in the register file.
5.3
Performance on PPE and SPEs
Figure 6 shows the performance plot of the block algorithm for the different APP instances on PPE and six SPEs; the experiments were carried out for the problem sizes in multiples of 64 (blocking factor). Three instances, all-pairs shortest paths (APSP), critical paths (CRP), and maximum reliability paths (MRP) problems, have the same time of solution on the Cell/B.E. The other two instances, maximum capacity paths (MCP) and minimum spanning tree (MST) problems, have also the same time of solution on the Cell/B.E. The performance for relatively small problem size on any number of SPEs is not high. This is due to the lower performance of PPE which updates the "black" block, in addition to the synchronization overhead between PPE and SPEs. On the other hand, when matrix size n is very large (in our case, n ≥ 6912), the performance degrades drastically due to the memory shortage that results in many page faults. In fact, the number of page faults increases to many hundred thousands for such large matrix size from almost zero for the smaller sizes³. 256 MB of the host memory of PS3 is known as a limiting factor for the high performance computing [33]. Table 2 shows the maximum performance using PPE and different number of SPEs. Comparing the performance of six SPEs with one SPE, the speedup ratio is around 5.98, i.e. we have almost linear speedup increasing. We have made additional experiments to see the effects of two kinds of optimization techniques, which are multi-buffering (MBUF) and computation with a block data layout (BDL) 4 . The multi-buffering technique enables us to
Problem size overlap the computation with the communication, so that the communication latency can be hidden. BDL was explained in Section 4.2. Table 3 shows the results of the experiments.  Both of the optimizations increase the performance.
5.4
Related Work
There has been a number of researches related to block algorithms for solving the APP. The researches has been done, in particular, to accelerate the per- formance for solving the all-pairs shortest paths (APSP) problem because the APSP is the most popular among all APP instances. Venkataraman et al. [9] presented a block FW algorithm.  Park et al. [10, 11] implemented a recursive version of the block FW algorithm. Han et al. [30] reported an implementation of block FW algorithm with auto-tuning technique. Bondhugula et al. [34, 35] proposed an FPGA-based method for a block FW algorithm.  Matsumoto and Sedukhin [7] presented a block algorithm for APSP on the Cell/B.E., and, independently, Vinjamuri and Prasanna [36] also published a research result for APSP with a blocking approach but a different scheduling strategy.  In addition, Tadonki [37] showed a ring pipelined algorithm for APSP on the Cell/B.E. The Graphical Processing Unit (GPU) computing is a popular solution to obtain high performance, and several APSP implementations exist for Nvidia GPUs [38, 39, 40, 41, 42, 43]. Buluç et al. implemented a recursive block FW algorithm [42]; Okuyama et al. compared the recursive block FW algorithm and their iterative block algorithm, and showed the implementation of the iter- ative algorithm runs around 4% faster than that of the recursive algorithm in the case of relatively small problem sizes of 256 ≤ n ≤ 1024. The performance of some implementations on GPUs is over 100 Gflop/s which is more than twice as high as our result on the Cell/B.E. It is because some GPUs contain hundreds of processing elements and, as the result, have higher computational power although the clock speed is lower than in Cell/B.E.. Compared with GPUs, a winning side of the Cell/B.E. is the programma- bility of the SPEs thanks to the unique embedded local store and the large register file. This programmability makes the performance of the computa- tional kernel (64 × 64 "MMA") very close to the peak performance whereas it is impossible to obtain such high efficiency on current GPUs.
6
Conclusion
In this paper, we have extended a block Floyd-Warshall algorithm to the in- stances of Algebraic Path Problem (APP) for idempotent semirings and have showed an efficient algorithm implementation on the Cell/B.E. The results of performance evaluation show that parallel implementation of the block APP algorithm, which is rich in matrix-matrix "multiply-add" operations, is well suited for the Cell/B.E., and the achieved performance is near the estimated peak. We have also shown the effectiveness of using the optimization techniques. The implementation of the block algorithm with block data layout is faster than that of row-major layout. The evaluation result also shows that overlapping computations and data transfers by utilizing multi-buffering technique is very effective. As future work, we will port and evaluate the high performance algorithms to multi-core CPUs and many-core GPUs, and cluster systems with a number of such processors. Alternately, our presented block APP algorithm computes only the distances for all-pairs of vertices; however, finding additionally the paths would also be important for some practical applications of the APP. Therefore, another future work includes a design of efficient methods to com- pute both the paths and the distances.