Differences between revisions 46 and 47

Deletions are marked like this. Additions are marked like this.
Line 424: Line 424:
Copy the contents of file `MMult13.c` into a file named `MMult14.c` and change the contents: Copy the contents of file `MMult14.c` into a file named `MMult14.c` and change the contents:

Contents

  1. Set Up
    1. Files
    2. Which Compiler?
    3. Initialization
    4. Optimization 1
  2. Optimization 2-9: unrolling i
    1. Optimization 2
    2. Optimization 3
    3. Optimization 4
    4. Optimization 5
    5. Optimization 6
    6. Optimization 7
    7. Optimization 8
    8. Optimization 9
  3. Optimization 10-16: unrolling j
    1. Optimization 10
    2. Optimization 11
    3. Optimization 12
    4. Optimization 13
    5. Optimization 14
    6. Optimization 15
    7. Optimization 16
Optimizing Matrix-Matrix Multiplication

This page leads one, step by step, through the optimization of matrix-matrix multiplication. For now, it is assumed that the reader has a Linux account on a Pentium4 based computer. We will use the gcc compiler as part of the exercise.

Set Up

  • Download OptimizeGemm.tar.gz

  • Uncompress by executing "gunzip OptimizeGemm.tar.gz"

  • Expand the tar file by executing "tar OptimizeGemm.tar"

  • Change into the directory that is created by executing "cd OptimizeGemm"

Files

In the directory OptimizeGemm you will find the following files that you will use to systematically optimize the matrix-matrix multiplication operation:

  • makefile

    makefile

    The makefile that describes how to compile, link, and execute the driver/implementations.

    Test driver

    parameters.h

    File that holds parameters that control what data is collected

    test_MMult.c

    driver routine that tests and times the different implementations. This routine executes a reference implementation and the current optimization to be timed. Parameters for this routine are initialized in parameters.h. In particular, in that file it is indicated how many times to repeat each experiment (problem size) and how each of the three dimensions ($ m $, $ n $, and $ k $ are tied to the problem size being timed.

    Matrix multiplication implementations

    REF_MMult.c

    Reference implementation used to check correctness

    MMult0.c

    Version 0: simplest implementation

    MMult1.c

    Optimization 1

    Utility routines

    compare_matrices.c

    Compares the contents of two matrices and returns the maximum absolute difference

    copy_matrix.c

    Copies one matrix to another

    dclock.c

    Returns elapsed time in seconds

    random_matrix.c

    Generates a random matrix

    transpose_matrix.c

    Transposes the contents of one matrix into another matrix

    print_matrix.c

    Prints the contents of a matrix

    Plotting the results

    PlotAll.m

    Plots graphs corresponding to the data in files output_old.m and output_new.m

    proc_parameters.m

    File in which parameters about the architecture are given

These last routines allow one to use Matlab to plot the performance of two different implementations.

Which Compiler?

In this exercise, we use the gcc (Gnu C) compiler with optimization level -O2. (See the makefile.) This is neither the best compiler nor the best optimization level. The reason is that with this compiler and level of optimization, we have a certain level of control:

  • Had we used the intel compiler, the "simple loops" in MMult0.c would probably have yielded the best performance: This compiler is very good at optimizing a triple loop.
  • Had we used -O3 (optimization level 3), the gnu compiler would have more aggressively optimized, making the step-by-step optimizations we demonstrate much less predictable.

Initialization

Make sure that the makefile starts with the following lines: 

  • OLD  := 0
NEW  := 0

This indicates that the performance of the version of matrix-matrix multiplication in MMult0.c is measured (by virtue of the statement OLD  :=0).

Next, to make sure that when plotting the graphs are properly scaled, set certain parameters in the file proc_parameters.m. See the comments in that file. (Setting these parameters will ensure that when plotting the y-axis ranges from 0 to the peak performance of the architecture.)

Execute

  • make run
    This will compile, link, and execute the test driver, linking to the implementation in MMult0.c. The performance data is saved in file output0.m.

  • more output0.m
    This will display the contents of the output file output0.m. It should look something like 

    • MY_MMult = [
20 4.953912e-01 0.000000e+00 
40 6.276695e-01 0.000000e+00 
60 6.324479e-01 0.000000e+00 
80 6.295030e-01 0.000000e+00 
  [ lines deleted ]
360 3.791921e-01 0.000000e+00 
380 3.228467e-01 0.000000e+00 
400 2.622753e-01 0.000000e+00 
];

    The first column equals the problem size. The second column the elapsed time (in seconds) when a matrix-matrix multiply with the indicated problem size ($m=n=k$ is executed. The last column reports the maximum absolute difference encountered between the implementation in REF_MMult.c and MMult0.c. It should be close to 0.00000e+00 although as different optimizations are added the difference may not be perfectly zero.

  • make plot
    This will execute the commands in PlotAll using Matlab, creating a graph in the file graph_0_vs_0.eps.

  • gv graph_0_vs_0.eps
    This displays the graph.

The performance graph will look something like

Notice that the two curves are right on top of each other because data for the same im plementation are being compared. From the fact that the top of the graph represents peak performance, it is obvious that this simple implementation achieves only a fraction of the ideal performance.

More performance graphs

Optimization 1

Change the first lines in the makefile to 

  • OLD  := 0
NEW  := 1

  • make run

  • make plot

  • gv graph_0_vs_1.eps

This time the performance graph will look something like

Compare the two implementations in MMult0.c and MMult1.c : 

  • /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( i=0; i<m; i++ ){
    for ( j=0; j<n; j++ ){
      for ( p=0; p<k; p++ ){
        C( i,j ) += A( i,p ) * B( p,j );
      }
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    for ( i=0; i<m; i++ ){
      register double c0=0.0;

      for ( p=0; p<k; p++ ){
        c0 += A( i,p ) * B( p,j );
      }
      
      C( i,j ) += c0;
    }
  }
}

The difference in performance can be attributed to the fact that the update to C( i,j ) is accumulated in a register, and the order of the loops indexed by i and j have been exchanged. (You may want to see what happens if you only make one of these changes, rather than both.)

It should be obvious now how setting OLD and NEW in the makefile allows one to compare different implementations.

Optimization 2-9: unrolling i

Optimization 2

Copy the contents of file MMult1.c into a file named MMult2.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    for ( i=0; i<m; i++ ){
      register double c0=0.0;

      for ( p=0; p<k; p++ ){
        c0 += A( i,p ) * B( p,j );
      }
      
      C( i,j ) += c0;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;

      for ( p=0; p<k; p++ ){
        c0 += A( i,p ) * B( p,j );
        c1 += A( i+1,p ) * B( p,j );
        c2 += A( i+2,p ) * B( p,j );
        c3 += A( i+3,p ) * B( p,j );
      }
      
      C( i,j )   += c0;
      C( i+1,j ) += c1;
      C( i+2,j ) += c2;
      C( i+3,j ) += c3;
    }
  }
}

