CS267 Assignment 1: Optimize Matrix Multiplication

Due Date: Tuesday February 14, 2012 at 11:59PM

Problem statement

Your task is to optimize matrix multiplication (matmul) code to run fast on a single processor core of NERSC's Franklin cluster.

C := C + A*B

where A, B, and C are n x n matrices. This can be performed using 2n³ floating point operations (n³ adds, n³ multiplies), as in the following pseudocode:

  for i = 1 to n
    for j = 1 to n
      for k = 1 to n
        C(i,j) = C(i,j) + A(i,k) * B(k,j)
      end
    end
  end

Instructions

• You will be paired up in teams of two. Let your instructors know if you still don't have a teammate after Feb 2.
• Your submission should be a gzipped tar archive, formatted (for Team 4) like: team4_hw1.tgz. It should contain:
  • dgemm-blocked.c, a C-language source file containing your implementation of the routine:

      void square_dgemm(int, double*, double*, double*);

    described in pseudocode above. We provide an example dgemm-blocked.c, below.
  • Makefile, only if you modified it. If you modified it, make sure it still correctly builds the provided benchmark.c, which we will use to grade your submission.
  • team4_hw1.pdf, your write-up.
  This link tells you how to use tar to make a .tgz file. Email your .tgz file to the GSIs.
• Your write-up should contain:
  • the names of the people in your group,
  • the optimizations used or attempted,
  • the results of those optimizations,
  • the reason for any odd behavior (e.g., dips) in performance, and
  • how the performance changed when running your optimized code on a different machine.
  For the last requirement, you may run your implementation on another NERSC machine, on your laptop/cellphone, on the cloud, etc.
• Please carefully read the notes for implementation details. Stay tuned to the CS267 Piazza page (sign up first) for updates and clarifications, as well as discussion.
• If you are new to optimizing numerical codes, we recommend reading the papers in the references section.

Notes:

