Why isn't my Matrix-Vector Multiplication Scaling?

Question

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:

Address Line    Assembly                CPU Time

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%
Overhead Time                       0us
Overhead Time:Total                 0.0%    
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:

• Pthreads implementation
• 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.

Answer 1

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?

Answer 2

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?

Answer 3

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.)