Change the first lines in the makefile to 

  • OLD  := 1
NEW  := 2

  • make run

  • make plot

  • gv graph_1_vs_2.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that the loop has been unrolled do that four elements of C, all in the same column, are computed. This decreases the overhead for computing the loop index p.

Optimization 3

Copy the contents of file MMult2.c into a file named MMult3.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;

      for ( p=0; p<k; p++ ){
        c0 += A( i,p ) * B( p,j );
        c1 += A( i+1,p ) * B( p,j );
        c2 += A( i+2,p ) * B( p,j );
        c3 += A( i+3,p ) * B( p,j );
      }
      
      C( i,j )   += c0;
      C( i+1,j ) += c1;
      C( i+2,j ) += c2;
      C( i+3,j ) += c3;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        c0 = c0 + ap[0] * *bp;
        c1 = c1 + ap[1] * *bp;
        c2 = c2 + ap[2] * *bp;
        c3 = c3 + ap[3] * *bp;

        bp++;
        ap += lda;
      }
      
      C( i,j )   += c0;
      C( i+1,j ) += c1;
      C( i+2,j ) += c2;
      C( i+3,j ) += c3;
    }
  }
}

Change the first lines in the makefile to 

  • OLD  := 2
NEW  := 3

  • make run

  • make plot

  • gv graph_2_vs_3.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimization 4

Copy the contents of file MMult3.c into a file named MMult4.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        c0 = c0 + ap[0] * *bp;
        c1 = c1 + ap[1] * *bp;
        c2 = c2 + ap[2] * *bp;
        c3 = c3 + ap[3] * *bp;

        bp++;
        ap += lda;
      }
      
      C( i,j )   += c0;
      C( i+1,j ) += c1;
      C( i+2,j ) += c2;
      C( i+3,j ) += c3;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for cogmputing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        register double bpj = *bp;

        c0 = c0 + ap[0] * bpj;
        c1 = c1 + ap[1] * bpj;
        c2 = c2 + ap[2] * bpj;
        c3 = c3 + ap[3] * bpj;

        bp++;
        ap+=lda;
      }
      
      C( i,j )   += c0;
      C( i+1,j ) += c1;
      C( i+2,j ) += c2;
      C( i+3,j ) += c3;
    }
  }
}

Change the first lines in the makefile to 

  • OLD  := 3
NEW  := 4

  • make run

  • make plot

  • gv graph_3_vs_4.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimization 5

