You will optimize a routine to multiply square matrices. Matrix multiplication is a basic building block in many scientific computations; and since it is an O(n3) algorithm, these codes often spend a lot of their time in matrix multiplication. The most naive code to multiply matrices is short, sweet, simple, and very slow:
for i = 1 to n
for j = 1 to m
for k = 1 to m
C(i,j) = C(i,j) + A(i,k) * B(k,j)
end
end
end
We're providing implementations of the trivial unoptimized matrix multiply code in C and Fortran, and a very simple blocked implementation in C. We're also providing a wrapper for the BLAS (Basic Linear Algebra Subroutines) library. The BLAS is a standard interface for building blocks like matrix multiplication, and many vendors provide optimized versions of the BLAS for their platforms. You may want to compare your routines to codes that others have optimized; if you do, think about the relative difficulty of using someone else's library to writing and optimizing your own routine. Knowing when and how to make use of standard libraries is an important skill in building fast programs!
The unblocked routines we provide do as poorly as you'd expect, and the blocked routine provided doesn't do much better.
The performance of the two routines is shown above. The naive unblocked routine is the "three loops" curve, and is just the standard three nested loops. The blocked code (block size of 56, selected as the "best" in an exhaustive search of square tile sizes up to 64) is shown in green.
Characterizing the peak performance is slightly tricky. The processors here are running at 1.3GHz and have two floating-point units, so the maximum rate of executing instructions is 2.6GHz or 2.6 GFlop/s. Each of those instructions can be a fused multiply-add, effectively running two operations at once. A matrix multiplication can take advantage of that, raising the peak performance to roughly 5.2 GFlop/s. These routines are achieving around 30% of peak.
However, the unblocked routine's performance completely dies on matrices larger than 511 x 511, while the blocked routine gives more consistent performance for all the tested sizes. Note the performance dips at powers of 2 (cache conflicts). The machine used for this diagram has a 3 MiB external cache, and 3 matrices * 512x512 entries * 8 bytes/entry is 6 MiB, leading to the first mystery.
You need to write a dgemm.c that contains a function with the following C signature:
void
square_dgemm (const unsigned M,
const double *A, const double *B, double *C)
If you would prefer, you can also write a Fortran function in a file fdgemm.f:
subroutine sdgemm(M, A, B, C)
c
c .. Parameters ..
integer M
double precision A(M,M)
double precision B(M,M)
double precision C(M,M)
Note that the matrices are stored in column-major order. Also, your program will actually be doing a multiply and add operation:
C := A*B + C
Look at the code in basic_dgemm.c if you find this confusing.
The necessary files are in tuning-matmul.tar. Included are the following:
We will be testing on the 1.3GHz Itanium2 nodes. The nodes listed in nodes-fast are currently 1.3GHz, dual-processor nodes. To restrict what servers are used by gexec, you can set the GEXEC_SVRS environment variable. For instance, to use c42 and c19, you set GEXEC_SVRS to "c42 c19":
bash% GEXEC_SVRS="c42 c19" gexec uname -n 0 c42 1 c19
You can use GEXEC_SVRS="`cat nodes-fast`" to pick from the faster nodes.
The y-axis in the gnuplot script plots is labeled MFlop/s, but really it should be "MFlop/s assuming you were using the standard algorithm." If you use something like Strassen's, we will still measure time/(2n3). However, the numerical error checking may haunt you.
You will get a lot of improvement from tiling, but you may want to play with other optimizations, too. Strassen-like algorithms, copy optimizations, recursive data layout, ... you can get some pretty elaborate optimization strategies. We'll cover some of the possibilities in discussion, but you might want to do some independent reading, too. The off-limit "optimizations" are those which weaken the floating point system: No flush-to-zero mode.
Your group needs to submit your group's dgemm.c (or fdgemm.c), Makefile (for compiler options) and a write-up. Your write-up should contain
To show the results of your optimizations, include a graph comparing your dgemm.c with the included basic_dgemm.c. For the last requirement, try your tuned dgemm.c on another hardware platform (like the x86 Millennium nodes, a Sun box, etc.) and explain why it performs poorly. Your explanations should rely heavily on knowledge of the memory hierarchy. (Benchmark graphs help.)
Please tar up your group's dgemm.c, write-up, and associated files and mail the TA either an encoded tar file (uuencode or Base-64) or a URL from which we can retrieve the tar file.
The race will be held during a discussion section, exact date to be announced.
Main CS267 page, and the TA's CS267 page
E. Jason Riedy