Fortran matrix multiplication performance in different optimization

Question

I'm reading the book "Scientific Software Development with Fortran", and there is an exercise in it I think very interesting: "Create a Fortran module called MatrixMultiplyModule. Add three subroutines to it called LoopMatrixMultiply, IntrinsicMatrixMultiply, and MixMatrixMultiply. Each routine should take two real matrices as argument, perform a matrix multiplication, and return the result via a third argument. LoopMatrixMultiply should be written entirely with do loops, and no array operations or intrinsic procedures; IntrinsicMatrixMultiply should be written utilizing the matmul intrinsic function; and MixMatrixMultiply should be written using some do loops and the intrinsic function dot_product. Write a small program to test the performance of these three different ways of performing the matrix multiplication for different sizes of matrices." I did some test of multiply of two rank 2 matrix and here are the results, under different optimization flags: ``` compiler:ifort version 13…

Accepted Answer

There are several things one should be aware of when looping over array elements:

Make sure the inner loop is over consecutive elements in memory. In your current 'loop' algorithm, the inner loop is over index k. Since matrices are laid out in memory as columns (first index varying most rapidly when going through the memory), accessing a new value of k might need to load a new page into cache. In this case, you could optimize your algorithm by reordering the loops as:

```
do j=1,n
    do k=1,size(mtx1,dim=2)
        do i=1,m
            mtx3(i,j)=mtx3(i,j)+mtx1(i,k)*mtx2(k,j)
        end do      
    end do
end do
```

now, the inner loop is over consecutive elements in memory (the `mtx2(k,j)` value will be probably be fetched only once before the inner loop by the compiler, if not you can store it in a temporary variable before the loop)

Make sure the entire loops can fit into the cache in order to avoid too much cache misses. This can be done by blocking the algorithm. In this case, a solution could be e.g.:

```
l=size(mtx1,dim=2)
ichunk=512 ! I have a 3MB cache size (real*4)
do jj=1,n,ichunk
    do kk=1,l,ichunk

do j=jj,min(jj+ichunk-1,n)
    do k=kk,min(kk+ichunk-1,l)
        do i=1,m
            mtx3(i,j)=mtx3(i,j)+mtx1(i,k)*mtx2(k,j)
        end do      
    end do
end do

end do
end do
```

in which case performance will depend in the size of `ichunk`, especially for large enough matrices (you could even block the inner loop, this is just an example).

Make sure the work needed to perform the loop is much smaller than the work inside the loop. This can be solved by 'loop unrolling', i.e. combining several statements in one iteration of the loop. Usually the compiler can do this by supplying the flag `-funroll-loops`.

If I use the above code and compile with the flags `-O3 -funroll-loops`, I get a slightly better performance than with `matmul`.

The important thing to remember of those three is the first point about loop ordering, since this is something that will affect performance in other use cases, and the compiler cannot usually fix that. The loop unrolling, you can leave to the compiler (but test it, as this does not always increase performance). As for the second point, since this is dependent on the hardware, you shouldn't (generally) try to implement a very efficient matrix multiplication yourself and instead consider using a library such as e.g. atlas, which can optimize for cache size, or a vendor library such as MKL or ACML.

Fortran matrix multiplication performance in different optimization

Solution

Problem