Copy the contents of file MMult4.c into a file named MMult5.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for cogmputing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        register double bpj = *bp;

        c0 = c0 + ap[0] * bpj;
        c1 = c1 + ap[1] * bpj;
        c2 = c2 + ap[2] * bpj;
        c3 = c3 + ap[3] * bpj;

        bp++;
        ap+=lda;
      }
      
      C( i,j )   += c0;
      C( i+1,j ) += c1;
      C( i+2,j ) += c2;
      C( i+3,j ) += c3;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    double *cp = &C( 0, j );

    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        register double bpj = *bp;

        c0 = c0 +  ap[0] * bpj;
        c1 = c1 +  ap[1] * bpj;
        c2 = c2 +  ap[2] * bpj;
        c3 = c3 +  ap[3] * bpj;

        bp++;
        ap+=lda;
      }
      
      cp[ 0 ] += c0;
      cp[ 1 ] += c1;
      cp[ 2 ] += c2;
      cp[ 3 ] += c3;
      
      cp += 4;
    }
  }
}

Change the first lines in the makefile to 

  • OLD  := 4
NEW  := 5

  • make run

  • make plot

  • gv graph_4_vs_5.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimization 6

Copy the contents of file MMult5.c into a file named MMult6.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    double *cp = &C( 0, j );

    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        register double bpj = *bp;

        c0 = c0 +  ap[0] * bpj;
        c1 = c1 +  ap[1] * bpj;
        c2 = c2 +  ap[2] * bpj;
        c3 = c3 +  ap[3] * bpj;

        bp++;
        ap+=lda;
      }
      
      cp[ 0 ] += c0;
      cp[ 1 ] += c1;
      cp[ 2 ] += c2;
      cp[ 3 ] += c3;
      
      cp += 4;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    double *cp = &C( 0, j );

    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        register double bpj = *bp;
        register double r0, r1, r2;

        r0 = ap[ 0 ];
        r1 = ap[ 1 ];
        r2 = ap[ 2 ];
        c0 = c0 +  r0 * bpj;

        r0 = ap[ 3 ];
        c1 = c1 +  r1 * bpj;
        c2 = c2 +  r2 * bpj;
        c3 = c3 +  r0 * bpj;

        bp++;
        ap+=lda;
      }
      
      cp[ 0 ] += c0;
      cp[ 1 ] += c1;
      cp[ 2 ] += c2;
      cp[ 3 ] += c3;
      
      cp += 4;
    }
  }
}

Change the first lines in the makefile to 

  • OLD  := 5
NEW  := 6

  • make run

  • make plot

  • gv graph_5_vs_6.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimization 7

Copy the contents of file MMult6.c into a file named MMult7.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    double *cp = &C( 0, j );

    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        register double bpj = *bp;
        register double r0, r1, r2;

        r0 = ap[ 0 ];
        r1 = ap[ 1 ];
        r2 = ap[ 2 ];
        c0 = c0 +  r0 * bpj;

        r0 = ap[ 3 ];
        c1 = c1 +  r1 * bpj;
        c2 = c2 +  r2 * bpj;
        c3 = c3 +  r0 * bpj;

        bp++;
        ap+=lda;
      }
      
      cp[ 0 ] += c0;
      cp[ 1 ] += c1;
      cp[ 2 ] += c2;
      cp[ 3 ] += c3;
      
      cp += 4;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#include "parameters.h"

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

#define min( x, y ) ( x < y ? x : y )

#define MB 48
#define KB 48

/* MY_MMult_inner prototype */
void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc );
                                                                        

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, ib, pb;

  for ( p=0; p<k; p+=KB ){
    pb = min( k-p, KB );
    for ( i=0; i<m; i+=MB ){
      ib = min( m-i, MB );

      MY_MMult_inner( ib, n, pb, 
                      &A( i, p ), lda, 
                      &B( p, 0 ), ldb, 
                      &C( i, 0 ), ldc );
    }
  }
}


  


/* Inner kernel for computing C = A * B + C */

void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    double *cp = &C( 0, j );

    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        register double bpj = *bp;

        c0 = c0 + ap[0] * bpj;
        c1 = c1 + ap[1] * bpj;
        c2 = c2 + ap[2] * bpj;
        c3 = c3 + ap[3] * bpj;

        bp++;
        ap+=lda;
      }
      
      cp[ 0 ] += c0;
      cp[ 1 ] += c1;
      cp[ 2 ] += c2;
      cp[ 3 ] += c3;
      
      cp += 4;
    }
  }
}

Change the first lines in the makefile to 

  • OLD  := 6
NEW  := 7

  • make run

  • make plot

  • gv graph_6_vs_7.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimization 8

