Dense Matrix Algorithms

C HAPTER 8 Dense Matrix Algorithms 1 The communication time of the algorithm for all-to-all personalized communication on a hypercube with cut-through routing (Section 4.4 (page 167)) is Tall to all per s = ts ( p − 1) + tw m( p − 1). With the given parameters, the communication time is 22428. Thus the algorithm for cut-through routing is faster for the given parameters. For a different set of parameters, the algorithm for store-and-forward routing may be faster. 4 Let us start from the central relation that must be satisfied in order to maintain a fixed efficiency. = K To √ = K (ts p log p + tw n p log p n 2 √ 0 = n − n(K t w p log p) − K t W 2 n = n2 = Kt √ w p log p ± K 2 t 2w p(log p)2 Kt 2 K 2t 2 s p log w p(log 2 p p)2 + 4K (ts p log p) + K t s p log p 2 w (log p) + 4K (ts log p) p K 2t 2 2 ¿From the preceding equation, it follows that if K t w2 (log p)2 is greater than (which is most likely to be the case for practical values of K , tw , and p), then the overall isoefficiency function is ( p(log p)2 ) due to the tw term of To . √ √ 5 The communication time remains 2(ts log p + tw n 2 / p). The algorithm performs p matrix multiplications √ √ 2.81 with submatrices of size n/ p × n/ p. Each of these multiplications take (n / p 1.41 ) time. The total √ parallel run time is (n 2.81/ p 0.91) + 2(ts log p + tw n 2 / p). The processor-time product is (n 2.81 p 0.09) + √ 2( pts log p + tw n 2 p). The (n 2.81 p 0.09) term of the processor-time product is greater than the (n 2.81) sequential complexity of the algorithm. Hence, the parallel algorithm is not cost-optimal with respect to a serial implementation of Strassen’s algorithm ± w p log 6 This variant of the DNS algorithm works with n 2 q processors, where 1 < q < n. The algorithm is similar to that for one element per processor except that a logical processor array of q 3 superprocessors (instead of n 3 processors) is used, each superprocessor containing (n/q)2 hypercube processors. In the second step, multiplication of blocks of (n/q) × (n/q) elements instead of individual elements is performed. Assume that this multiplication of (n/q) × (n/q) blocks is performed by using Cannon’s algorithm for one element per processor. This step requires a communication time of 2(ts + tw )n/q . In the first stage of the algorithm, each data element is broadcast over q processors. In order to place the elements of matrix A in their respective positions, first A[ j, k] stored at processor P(0, j,k) is sent to processor P(k , j,k ) in log q steps. This data is then broadcast from processor P(k , j,k ) to processors P(k, j ,l) , 0 ≤ l < q, again in log q steps. By following a similar procedure, the elements of matrix B are transmitted to their respective processors. In the second stage, groups of (n/q)2 processors multiply blocks of (n/q) × (n /q ) elements each processor performing n/q computations and 2n/q communications. In the final step, the elements of matrix C are restored to their designated processors in log q steps. The total communication time is (ts + tw )(5 log q + 2n/q) resulting in the parallel run time given by the following equation: TP = n3 p n3 + (ts + tw )(5 log( 2 ) + 2 ) p p n The communication time of this variant of the DNS algorithm depends on the choice of the parallel matrix multiplication algorithm used to multiply the (n/q) × (n/q) submatrices. The communication time can be reduced if the simple algorithm is used instead of Cannon’s algorithm. 7 Let the division operation on the i th row (0 ≤ i < n) be performed during time step ti . Then in time step ti + 1, the active part of the i th row after division is sent to the processor storing the (i + 1)st row. In time step ti + 2, the processor storing the (i + 1)st row sends this data to the processor storing the (i + 2)nd row. In time step ti + 3, the elimination operation is performed on the (i + 1)st row. and in time step ti + 4, the division operation on the (i + 1)st row. Thus, ti +1 = ti + 4. In other words, four time steps separate the division operations on any two consecutive rows. Hence, in 4n steps, all the n division operations can be performed and Gaussian elimination can be completed on an n × n matrix. However, the communication operations of time steps tn−2 + 3, tn−1 + 1, tn−1 + 2, and tn−1 + 3 are not performed. Thus, the actual number of steps is 4(n − 1). 8 See papers by Geist and Romine [GR88] and Demmel et al. [DHvdV93] for solutions to Problems 8.8 (page 373), 8.9 (page 374), and 8.10 (page 374). √ 11 The communication time of each one-to-all broadcast is ( n). Therefore, the total time spent in performing √ these broadcast in all the iterations is (n n). Since the number of processors is n 2 , the communication √ time contributes (n 3 n) to the processor-time product, which is greater than the (n 3 ) complexity of serial Gaussian elimination. 12 13 As shown in Section 8.3.1 (page 353), disregarding the communication time, the parallel run time of Gaussian elimination with block-striped and block-checkerboard partitioning is n 3 / p and 2n 3 / p, respectively. Since the sequential run time is 2n 3 /3, the total overhead due to processor idling in the two cases is n 3 /3 and 4n 3 /3, respectively. With cyclic-striped mapping, the difference between the maximally loaded processor and the minimally loaded processor in any iteration is, at most, that of one row. Thus, the overhead due to idle time per iteration is O(np), and the overhead due to idle time over all the iterations s O(n 2 p). With cyclic-checkerboard mapping, the difference between the maximally loaded processor and the minimally loaded processor in any iteration is, at most, that of one row and one column. The maximum overhead due to √ idle time per iteration is O( p × (n 2 / p − (n/sqrt p − 1)2 )), which is O(n p). The overhead due to idle time √ 2 over all the iterations s O(n p). As shown in Section 8.3.1 (page 353), the parallel run time of the asynchronous (pipelined) version of Gaussian elimination with checkerboard partitioning is (n 3 / p). This implies that for all values of n and p , the order of magnitude of the communication time is either the same as or less than that of the time spent in useful computation. The reason is that the sequential complexity of Gaussian elimination is (n 3 ), and the parallel run time on p processors has to be at least (n 3 / p). Hence, the overall isoefficiency function is the same as 16 that due to concurrency. Since p ≤ n 2 , W = n 3 = ( p 3/2 ) is the isoefficiency function due to concurrency. The communication pattern (and hence, the parallel run time) of pipelined Cholesky factorization with checkerboard partitioning is very similar to that of pipelined Gaussian elimination without pivoting. Figure 8.1 illustrates the first 16 steps of Cholesky factorization on a 5 × 5 matrix on 25 processors. A comparison of Figures 8.11 (page 364) and 8.1 shows that the communication and and computation steps in the upper triangular halves of coefficient matrices are the same in both cases. For simplicity we remove the step of line 5 √ of Program 8.6 (page 375) and replace line 7 by A[k, j] := A[k, j ]/ A[k , k ]. The original A[k , k ] can be √ replaced by A[k , k ] later, whenever Pk,k is free, and this step is not shown in Figure 8.1. The actions in first eight steps of Figure 8.1 are as follows: (a) (b) A[0, 0] sent to P0,1 from P0,0 (iteration k = 0 starts). √ A[0, 1] := A[0, 1]/ A[0, 0]. (c) A[0, 0] sent to P0,2 from P0,1 ; (d) A[0, 1] (after division) sent to P1,1 from P0,1 . √ A[0, 2] := A[0, 2]/ A[0, 0]; A[1, 1] := A[1, 1] − A[0, 1] × A[0, 1]. (e) A[0, 0] sent to P0,3 from P0,2 ; A[0, 2] (after division) sent to P1,2 from P0,2 ; (f) (g) A[0, 1] sent to P1,2 from P1,1 . √ A[0, 3] := A[0, 3]/ A[0, 0]; A[1, 2] := A[1, 2] − A[0, 1] × A[0, 2]. A[1, 1] sent to P1,2 from P1,1 (iteration k = 1 starts); A[0, 0] sent to P0,4 from P0,3 ; A[0, 3] (after division) sent to P1,3 from P0,3 ; A[0, 1] sent to P1,3 from P1,2 ; (h) A[0, 2] sent to P2,2 from P1,2 . √ A[1, 2] := A[1, 2]/ A[1, 1]; √ A[0, 4] := A[0, 4]/ A[0, 0]; A[1, 3] := A[1, 3] − A[0, 1] × A[0, 3]; 17 A[2, 2] := A[2, 2] − A[0, 2] × A[0, 2]. Plots (a) and (b) in Figure 8.2 represent the standard and scaled speedup curves, respectively. 18 See Plot (c) in Figure 8.2. 19 Figure 8.3 plots the efficiency curves corresponding to the speedup curves in Figure 8.2. 20 Scaled speedup, in general, does not solve the problem of declining efficiency as the number of processors is increased. Only if the isoefficiency function of a parallel system is linear (that is, ( p)), then the efficiency corresponding to the scaled speedup curve will be constant. 21 The following table gives the largest n corresponding to a p such that two n × n matrices can be multiplied in time less than or equal to 512 units. However, it should be noted that these values are based on the expression for TP given by Equation 8.14 (page 347). Equation 8.14 (page 347) gives an expression for the parallel run time when the matrices are partitioned uniformly among the processors. For the values of n given in the following table, it may not be possible to partition the matrices uniformly. p n 1 8 4 13 16 22 64 35 256 58 1024 97 4096 167 It is not possible to solve an arbitrarily large problem in a fixed amount of time, even if an unlimited number of processors are available. For any parallel system with an isoefficiency function greater than ( p) (such as, the one in Problem 8.17 (page 375)), a plot between p and the size of the largest problem that can be solved in a given time using p processors will reach a maximum. See the solution to Problem 5.9 (page 231) for more details. 23 Steps 3 and 6 of Program 8.7 (page 376) involve the initialization of the resultant matrix and can be performed in constant time. Step 8 requires processor (i, j) to read the values of elements A[i, k] and B[k, j] for k = 0 to n − 1. A CREW PRAM allows these read operations to be performed concurrently. There are 3 memory read operations and one memory write operation and one multiply and add operation. Therefore the total time for step 8 is 4tlocal + tc . For n iterations of this loop, the total time taken is n(4tlocal + tc ). The formulation is cost optimal because the processor time product of (n 3 ) is equal to the serial runtime of the algorithm. 24 In absence of the concurrent read operation, memory reads in step 8 will be serialized. During iteration 1, (k = 0), a processor (i, j) accesses values C[i, j ], A[i, 0], and B[0, j]. There is no conflict in accessing C[i, j]. However, there are n processors trying to access a single element of A and n processors trying to access a single element of B. Similarly, in any iteration, there are n accesses to any memory location. Serializing these accesses, we have a single memory read and write operation (for C[i, j ]), n read operations for A[i, k ], and n read operations for B[k, j ]. Therefore, the time taken for this step is 2tlocal + 2ntlocal + tc . For n iterations of the loop, the total time taken is n(2tlocal + 2ntl ocal + tc ). The algorithm is not cost optimal because the processor-time product is (n 4 ), which is higher than the serial runtime of (n 3 ). 25 26 Let us consider the time taken in step 8 of Program 8.7 (page 376). Elements C[i, j ], A[i, j ], and B[i, j ] are available locally. The remaining values have to be fetched from other processors. We assume that at any point of time, only one processor can access a processors’ memory. In any iteration, there are n − 1 nonlocal accesses to each element of matrices A and B, and one local reference. This takes time 2(tlocal + (n − 1)tnonlocal ). Furthermore, 2tlocal time is taken for reading and writing the values of C[i, j]. Therefore, step 8 takes time tc + 2tlocal + 2(tlocal + (n − 1)tnonlocal ). For n iterations, the total time is n(tc + 2tlocal + 2(tlocal + (n − 1)tnonlocal )), which is equal to ntc + 4ntlocal + 2n(n − 1)tnonlocal . (a) Access times described in Problem 5.31: In Program 8.8 (page 377), step 8 staggers access to elements of matrices A and B such that there are no concurrent read operations. Therefore, step 8 can be performed in time 4tlocal + tc . For n iterations, the total time taken is n(4tlocal + tc ). The formulation is cost optimal because the processor time product of (n 3 ) is equal to the serial runtime of the algorithm. (b) 27 28 Access times described in Problem 5.32: During any iteration of step 8, one processor is accessing an element of A locally and the remaining processors are accessing n−1 different elements of A non-locally. This is also the case with matrix B. Therefore, the time taken for executing step 8 is tc +2tlocal +2tnonlocal . The formulation is still cost optimal because the processor time product of (n 3 ) is equal to the serial runtime of the algorithm. √ √ Consider the case in which p processors are organized in a logical mesh of p × p. Each processor is √ √ assigned a block of n/ p × n/ p elements. In the k th step, the processor responsible for computing C[i, j ] reads the values of A[i, k] and B[k, j] and keeps a running sum C[i, j] = C [i, j] + A[i, (i + j + k)modn] × B[(i + j + k)modn, j ]. The processor then computes the values of the other product matrix elements √ assigned to it. To compute each value of C , the processors read 2(n − n/ p) values from non-local memories √ and 2n/ p values locally. The time taken for computing each value of the resulting matrix is ntc + 2tlocal + √ √ 2tlocal n/ p + 2tnonlocal (n − n/ p). The parallel runtime of this formulation is n2 n n3 tnonlocal + √ × tlocal . TP =   p p n 2p n × tc + 2 n − √ p p first receives all the data it needs to compute its n 2 / p elements and then In this algorithm, each processor √ performs the computation. The total amount of non-local data required by a processor is 2n 2 / p − 2n 2 / p. The time taken to compute all n 2 / p elements assigned to a processor is TP = n3 p tc + 2tlocal + 2tlocal n2 n2 n2 + 2tnonlocal (2 √ − 2 ) p p p For large values of n and p, this expression becomes TP = n3 n2 tc + 2 √ tnonlocal . p p The relative performance of the two formulations of Problems 8.27 (page 377) and 8.28 (page 377) becomes obvious when we compute the speedups obtained from the two. The speedup from the previous formulation (Problem 5.34) is p S = . 1 + 2tnonlocal / tc ¿From the above expression, the algorithm’s maximum speedup is determined by the ratio 2tnonlocal /tc . Typically the time for a non-local memory access tnonlocal is much higher than tc ; consequently, the speedup may be poor. The seedup of the latter formulation is S = p √ 1 + p/n × 2tnonlocal / tc . As the size of the matrices increase for a fixed value of p, the contribution of essential computation grows faster than that of the overhead. It is thus possible to obtain parallel run times close to optimal by solving larger problem instances. Furthermore, the speedup approaches p when n/ p ≫ 2tnonlocal /tc . This algorithm is therefore better than the first. 29 Problems 8.23 (page 376)–8.28 (page 377) illustrate the inadequacy of PRAM models for algorithm design. They fail to model the communication costs of parallel algorithms on practical computers. Consequently, an algorithm designed the PRAM model may exhibit very poor performance, even on well connected practical architectures like shared address space computers. In general, for an architecture in which the cost of non-local data access is more expensive than local data access, maximizing locality is critical. All practical architectures fall into this category. Furthermore, these problems also illustrate that shared address space computers should in fact be programmed with a view to maximizing data locality, much the same way message passing computers are programmed. Ignoring data locality in a shared address space computer leads to poor performance. (0,0) (0,1) (0,2) (0,3) (0,4) (0,1) (0,2) (0,3) (0,4) 1 (0,1) (0,2) (0,3) (0,4) (1,0) (1,1) (1,2) (1,3) (1,4) (1,0) (1,1) (1,2) (1,3) (1,4) 1 (0,1) (0,2) (0,3) (0,4) (1,0) (1,1) (1,2) (1,3) (1,4) 0 (1,1) (1,2) (1,3) (1,4) (2,0) (2,1) (2,2) (2,3) (2,4) (2,0) (2,1) (2,2) (2,3) (2,4) (2,0) (2,1) (2,2) (2,3) (2,4) (2,0) (2,1) (2,2) (2,3) (2,4) (3,0) (3,1) (3,2) (3,3) (3,4) (3,0) (3,1) (3,2) (3,3) (3,4) (3,0) (3,1) (3,2) (3,3) (3,4) (3,0) (3,1) (3,2) (3,3) (3,4) (4,0) (4,1) (4,2) (4,3) (4,4) (4,0) (4,1) (4,2) (4,3) (4,4) (4,0) (4,1) (4,2) (4,3) (4,4) (4,0) (4,1) (4,2) (4,3) (4,4) (a) Iteration k = 0 starts 1 (b) (c) (d) 1 (0,1) (0,2) (0,3) (0,4) 1 (0,1) (0,2) (0,3) (0,4) 1 (0,1) (0,2) (0,3) (0,4) 1 0 (1,1) (1,2) (1,3) (1,4) 0 (1,1) (1,2) (1,3) (1,4) 0 (1,1) (1,2) (1,3) (1,4) 0 (2,0) (2,1) (2,2) (2,3) (2,4) 0 (2,1) (2,2) (2,3) (2,4) (2,1) (2,2) (2,3) (2,4) 0 (2,1) (2,2) (2,3) (2,4) 0 (3,1) (3,2) (3,3) (3,4) 0 (3,0) (3,1) (3,2) (3,3) (3,4) (3,0) (3,1) (3,2) (3,3) (3,4) (3,0) (3,1) (3,2) (3,3) (3,4) (4,0) (4,1) (4,2) (4,3) (4,4) (4,0) (4,1) (4,2) (4,3) (4,4) (4,0) (4,1) (4,2) (4,3) (4,4) (e) 1 (0,1) (0,2) (0,3) (0,4) (f) 1 (0,1) (0,2) (0,3) (0,4) 1 (4,0) (4,1) (4,2) (4,3) (4,4) (g) Iteration k = 1 starts (0,1) (0,2) (0,3) (0,4) 1 (0,1) (0,2) (0,3) (0,4) (1,2) (1,3) (1,4) (h) 1 (0,1) (0,2) (0,3) (0,4) (1,2) (1,3) (1,4) 0 1 (1,2) (1,3) (1,4) 0 1 (1,2) (1,3) (1,4) 0 1 (1,2) (1,3) (1,4) 0 (2,1) (2,2) (2,3) (2,4) 0 0 (2,2) (2,3) (2,4) 0 0 (2,2) (2,3) (2,4) 0 0 (2,2) (2,3) (2,4) 0 (3,1) (3,2) (3,3) (3,4) 0 (3,1) (3,2) (3,3) (3,4) 0 (3,1) (3,2) (3,3) (3,4) 0 0 (3,2) (3,3) (3,4) (4,0) (4,1) (4,2) (4,3) (4,4) 0 (4,1) (4,2) (4,3) (4,4) 0 (4,1) (4,2) (4,3) (4,4) 0 (4,1) (4,2) (4,3) (4,4) 1 (0,1) (0,2) (0,3) (0,4) 1 (0,1) (0,2) (0,3) (0,4) 1 (0,1) (0,2) (0,3) (0,4) 0 1 (i) 1 (0,1) (0,2) (0,3) (0,4) (j) 0 1 (1,2) (1,3) (1,4) 0 1 0 0 (2,2) (2,3) (2,4) 0 0 0 0 (3,2) (3,3) (3,4) 0 (4,1) (4,2) (4,3) (4,4) 0 0 (m) Iteration k = 2 starts (k) (1,2) (1,3) (1,4) 0 1 (2,3) (2,4) 0 0 0 (3,2) (3,3) (3,4) 0 0 (4,2) (4,3) (4,4) 0 1 (n) (l) (1,2) (1,3) (1,4) 0 1 (2,3) (2,4) 0 0 1 (2,3) (2,4) 0 (3,2) (3,3) (3,4) 0 0 0 (3,3) (3,4) 0 (4,2) (4,3) (4,4) 0 0 1 (o) (1,2) (1,3) (1,4) (4,2) (4,3) (4,4) (p) Iteration k = 0 ends Communication for k = 0 Computation for k = 0 Communication for k = 1 Computation for k = 1 Communication for k = 2 Computation for k = 2 Figure 8.1 Pipelined Cholesky factorization on a 5 × 5 matrix on 25 processors. 250 Plot (c) 200 150 Plot (b) 100 S 50 Plot (a) 0 0 50 100 150 200 250 300 p Figure 8.2 Standard, scaled, and isoefficient speedup curves for Problems 8.17 (page 375) and 8.18 (page 375). 0.9 Plot (c) 0.8 Plot (b) 0.7 0.6 0.5 0.4 E 0.3 0.2 0.1 0 Plot (a) 0 50 100 150 200 250 300 p Figure 8.3 Efficiencies corresponding to the speedups in Figure 8.2.