COMP4300/6430 2013

Assignment 1 sample solution

Setup (not marked)

The naive implementation is significantly slower than MKL’s DGEMM implementation. MKL is highly optimised to use vector instructions and blocking for registers and cache.

The number of FLOPs required for matrix multiplication (both naive and DGEMM) is approximately 2mnk. Therefore we calculate GFLOP/s as (2mnk /(1.0e9 × time)). Peak FLOP/s for a single core of XE is 12.0 GFLOP/s (3.0GHz × 4 FLOPs per cycle). DGEMM gets quite close to this for matrices larger than about 200x200 - see figure 1. In comparison, the naive implementation achieves less than 1/8 of peak FLOPS for small matrices, and performs even worse for large matrices due to cache effects.

Figure 1: Performance of Naive vs. Intel MKL DGEMM on XE

Note on MPI error handling
The default error handler for MPI programs is MPI_ERRORS_ARE_FATAL, which means that an MPI error will abort the entire program. This means that it is allowable for user programs to call MPI functions without checking for errors after each call. See MPI 2.2 Error Handling.

Section 1: Replicated Data

1. See mult_replicated and mult_replicated2 in attached ass1.c.
   (Note: these are only two possible solutions - there are many other possibilities.)

2. mult_replicated divides p processes into a P×Q grid. Process 0 sends a row panel $\frac{m}{P}\times k$ of A and a column panel $k\times \frac{n}{Q}$ of B to each process in turn. Each process calculates the corresponding $\frac{m}{P}\times \frac{n}{Q}$ block of matrix C. Each process then sends its block of C to process 0.
   Assuming square matrices and a square grid of processes, the execution times for each phase of this algorithm are approximately:

• Communication 1 (distribute A and B): $PQ\left(t_{\mathrm{stup}}+t_{\mathrm{data}}\frac{mk}{P}\right)+PQ\left(t_{\mathrm{stup}}+t_{\mathrm{data}}\frac{kn}{Q}\right)=O\left(p+\sqrt{p}{N}^2\right)$
• Computation: $\frac{2mkn}{PQ}=O\left(N^3∕p\right)$
• Communication 2 (gather C): $PQ\left(t_{\mathrm{stup}}+t_{\mathrm{data}}\frac{mn}{PQ}\right)=O\left(p+N^2\right)$
• Total: $PQ\left(3t_{\mathrm{stup}}+t_{\mathrm{data}}\left[\frac{mk}{P}+\frac{kn}{Q}+\frac{mn}{PQ}\right]\right)+\frac{2mkn}{PQ}=O\left(p+\sqrt{p}{N}^2+N^3∕p\right)$

Computation and communication are not overlapped, so the total time is the sum of both. Given that, on XE, $t_{\mathrm{stup}}>t_{\mathrm{data}}\gg t_{\mathrm{flop}}$, we would expect communication to dominate for small matrix sizes. Communication is $O\left(p+\sqrt{p}{N}^2\right)$ and computation is $O\left(N^3∕p\right)$ , so the algorithm is computation-limited for large matrix sizes. As communication time increases with both the problem size and the square root of the number of processes, weak scaling (increasing the problem size with the number of processes so that each process has a constant amount of computation) will be increasingly dominated by communication.

mult_replicated2 broadcasts A to all p processes (using MPI_Bcast) and sends a column panel $k\times \frac{n}{p}$ of B to each process in turn. Each process calculates the corresponding $m\times \frac{n}{p}$ column-panel of matrix C. Each process then sends its column-panel of C to process 0.
The execution times for each phase of this algorithm are approximately:

• Communication 1 (broadcast A and distribute B): $\mathrm{log(p)}\left(t_{\mathrm{stup}}+t_{\mathrm{data}}mk\right)+p\left(t_{\mathrm{stup}}+t_{\mathrm{data}}\frac{kn}{p}\right)=O\left(p+\mathrm{log(p)}{N}^2\right)$
• Computation: $\frac{2mkn}{p}=O\left(N^3∕p\right)$
• Communication 2 (gather C): $p\left(t_{\mathrm{stup}}+t_{\mathrm{data}}\frac{mn}{p}\right)=O\left(p+N^2\right)$
• Total: $\mathrm{log(p)}\left(t_{\mathrm{stup}}+t_{\mathrm{data}}mk\right)+p\left(2t_{\mathrm{stup}}+t_{\mathrm{data}}\left[\frac{kn}{p}+\frac{mn}{p}\right]\right)+\frac{2mkn}{p}=O\left(p+\mathrm{log(p)}{N}^2+N^3∕p\right)$

assuming square matrices and a square grid of processes. Similar considerations apply as for mult_replicated above, except that as communication is $O\left(\mathrm{log(p)}{N}^2\right)$ rather than $O\left(\sqrt{p}{N}^2\right)$ we would expect the communication time to increase more slowly with weak scaling.