Copy the contents of file MMult7.c into a file named MMult8.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#include "parameters.h"

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

#define min( x, y ) ( x < y ? x : y )

#define MB 48
#define KB 48

/* MY_MMult_inner prototype */
void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc );
                                                                        

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, ib, pb;

  for ( p=0; p<k; p+=KB ){
    pb = min( k-p, KB );
    for ( i=0; i<m; i+=MB ){
      ib = min( m-i, MB );

      MY_MMult_inner( ib, n, pb, 
                      &A( i, p ), lda, 
                      &B( p, 0 ), ldb, 
                      &C( i, 0 ), ldc );
    }
  }
}


  


/* Inner kernel for computing C = A * B + C */

void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    double *cp = &C( 0, j );

    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        register double bpj = *bp;

        c0 = c0 + ap[0] * bpj;
        c1 = c1 + ap[1] * bpj;
        c2 = c2 + ap[2] * bpj;
        c3 = c3 + ap[3] * bpj;

        bp++;
        ap+=lda;
      }
      
      cp[ 0 ] += c0;
      cp[ 1 ] += c1;
      cp[ 2 ] += c2;
      cp[ 3 ] += c3;
      
      cp += 4;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#include "parameters.h"

#define A(i,j) a[ (j)*lda + (i) ]
#define AT(i,j) at[ (j)*(k) + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

#define min( x, y ) ( x < y ? x : y )

#define MB 48
#define KB 48

/* MY_MMult_inner prototype */
void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc );
                                                                        

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, ib, pb;

  for ( p=0; p<k; p+=KB ){
    pb = min( k-p, KB );
    for ( i=0; i<m; i+=MB ){
      ib = min( m-i, MB );

      MY_MMult_inner( ib, n, pb, 
                      &A( i, p ), lda, 
                      &B( p, 0 ), ldb, 
                      &C( i, 0 ), ldc );
    }
  }
}


  


/* Inner kernel for computing C = A * B + C */

void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;
  double at[ MB * KB ];

  transpose_matrix( m, k, a, lda, at, k );

  for ( j=0; j<n; j++ ){
    double *cp = &C( 0, j );

    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *a0p = &AT(0,i  ), *a1p = &AT(0,i+1), 
             *a2p = &AT(0,i+2), *a3p = &AT(0,i+3), 
             *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        register double bpj = *bp;

        c0 = c0 +  *a0p++ * bpj;
        c1 = c1 +  *a1p++ * bpj;
        c2 = c2 +  *a2p++ * bpj;
        c3 = c3 +  *a3p++ * bpj;

        bp++;
      }
      
      cp[ 0 ] += c0;
      cp[ 1 ] += c1;
      cp[ 2 ] += c2;
      cp[ 3 ] += c3;
      
      cp += 4;
    }
  }
}

Change the first lines in the makefile to 

  • OLD  := 7
NEW  := 8

  • make run

  • make plot

  • gv graph_7_vs_8.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimizing Optimization 8

By transposing the matrix, it turns out that a larger block size can be used. (The reason for this is that the block of A can be increased beyond the size of the L1 cache.) Try changing MB = 144 and KB = 144 in file MMult8.c and reexecute make run, make plot, and gv graph_7_vs_8.eps. On my machine this yields

Optimization 9

Copy the contents of file MMult8.c into a file named MMult9.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#include "parameters.h"

#define A(i,j) a[ (j)*lda + (i) ]
#define AT(i,j) at[ (j)*(k) + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

#define min( x, y ) ( x < y ? x : y )

#define MB 48
#define KB 48

/* MY_MMult_inner prototype */
void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc );
                                                                        

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, ib, pb;

  for ( p=0; p<k; p+=KB ){
    pb = min( k-p, KB );
    for ( i=0; i<m; i+=MB ){
      ib = min( m-i, MB );

      MY_MMult_inner( ib, n, pb, 
                      &A( i, p ), lda, 
                      &B( p, 0 ), ldb, 
                      &C( i, 0 ), ldc );
    }
  }
}


  


/* Inner kernel for computing C = A * B + C */

