# Why isn't my Matrix-Vector Multiplication Scaling?

Sorry for the long post but I wanted to include everything that I thought was relevant in the first go.

# What I want

I am implementing a parallel version of Krylov Subspace Methods for Dense Matrices. Mainly GMRES, QMR and CG. I realized (after profiling) that my DGEMV routine was pathetic. So I decided to concentrate on that by isolating it. I have tried running it on a 12 core machine but the results below are for a 4 core Intel i3 Laptop. There is not much difference in the trend.

My KMP_AFFINITY=VERBOSE output is available here.

I wrote out a small code :

size_N = 15000
A = randomly_generated_dense_matrix(size_N,size_N); %Condition Number is not bad
b = randomly_generated_dense_vector(size_N);
for it=1:n_times %n_times I kept at 50
x = Matrix_Vector_Multi(A,b);
end

I believe this simulates the behaviour of CG for 50 iterations.

# What I have tried:

### Translation

I had originally written the code in Fortran. I translated it to C, MATLAB and Python (Numpy). Needless to say, MATLAB and Python were horrible. Surprisingly, C was better than FORTRAN by a second or two for the above values. Consistently.

### Profiling

I profiled my code to run and it ran for 46.075 seconds. This was when MKL_DYNAMIC was set to FALSE and all cores were used. If I used MKL_DYNAMIC as true, only (approx.) half the number of cores were in use at any given point of time. Here are a few details:

0x5cb51c        mulpd %xmm9, %xmm14     36.591s

The most time consuming process seems to be :

Call Stack                          LAX16_N4_Loop_M16gas_1
CPU Time by Utilization             157.926s
CPU Time:Total by Utilization       94.1%
Module                              libmkl_mc3.so

Here are a few pictures:

# Conclusions:

I'm a real beginner at profiling but I realize that the speed up is still not good. The sequential (1 Core) code finishes in 53 seconds. That is a speed up of less than 1.1 !

Real Question : What should I do to improve my speedup?

## Stuff that I think might help but I can't be sure:

• MPI (ScaLapack) implementation
• Manual Tuning (I don't know how. Please recommend a resource if you suggest this)

If anyone needs more (especially regarding memory) details, please let me know what I should run and how. I have never memory profiled before.

Your matrix is of size 15,000 x 15,000, so you have 225M elements in the matrix. This makes for roughly 2GB of memory. This is much more than the cache size of your processor, so it has to be loaded completely from main memory in every matrix multiplication, making for approximately 100GB of data transfers, plus what you need for the source and destination vectors.

The maximum memory bandwith of the i3 is approximately 21 GB/s based on the Intel specs, but if you look around the web, you'll find that at most half of that is really available in reality. Thus, at the very least, you'd expect your benchmark to last 10 seconds, and your actual measurement of 45 seconds isn't so far off that mark.

At the same time, you are also doing some 10 billion floating point multiplies and adds. Considering, say, 10 clock cycles for the combination, and 3 GHz clock rate, you'll come out at ~30 seconds. Of course, they can run concurrently with speculative memory loads if the cache is clever.

All in all, I'd say you're not too far off the mark. What would you have expected?

• Isn't there a way of getting atleast a speedup of 2-3? Mar 23, 2012 at 10:11
• @Nunoxic - You may want to benchmark memory performance on your system using a tool like SiSoftware Sandra. Wolfgangs analysis looks spot on to me, if your application is memory bandwidth bound, parallelisation is going to help little if at all. Also, look at any power saving options you may have, they may be throttling memory performance. Also, consider replacing your memory with higher quality memory, a lower CAS latency for example could make a big difference to your wall time. Mar 23, 2012 at 16:55

How are you doing the matrix-vector multiply? A double-loop by hand? Or calls to BLAS? If you're using MKL, I would strongly recommend using the BLAS routines of the threaded version.

Out of curiosity, you might also want to compile your own tuned version of ATLAS and see how that does on your problem.

Update

Following the discussion in the comments below, it turns out that your Intel Core i3-330M only has two "real" cores. The two missing cores are emulated with hyperthreading. Since in hyperthreaded cores both the memory bus and the floating-point units are shared, you won't get any speedup if any of the two are a limiting factor. In fact, using four cores will probably even slow things down.

What kind of results do you get on "only" two cores?

• I've tried ATLAs, GoTo and Netlib BLAS. All are weaker than MKL in performance. Is this expected or am I doing something wrong? I compiled ATLAS as mentioned in the handbook. Further, I have pasted my (exact) code here. Its calling MKL's BLAS. Mar 23, 2012 at 10:15
• Ok, and for the scaling, are you sure that in your baseline case, the code is only running on a single CPU? E.g. if you benchmark it, does the CPU usage histogram show only a single core? Mar 23, 2012 at 14:11
• Yes. The CPU histogram shows 1 core. Mar 23, 2012 at 14:59
• Just out of curiosity again, what do you get for two or three cores? Does your machine actually have four physical cores, or just two cores with hyperthreading? Mar 23, 2012 at 15:12
• How do I find that out? I havev included my KMP_AFFINITY in the main. Mar 23, 2012 at 15:14

I have the impression that row-major ordering is optimal for this problem with respect to memory access times, cache lines usage and TLB misses. I guess your FORTRAN version used column-major ordering instead, which could explain why it is consistently slower than the C version.

As already pointed out before, you're memory bandwidth limited here. What could possibly go wrong is that the vector $b$ isn't kept in the cache. You could test whether you observe the same (effective) memory bandwidth for size_N = 15000 than for size_N = 5000. If you do, there is a pretty good chance that the code is already optimal, and that the memory bandwidth of your system is simply not that great.

You could also test the speed if you just sum up all elements of the matrix in a single loop instead of the matrix vector multiplication. (You might want to unroll the loop by a factor 4, because non-associativity of addition could prevent the compiler from doing this optimization for you.)