3. Performance on a single process is less than for DGEMM alone: for square matrices of size N=1000, mult_replicated and mult_replicated2 achieve only 9.40 and 9.98 GFLOP/s respectively, compared to 10.8 GFLOP/s for DGEMM. This is due to the extra time taken to allocate memory and move data between the provided matrices A, B, C and the temporary storage for these matrices.
   Peak FLOPS on XE is 12.0 GFLOP/s - neither algorithm gets as close to this number as does DGEMM. It's very rare for any algorithm to achieve peak FLOPS as this requires an exact match between the instructions issued (load, store, multiply, add) and the capabilities of the architecture, as well the hiding of all memory latencies.

   Weak scaling

   Table 1: weak scaling for matrix multiplication with mult_replicated

   No. Processes	N	Mean Time (s)	GFLOP/s	GFLOP/s per node
   1	1000	0.213	9.40	9.40
   2	1260	0.250	16.0	8.00
   4	1587	0.298	26.8	6.70
   8	2000	0.520	30.8	3.84
   16	2520	0.720	44.5	2.78
   32	3175	1.44	44.3	1.38

   
   

   Table 2: weak scaling for matrix multiplication with mult_replicated2

   No. Processes	N	Mean Time (s)	GFLOP/s	GFLOP/s per node
   1	1000	0.200	9.98	9.98
   2	1260	0.239	16.8	8.40
   4	1587	0.278	28.7	7.18
   8	2000	0.398	40.3	5.03
   16	2520	0.532	60.2	3.76
   32	3175	0.760	84.2	2.63

   
   mult_replicated scales poorly with the number of processes - while total GFLOP/s increases up to 16 processes, GFLOP/s per node drops rapidly. As predicted, weak scaling becomes dominated by communication time, with 32 processes actually achieving lower total GFLOP/s than 16 processes.
   mult_replicated2 scales better, continuing to increase total GFLOP/s to 32 processes. Efficiency still drops rapidly, although less than for mult_replicated as predicted.

4. Strong scaling

   Table 3: strong scaling for matrix multiplication with mult_replicated

   No. Processes	N	Mean Time (s)	GFLOP/s	GFLOP/s per node
   1	4096	12.5	11.0	11.0
   2	4096	6.96	19.7	9.85
   4	4096	4.01	34.2	8.55
   8	4096	2.82	48.6	6.08
   16	4096	2.16	63.6	3.98
   32	4096	2.42	56.8	1.78

   Strong scaling is simpler to compare with the cost model as N is fixed, therefore mult_replicated scales as $O\left(p+1∕p\right)$.
   This suggests a crossover point, a value of p for which increasing communication time outweighs the 1/p reduction in computation. After this point, adding further processes will increase the total time. We see such a crossover point in the measured timings between 16 and 32 processes.

   Table 4: strong scaling for matrix multiplication with mult_replicated2

   No. Processes	N	Mean Time (s)	GFLOP/s	GFLOP/s per node
   1	4096	12.3	11.2	11.2
   2	4096	6.48	21.2	10.6
   4	4096	3.76	36.6	9.14
   8	4096	2.45	56.1	7.01
   16	4096	1.72	80.1	5.01
   32	4096	1.32	104	3.25

   The communication time for mult_replicated2 grows more slowly with the number of processes; it scales as $O\left(p+1∕\mathrm{log(p)}\right)$. This also suggests a crossover point, however it will be at a higher value of p than for mult_replicated. We do not observe any crossover in the measurements: the total time continues to decrease to p=32.

Section 2: SUMMA

1. See mult_summa in attached ass1.c.

2. From the SUMMA paper, the execution times for the algorithm are approximately:

   • Communication (pipelined ring broadcast of A and B): $2\left(N+2\sqrt{p}-3\right)\left(t_{\mathrm{stup}}+\frac{N}{\sqrt{p}}{t}_{\mathrm{data}}\right)=O\left(\sqrt{p}+\frac{N^2}{\sqrt{p}}+N\right)$
   • Computation: $\frac{2mkn}{p}=O\left(N^3∕p\right)$

   In addition, for our code, communications are required to distribute matrices A and B and gather C at the end.

   • $t_{\mathrm{stup}}\left(3p\right)+t_{\mathrm{data}}\left(mk+kn+mn\right)=O\left(p+N^2\right)$

3. Measured times for SUMMA as follows:

Note: times from here on are from the Vayu system

Weak scaling is better than strong scaling - which is surprising as SUMMA is meant to exhibit linear weak scaling! The reason is that the problem size is small, so execution time is dominated by $O\left(p\right)$ scattering and gathering of data not included in the SUMMA algorithm.

Section 3: SUMMA, non-pipelined

1. See mult_summa in attached ass1.c. Note the use of the function do_summa with a function parameter for both pipelined and non-pipelined versions of SUMMA.

2. Using a tree broadcast instead of the pipelined ring broadcast does change the cost model, multiplying SUMMA communications by a factor of log(p).

3. Measured times for non-pipelined SUMMA as follows:

Performance is initially better than mult_summa for both weak and strong scaling tests, probably due to the efficiency of calling MPI_Bcast compared to a series of blocking MPI_Send calls. However the additional factor of log(p) soon takes over and tree broadcast is less efficient for 8 or more processes.

Section 4: (COMP6430 students only) PUMMA

PUMMA uses a block-scattered (also called block-cylic) distribution, in which a matrix is decomposed into blocks of size r × s. SUMMA uses a block distribution, which is a special case of the block-scattered distribution for which $r=⌈m/P⌉$ and $s=⌈n/Q⌉$.
PUMMA uses standard tree broadcast for blocks of the matrix B, while rolling blocks of A between processes in a row, similar to the ring broadcast in the SUMMA algorithm. PUMMA requires more workspace (memory) per node. The block-scattered distribution may be advantageous for applications that use matrix multiplication with other linear algebra operations, as it provides better load balance for some operations, for example LU factorization. While it is possible to redistribute matrices between operations to use the best distribution for each operation, it is preferable to avoid the overhead of redistribution. PUMMA avoids the need for redistribution by using a decomposition that provides a good load balance for many linear algebra operations.