void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;
  double at[ MB * KB ];

  transpose_matrix( m, k, a, lda, at, k );

  for ( j=0; j<n; j++ ){
    double *cp = &C( 0, j );

    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *a0p = &AT(0,i  ), *a1p = &AT(0,i+1), 
             *a2p = &AT(0,i+2), *a3p = &AT(0,i+3), 
             *bp=&B(0,j);

      for ( p=0; p<k; p++ ){
        register double bpj = *bp;

        c0 = c0 +  *a0p++ * bpj;
        c1 = c1 +  *a1p++ * bpj;
        c2 = c2 +  *a2p++ * bpj;
        c3 = c3 +  *a3p++ * bpj;

        bp++;
      }
      
      cp[ 0 ] += c0;
      cp[ 1 ] += c1;
      cp[ 2 ] += c2;
      cp[ 3 ] += c3;
      
      cp += 4;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#include "parameters.h"

#define A(i,j) a[ (j)*lda + (i) ]
#define AT(i,j) at[ (j)*(k) + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

#define min( x, y ) ( x < y ? x : y )

#define MB 144
#define KB 144

/* MY_MMult_inner prototype */
void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc );
                                                                        

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, ib, pb;

  for ( p=0; p<k; p+=KB ){
    pb = min( k-p, KB );
    for ( i=0; i<m; i+=MB ){
      ib = min( m-i, MB );

      MY_MMult_inner( ib, n, pb, 
                      &A( i, p ), lda, 
                      &B( p, 0 ), ldb, 
                      &C( i, 0 ), ldc );
    }
  }
}


  


/* Inner kernel for computing C = A * B + C */

void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;
  double at[ MB * KB ];

  transpose_matrix( m, k, a, lda, at, k );

  for ( j=0; j<n; j++ ){
    double *cp = &C( 0, j );

    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *a0p = &AT(0,i), *a1p = &AT(0,i+1), 
             *a2p = &AT(0,i+2), *a3p = &AT(0,i+3), 
             *bp=&B(0,j);

      for ( p=0; p<k; p+=4 ){
        register double bpj = bp[0];

        c0 = c0 +  a0p[0] * bpj;
        c1 = c1 +  a1p[0] * bpj;
        c2 = c2 +  a2p[0] * bpj;
        c3 = c3 +  a3p[0] * bpj;

        bpj = bp[1];

        c0 = c0 +  a0p[1] * bpj;
        c1 = c1 +  a1p[1] * bpj;
        c2 = c2 +  a2p[1] * bpj;
        c3 = c3 +  a3p[1] * bpj;

        bpj = bp[2];

        c0 = c0 +  a0p[2] * bpj;
        c1 = c1 +  a1p[2] * bpj;
        c2 = c2 +  a2p[2] * bpj;
        c3 = c3 +  a3p[2] * bpj;

        bpj = bp[3];

        c0 = c0 +  a0p[3] * bpj;
        c1 = c1 +  a1p[3] * bpj;
        c2 = c2 +  a2p[3] * bpj;
        c3 = c3 +  a3p[3] * bpj;

        bp+=4;

        a0p+=4;
        a1p+=4;
        a2p+=4;
        a3p+=4;
      }
      
      cp[ 0 ] += c0;
      cp[ 1 ] += c1;
      cp[ 2 ] += c2;
      cp[ 3 ] += c3;
      
      cp += 4;
    }
  }
}

Change the first lines in the makefile to 

  • OLD  := 8
NEW  := 9

  • make run

  • make plot

  • gv graph_8_vs_9.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimization 10-16: unrolling j

Optimization 10

Copy the contents of file MMult2.c into a file named MMult10.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j++ ){
    for ( i=0; i<m; i+=4 ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;

      for ( p=0; p<k; p++ ){
        c0 += A( i,p ) * B( p,j );
        c1 += A( i+1,p ) * B( p,j );
        c2 += A( i+2,p ) * B( p,j );
        c3 += A( i+3,p ) * B( p,j );
      }
      
      C( i,j )   += c0;
      C( i+1,j ) += c1;
      C( i+2,j ) += c2;
      C( i+3,j ) += c3;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j+=4 ){
    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;

      for ( p=0; p<k; p++ ){
        c0 += A( i,p ) * B( p,j );
        c1 += A( i,p ) * B( p,j+1 );
        c2 += A( i,p ) * B( p,j+2 );
        c3 += A( i,p ) * B( p,j+3 );
      }
      
      C( i,j )   += c0;
      C( i,j+1 ) += c1;
      C( i,j+2 ) += c2;
      C( i,j+3 ) += c3;
    }
  }
}

Change the first lines in the makefile to 

  • OLD  := 2
NEW  := 10

  • make run

  • make plot

  • gv graph_2_vs_10.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimization 11

Copy the contents of file MMult10.c into a file named MMult11.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j+=4 ){
    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;

      for ( p=0; p<k; p++ ){
        c0 += A( i,p ) * B( p,j );
        c1 += A( i,p ) * B( p,j+1 );
        c2 += A( i,p ) * B( p,j+2 );
        c3 += A( i,p ) * B( p,j+3 );
      }
      
      C( i,j )   += c0;
      C( i,j+1 ) += c1;
      C( i,j+2 ) += c2;
      C( i,j+3 ) += c3;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j+=4 ){
    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p++ ){
        c0 = c0 + *ap * *b0p++;
        c1 = c1 + *ap * *b1p++;
        c2 = c2 + *ap * *b2p++;
        c3 = c3 + *ap * *b3p++;

        ap += lda;
      }
      
      C( i,j )   += c0;
      C( i,j+1 ) += c1;
      C( i,j+2 ) += c2;
      C( i,j+3 ) += c3;
    }
  }
}