• Your grade will mostly depend on two factors:
  • performance sustained by your codes on Franklin,
  • explanations of the performance features you observed (including what didn't work)
• There are other formulations of matmul (eg, Strassen) that are mathematically equivalent, but perform asymptotically fewer computations - we will not grade submissions that do fewer computations than the 2n³ algorithm. This is actually an optional part of HW1.
• Your code must use double-precision to represent real numbers. A common reference for double-precision matrix multiplication is the dgemm (double-precision general matrix-matrix multiply) routine in the level-3 BLAS. We will compare your implementation with the tuned dgemm implementation in the vendor-provided BLAS library - on Franklin (Cray XT4), we will compare with the Cray LibSci implementation of dgemm. Note that dgemm has a more general interface than square_dgemm - an optional part of HW1 encourages you to explore this richer tuning space.
• You may use any compiler available. We recommend starting with the GNU C compiler (gcc). If you use a compiler other than gcc, you will have to change the provided Makefile, which uses gcc-specific flags. Note that the default compilers, every time you open a new terminal, are PGI - you will have to type "module swap PrgEnv-pgi PrgEnv-gnu" to switch to, eg, GNU compilers. You can type "module list" to see which compiler wrapper you have loaded.
• You may use C99 features when available. The provided benchmark.c uses C99 features, so you must compile with C99 support - for gcc, this means using the flag -std=gnu99 (see the Makefile). Here is the status of C99 functionality in gcc - note that C90 (ANSI C) is fully implemented.
• Besides compiler intrinsic functions and built-ins, your code (dgemm-blocked.c) must only call into the C standard library.
• You should try to use your compiler's automatic vectorizer before manually vectorizing.
  • GNU C provides many extensions, which include intrinsics for vector (SIMD) instructions and data alignment. (Other compilers may have different interfaces.)
  • Ideally your compiler injects the appropriate intrinsics into your code automatically (eg, automatic vectorization and/or automatic data alignment). gcc's auto-vectorizer reports diagnostics that may help you identify if manual vectorization is required.
  • To manually vectorize, you must add compiler intrinsics to your code.
  • Consult your compiler's documentation.
• You may assume that A and B do not alias C; however, A and B may alias each other. It is semantically correct to qualify C (the last argument to square_dgemm) with the C99 restrict keyword. There is a lot online about restrict and pointer-aliasing - this is a good article to start with.
• The matrices are all stored in column-major order, i.e. C_i,j == C(i,j) == C[(i-1)+(j-1)*n], for i=1:n, where n is the number of rows in C. Note that we use 1-based indexing when using mathematical symbols and MATLAB index notation (parentheses), and 0-based indexing when using C index notation (square brackets).
• We will check correctness by the following componentwise error bound:

   |square_dgemm(n,A,B,0) - A*B| < eps*n*|A|*|B|. 

  where eps := 2^–52 = 2.2 * 10^–16 is the machine epsilon.
• The target processor is a quad-core AMD Opteron 1356 (see AMD documentation). Its cores have separate 128-bit wide pipelines for addition and multiplication clocked at 2.3GHz. This yields 2 double-precision (ie, 64-bit) flops per pipeline * 2 pipelines * 2.3 GHz = 9.2 Gflops/s. The following graph shows the performance of the supplied BLAS, blocked, and naive implementations (compiled using gcc with optimizations: -O3 -ffast-math -funroll-loops -march=amdfam10)
  
  
  
  Part of your grade depends on how fast your code runs, relative to the supplied codes. In previous years, some students actually beat the vendor-optimized BLAS.

Optional:

These parts are not graded. You should be satisfied with your square_dgemm results and write-up before beginning an optional part.

• Implement Strassen matmul. Consider switching over to the three-nested-loops algorithm when the recursive subproblems are small enough.
• Support the dgemm interface (ie, rectangular matrices, transposing, scalar multiples).
• Try float (single-precision). This means you can use 4-way SIMD parallelism on Franklin.
• Try complex numbers (single- and double-precision) - note that complex numbers are part of C99 and supported in gcc. This forum thread gives advice on vectorizing complex multiplication with the conventional approach - but note that there are other algorithms for this operation.
• Optimize your matmul for the case when the inputs are symmetric. Consider conventional and packed symmetric storage.

Source files

We provide two simple implementations for you to start with: a naive three-loop implementation similar to the pseudocode above, and a more cache-efficient blocked implementation.

The necessary files are in cs267_hw1.tgz. Included are the following:

dgemm-naive.c A naive implementation of matrix multiply using three nested loops, dgemm-blocked.c A simple blocked implementation of matrix multiply, dgemm-blas.c A wrapper for the vendor's optimized BLAS implementation of matrix multiply (default: Cray LibSci), benchmark.c The driver program that measures the runtime and verifies the correctness by comparing with the vendor's implementation, Makefile A simple makefile to build the executables, job-blas, job-blocked, job-naive Scripts to run the executables on Franklin compute nodes. For example, type "qsub job-blas" to benchmark the BLAS version.

The documentation for Franklin's programming environment can be found below.

References

• Goto, K., and van de Geijn, R. A. 2008. Anatomy of High-Performance Matrix Multiplication, ACM Transactions on Mathematical Software 34, 3, Article 12.
  (Note: explains the design decisions for the GotoBLAS dgemm implementation, which also apply to your code.)
• Chellappa, S., Franchetti, F., and Püschel, M. 2008. How To Write Fast Numerical Code: A Small Introduction, Lecture Notes in Computer Science 5235, 196–259.
  (Note: how to write C code for modern compilers and memory hierarchies, so that it runs fast. Recommended reading, especially for newcomers to code optimization.)
• Bilmes, et al. The PHiPAC (Portable High Performance ANSI C) Page for BLAS3 Compatible Fast Matrix Matrix Multiply.
  (Note: PHiPAC is a code-generating autotuner for matmul that started as a submission for this HW in a previous semester of CS267. Also see ATLAS; both are good examples if you are considering code generation strategies.)
• Lam, M. S., Rothberg, E. E, and Wolf, M. E. 1991. The Cache Performance and Optimization of Blocked Algorithms, ASPLOS'91, 63–74.
  (Note: clearly explains cache blocking, supported by with performance models.)
• Notes on vectorizing with SSE intrinsics, from lecture 2/9/12, here

Documentation:

• Franklin's programming environment documentation
• AMD architecture documentation.
• GCC documentation - Franklin's default version currently is GCC 4.5.3, but offers older versions too.
• Intel Intrinsic Guide (scroll down). This Java program has a complete reference of all SIMD intrinsics on Intel architectures. Note that Franklin supports SSE, SSE2, SSE3, and SSE4a (not SSE4) instruction sets.

Results

Below are results recorded on Franklin using the provided benchmark. Performance was reproducible to within 5%, so if you feel your performance is misrepresented, please re-run your submitted code to make sure, and then contact the GSIs (cs267.sp12@gmail.com) with this data.

Note that

• Team -1 = naive blocked code (b=41) compiled with GNU, "-O1" optimization
• Team 0 = Cray LibSci DGEMM
• Team 17 = naive blocked code (b=41) compiled with GNU, "-O3 -ffast-math -funroll-loops -march=amdfam10" optimization