Compare also Optimization 3 (above) to this Optimization 11. (They are similar, except that the same optimizations are applied to Optization 2 and 10, respectively.)

Change the first lines in the makefile to 

  • OLD  := 3
NEW  := 11

  • make run

  • make plot

  • gv graph_3_vs_11.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimization 12

Copy the contents of file MMult11.c into a file named MMult12.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j+=4 ){
    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p++ ){
        c0 = c0 + *ap * *b0p++;
        c1 = c1 + *ap * *b1p++;
        c2 = c2 + *ap * *b2p++;
        c3 = c3 + *ap * *b3p++;

        ap += lda;
      }
      
      C( i,j )   += c0;
      C( i,j+1 ) += c1;
      C( i,j+2 ) += c2;
      C( i,j+3 ) += c3;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j+=4 ){
    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p++ ){
        register double aip = *ap;

        c0 = c0 + aip * *b0p++;
        c1 = c1 + aip * *b1p++;
        c2 = c2 + aip * *b2p++;
        c3 = c3 + aip * *b3p++;

        ap += lda;
      }
      
      C( i,j )   += c0;
      C( i,j+1 ) += c1;
      C( i,j+2 ) += c2;
      C( i,j+3 ) += c3;
    }
  }
}

Compare also Optimization 4 (above) to this Optimization 12. (They are similar, except that the same optimizations are applied to Optization 3 and 11, respectively.)

Change the first lines in the makefile to 

  • OLD  := 4
NEW  := 12

  • make run

  • make plot

  • gv graph_4_vs_12.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimization 13

Copy the contents of file MMult12.c into a file named MMult13.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j+=4 ){
    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p++ ){
        register double aip = *ap;

        c0 = c0 + aip * *b0p++;
        c1 = c1 + aip * *b1p++;
        c2 = c2 + aip * *b2p++;
        c3 = c3 + aip * *b3p++;

        ap += lda;
      }
      
      C( i,j )   += c0;
      C( i,j+1 ) += c1;
      C( i,j+2 ) += c2;
      C( i,j+3 ) += c3;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j+=4 ){
    double *c0p = &C(0,j), *c1p = &C(0,j+1), 
           *c2p = &C(0,j+2), *c3p = &C(0,j+3);

    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p++ ){
        register double aip = *ap;

        c0 = c0 + aip * *b0p++;
        c1 = c1 + aip * *b1p++;
        c2 = c2 + aip * *b2p++;
        c3 = c3 + aip * *b3p++;

        ap += lda;
      }
      
      *c0p++ += c0;
      *c1p++ += c1;
      *c2p++ += c2;
      *c3p++ += c3;
    }
  }
}

Compare also Optimization 5 (above) to this Optimization 13. (They are similar, except that the same optimizations are applied to Optization 4 and 12, respectively.)

Change the first lines in the makefile to 

  • OLD  := 5
NEW  := 13

  • make run

  • make plot

  • gv graph_5_vs_13.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimization 14

Copy the contents of file MMult13.c into a file named MMult14.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j+=4 ){
    double *c0p = &C(0,j), *c1p = &C(0,j+1), 
           *c2p = &C(0,j+2), *c3p = &C(0,j+3);

    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p++ ){
        register double aip = *ap;

        c0 = c0 + aip * *b0p++;
        c1 = c1 + aip * *b1p++;
        c2 = c2 + aip * *b2p++;
        c3 = c3 + aip * *b3p++;

        ap += lda;
      }
      
      *c0p++ += c0;
      *c1p++ += c1;
      *c2p++ += c2;
      *c3p++ += c3;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

#define min( x, y ) ( x < y ? x : y )

#define MB 48
#define KB 48

/* MY_MMult_inner prototype */
void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc );
                                                                        

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, ib, pb;

  for ( p=0; p<k; p+=KB ){
    pb = min( k-p, KB );
    for ( i=0; i<m; i+=MB ){
      ib = min( m-i, MB );

      MY_MMult_inner( ib, n, pb, 
                      &A( i, p ), lda, 
                      &B( p, 0 ), ldb, 
                      &C( i, 0 ), ldc );
    }
  }
}


  


/* Inner kernel for computing C = A * B + C */

void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                          double *b, int ldb,
                                          double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j+=4 ){
    double *c0p = &C(0,j), *c1p = &C(0,j+1), *c2p = &C(0,j+2), *c3p = &C(0,j+3);

    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p++ ){
        register double aip = *ap;

        c0 = c0 + aip * *b0p++;
        c1 = c1 + aip * *b1p++;
        c2 = c2 + aip * *b2p++;
        c3 = c3 + aip * *b3p++;

        ap += lda;
      }
      
      *c0p++ += c0;
      *c1p++ += c1;
      *c2p++ += c2;
      *c3p++ += c3;
    }
  }
}

Compare also Optimization 7 (above) to this Optimization 14. (They are similar, except that the same optimizations are applied to Optization 4 and 13, respectively.)

Change the first lines in the makefile to 

  • OLD  := 7
NEW  := 14

  • make run

  • make plot

  • gv graph_7_vs_14.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimization 15

Copy the contents of file MMult14.c into a file named MMult14.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

#define min( x, y ) ( x < y ? x : y )

#define MB 48
#define KB 48

/* MY_MMult_inner prototype */
void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc );
                                                                        

/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, ib, pb;

  for ( p=0; p<k; p+=KB ){
    pb = min( k-p, KB );
    for ( i=0; i<m; i+=MB ){
      ib = min( m-i, MB );

      MY_MMult_inner( ib, n, pb, 
                      &A( i, p ), lda, 
                      &B( p, 0 ), ldb, 
                      &C( i, 0 ), ldc );
    }
  }
}


  


/* Inner kernel for computing C = A * B + C */

void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                          double *b, int ldb,
                                          double *c, int ldc )
{
  int i, j, p;

  for ( j=0; j<n; j+=4 ){
    double *c0p = &C(0,j), *c1p = &C(0,j+1), *c2p = &C(0,j+2), *c3p = &C(0,j+3);

    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &A(i,0), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p++ ){
        register double aip = *ap;

        c0 = c0 + aip * *b0p++;
        c1 = c1 + aip * *b1p++;
        c2 = c2 + aip * *b2p++;
        c3 = c3 + aip * *b3p++;

        ap += lda;
      }
      
      *c0p++ += c0;
      *c1p++ += c1;
      *c2p++ += c2;
      *c3p++ += c3;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define AT(i,j) at[ (j)*(k) + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

#define min( x, y ) ( x < y ? x : y )

#define MB 48
#define KB 48

/* MY_MMult_inner prototype */
void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc );


/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, ib, pb;

  for ( p=0; p<k; p+=KB ){
    pb = min( k-p, KB );
    for ( i=0; i<m; i+=MB ){
      ib = min( m-i, MB );

      MY_MMult_inner( ib, n, pb, 
                      &A( i, p ), lda, 
                      &B( p, 0 ), ldb, 
                      &C( i, 0 ), ldc );
    }
  }
}


  


/* Inner kernel for computing C = A * B + C */

void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                          double *b, int ldb,
                                          double *c, int ldc )
{
  int i, j, p;
  double at[ MB * KB ];

  transpose_matrix( m, k, a, lda, at, k );

  for ( j=0; j<n; j+=4 ){
    double *c0p = &C(0,j), *c1p = &C(0,j+1), *c2p = &C(0,j+2), *c3p = &C(0,j+3);

    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &AT(0,i), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p++ ){
        register double aip = *ap;

        c0 = c0 + aip * *b0p++;
        c1 = c1 + aip * *b1p++;
        c2 = c2 + aip * *b2p++;
        c3 = c3 + aip * *b3p++;

        ap++;
      }
      
      *c0p++ += c0;
      *c1p++ += c1;
      *c2p++ += c2;
      *c3p++ += c3;
    }
  }
}

Compare also Optimization 8 (above) to this Optimization 15. (They are similar, except that the same optimizations are applied to Optization 7 and 14, respectively.)

Change the first lines in the makefile to 

  • OLD  := 8
NEW  := 15

  • make run

  • make plot

  • gv graph_8_vs_15.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

Optimizing Optimization 15

Again (as shown for Optimization 8), by transposing the matrix, it turns out that a larger block size can be used. (The reason for this is that the block of A can be increased beyond the size of the L1 cache.) Try changing MB = 144 and KB = 144 in file MMult15.c and reexecute make run, make plot, and gv graph_8_vs_15.eps. On my machine this yields

Optimization 16

Copy the contents of file MMult15.c into a file named MMult16.c and change the contents:

  • from

    to 

    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define AT(i,j) at[ (j)*(k) + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

#define min( x, y ) ( x < y ? x : y )

#define MB 48
#define KB 48

/* MY_MMult_inner prototype */
void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc );


/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, ib, pb;

  for ( p=0; p<k; p+=KB ){
    pb = min( k-p, KB );
    for ( i=0; i<m; i+=MB ){
      ib = min( m-i, MB );

      MY_MMult_inner( ib, n, pb, 
                      &A( i, p ), lda, 
                      &B( p, 0 ), ldb, 
                      &C( i, 0 ), ldc );
    }
  }
}


  


/* Inner kernel for computing C = A * B + C */

void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                          double *b, int ldb,
                                          double *c, int ldc )
{
  int i, j, p;
  double at[ MB * KB ];

  transpose_matrix( m, k, a, lda, at, k );

  for ( j=0; j<n; j+=4 ){
    double *c0p = &C(0,j), *c1p = &C(0,j+1), *c2p = &C(0,j+2), *c3p = &C(0,j+3);

    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &AT(0,i), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p++ ){
        register double aip = *ap;

        c0 = c0 + aip * *b0p++;
        c1 = c1 + aip * *b1p++;
        c2 = c2 + aip * *b2p++;
        c3 = c3 + aip * *b3p++;

        ap++;
      }
      
      *c0p++ += c0;
      *c1p++ += c1;
      *c2p++ += c2;
      *c3p++ += c3;
    }
  }
}


    /* Create macros so that the matrices are stored in column-major order */

#define A(i,j) a[ (j)*lda + (i) ]
#define AT(i,j) at[ (j)*(k) + (i) ]
#define B(i,j) b[ (j)*ldb + (i) ]
#define C(i,j) c[ (j)*ldc + (i) ]

#define min( x, y ) ( x < y ? x : y )

#define MB 144
#define KB 144

/* MY_MMult_inner prototype */
void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc );


/* Routine for computing C = A * B + C */

void MY_MMult( int m, int n, int k, double *a, int lda, 
                                    double *b, int ldb,
                                    double *c, int ldc )
{
  int i, p, ib, pb;

  for ( p=0; p<k; p+=KB ){
    pb = min( k-p, KB );
    for ( i=0; i<m; i+=MB ){
      ib = min( m-i, MB );

      MY_MMult_inner( ib, n, pb, 
                      &A( i, p ), lda, 
                      &B( p, 0 ), ldb, 
                      &C( i, 0 ), ldc );
    }
  }
}


  


/* Inner kernel for computing C = A * B + C */

void MY_MMult_inner( int m, int n, int k, double *a, int lda, 
                                          double *b, int ldb,
                                          double *c, int ldc )
{
  int i, j, p;
  double at[ MB * KB ];

  transpose_matrix( m, k, a, lda, at, k );

  for ( j=0; j<n; j+=4 ){
    double *c0p = &C(0,j), *c1p = &C(0,j+1), *c2p = &C(0,j+2), *c3p = &C(0,j+3);

    for ( i=0; i<m; i++ ){
      register double c0=0.0, c1=0.0, c2=0.0, c3=0.0;
      double *ap = &AT(0,i), *b0p=&B(0,j  ), *b1p=&B(0,j+1), 
                            *b2p=&B(0,j+2), *b3p=&B(0,j+3);

      for ( p=0; p<k; p+=4 ){
        register double aip;

        aip = *ap++;

        c0 = c0 + aip * b0p[0];
        c1 = c1 + aip * b1p[0];
        c2 = c2 + aip * b2p[0];
        c3 = c3 + aip * b3p[0];

        aip = *ap++;

        c0 = c0 + aip * b0p[1];
        c1 = c1 + aip * b1p[1];
        c2 = c2 + aip * b2p[1];
        c3 = c3 + aip * b3p[1];

        aip = *ap++;

        c0 = c0 + aip * b0p[2];
        c1 = c1 + aip * b1p[2];
        c2 = c2 + aip * b2p[2];
        c3 = c3 + aip * b3p[2];

        aip = *ap++;

        c0 = c0 + aip * b0p[3];
        c1 = c1 + aip * b1p[3];
        c2 = c2 + aip * b2p[3];
        c3 = c3 + aip * b3p[3];

        b0p+=4;
        b1p+=4;
        b2p+=4;
        b3p+=4;
      }
      
      *c0p++ += c0;
      *c1p++ += c1;
      *c2p++ += c2;
      *c3p++ += c3;
    }
  }
}

Compare also Optimization 9 (above) to this Optimization 16. (They are similar, except that the same optimizations are applied to Optization 8 and 15, respectively.)

Change the first lines in the makefile to 

  • OLD  := 9
NEW  := 16

  • make run

  • make plot

  • gv graph_9_vs_16.eps

This time the performance graph will look something like

The difference in performance can be attributed to the fact that ...

LinearAlgebraWiki: OptimizingGemm (last edited 2009-11-02 13:37:13 by RobertVanDeGeijn)