Netlib BLAS, ACML BLAS and Goto BLAS, which one is faster?

[neodreamer]

I tested the following simple matrix multiplication program on my laptop, Inspiron 1501, the operating system Ubuntu Linux, the program did matrix multiplication C=AxB for 1000 times, A,B,C are all 300x300 matrix. Use time command, I get the following result,

Netlib BLAS:
real 0m2.150s
user 0m2.128s
sys 0m0.016s

ACML BLAS:
real 1m24.560s
user 1m22.985s
sys 0m1.556s


Goto BLAS:
real 1m22.421s
user 1m21.025s
sys 0m0.116s


Why Netlib BLAS are much much faster than those so-called optimized BLAS? I tested the same program under different operating systems and on different computers and got the similar results: plain Netlib BLAS is much faster. Can anybody tell me why?



CPP program
/////////////////////////////////////////////////////////////////////////////////////////
#include<iostream>
#include<complex>


// define complex type
typedef std::complex<double> zomplex;

const zomplex Im = zomplex(0.0, 1.0);

/* ZGEMM */
extern "C" void zgemm_(const char*, const char*, const int*, const int*, const int*,
const double*, void*, const int*, void*, const int*, const double*, void*, const int*);


const int m = 300; // the dimension of the matrices in the test
const int n = m*m; // the number of elements in a matrix
const int ntimes = 1000; // the number of times we do the computation C = AxB;


int main(void)
{
zomplex* A = new zomplex[m*m];
zomplex* B = new zomplex[m*m];
zomplex* C = new zomplex[m*m];


const char trans[] = {'N'};
const double alpha[] = {1.0, 0.0};
const double beta[] = {0.0, 0.0};

int i;

for (i = 0; i < ntimes; ++i)
{
zgemm_(trans, trans, &m, &m, &m, alpha, A, &m, B, &m, beta, C, &m);
}


delete[] A;
delete[] B;
delete[] C;
}

//////////////////////////////////////////////////////////////////////////////////////////////////////


Makefile
//////////////////////////////

CXX = g++
CC = g++
CXXFLAGS = -Wall -W
#LOADLIBES= -lf77blas -latlas -lg2c
LOADLIBES= -lgotoblas -lg2c -lpthread
#LOADLIBES= -lblas -lg2c
#LOADLIBES= -lacml -lg2c

TARGETS = zgemm_test

all: $(TARGETS)


.PHONY: clean
clean:
rm -rf *.o *~ *orig *exe $(TARGETS)


////////////////////////////////////////////////////////

[buttari]

Neo,
the answer is actually quite easy even if may not be obvious. The reason for the strange behavior you're seeing is that you're multiplying two empty matrices, i.e., two matrices full of zero values. Compilers usually generate custom code for these cases and, as a result, no floating point operation is actually performed. In fact, if you make some quick computations, it turns out that you progam, when linked to reference BLAS, is spinning at 200 Gflop/s on a processor can can barely do 2 Gflop/s.
Now if, instead you fill your matrices with random values, here's what you get (this is what I measured on my laptop):

Netlib BLAS
real 1m1.237s
user 0m58.040s
sys 0m0.068s

MKL (should be the same with Goto, ACML...)
real 0m7.220s
user 0m6.952s
sys 0m0.024s

Note that I set the number of repetitions to 100.
Here's how I modified your code:

Code: Select all

#include<iostream>
#include<complex>


// define complex type
typedef std::complex<double> zomplex;

const zomplex Im = zomplex(0.0, 1.0);

/* ZGEMM */
extern "C" void zgemm_(const char*, const char*, const int*, const int*, const int*,
const double*, void*, const int*, void*, const int*, const double*, void*, const int*);


const int m = 300; // the dimension of the matrices in the test
const int n = m*m; // the number of elements in a matrix
const int ntimes = 100; // the number of times we do the computation C = AxB;


int main(void)
{
zomplex* A = new zomplex[m*m];
zomplex* B = new zomplex[m*m];
zomplex* C = new zomplex[m*m];


int i;
 for(i=0; i<n; i++){
   A[i]=zomplex((double) (rand()-1)/RAND_MAX, (double) (rand()-1)/RAND_MAX);
   B[i]=zomplex((double) (rand()-1)/RAND_MAX, (double) (rand()-1)/RAND_MAX);
   C[i]=zomplex((double) (rand()-1)/RAND_MAX, (double) (rand()-1)/RAND_MAX);
 }

const char trans[] = {'N'};
const double alpha[] = {1.4, 0.5};
const double beta[] = {2.3, 0.2};


for (i = 0; i < ntimes; ++i)
{
zgemm_(trans, trans, &m, &m, &m, alpha, A, &m, B, &m, beta, C, &m);
}


delete[] A;
delete[] B;
delete[] C;
}


[neodreamer]

Thank you very much for you detailed explanation. But I have another question.


I changed the following constants in your code:

m = 10;
ntimes = 1000000;

and I got the following results;


Netlib BLAS
real 0m8.820s
user 0m8.813s
sys 0m0.004s


ACML BLAS
real 0m11.654s
user 0m11.621s
sys 0m0.032s

Goto BLAS
real 0m6.246s
user 0m6.232s
sys 0m0.012s


When m = 5, Netlib BLAS is a little faster than Goto BLAS, ACML is the slowest. Is there any reason for that?

[buttari]

Thank you very much for you detailed explanation. But I have another question.


I changed the following constants in your code:

m = 10;
ntimes = 1000000;

and I got the following results;


Netlib BLAS
real 0m8.820s
user 0m8.813s
sys 0m0.004s


ACML BLAS
real 0m11.654s
user 0m11.621s
sys 0m0.032s

Goto BLAS
real 0m6.246s
user 0m6.232s
sys 0m0.012s


When m = 5, Netlib BLAS is a little faster than Goto BLAS, ACML is the slowest. Is there any reason for that?

Neo,
yes, the problem is that, repeating the operation 10000 time on such a small data doesn't provide any useful information. A 5x5 matrix is just 2KB and, thus, it entirely fits in L1 cache (in some architectures a 5x5 matrix may even fit in registers!). All the highly optimized BLAS libraries like Goto, ACML etc mostly try to optimize cache behavior; thus, in the case of such a small matrix there's nothing left to optimize because it's already an ideal case.
Regards

Alfredo

[PHinker]

I would like to also suggest using the Sun Performance Library (aka Perflib or sunperf) which is included free with the Sun Studio compilers. Taking the examples in this thread and running on a 2.2 Ghz opteron which was running Solaris 10 I got the following :

% CC -o driver driver.cc -library=sunperf -m64
% time driver
58.0u 0.0s 0:58 99% 0+0k 0+0io 0pf+0w

% CC -o driver2 driver2.cc -library=sunperf -m64
% time driver2
5.0u 0.0s 0:05 84% 0+0k 0+0io 0pf+0w

Most of the optimized versions run at very high percentages of machine peak for the xGEMM routines (typically 88-92% of peak). Goto, MKL, ACML, and Sunperf are typically within 2-4% of each other. ATLAS comes in a couple % below them and then the reference (netlib) implementation is usually around 40-50% of peak if you use a very good optimizing compiler *and* the problems are small enough to fit in L2 cache.

buttari is correct in saying that the original timer isn't giving a good indication of the performance of the various libraries. When you set alpha = (1.0, 0.0) and beta = (0.0, 0.0) and not initializing the A and B matrices, you're likely multiplying 0.0 * 0.0. If you look in the netlib implementation of the xGEMM routines, you'll see something like the following :

Code: Select all

                  DO 80 L = 1,K
                      IF (B(L,J).NE.ZERO) THEN
                          TEMP = ALPHA*B(L,J)
                          DO 70 I = 1,M
                              C(I,J) = C(I,J) + TEMP*A(I,L)
   70                     CONTINUE
                      END IF
   80             CONTINUE


This allows you to skip alot of operations when the matrices have zeros but it's an 'optimization' that destroys performance if there are few zeros in the matrices.

Also, the optimized versions of the routines do alot of work to insure that portions of the A and B matrices are in cache. This does add some overhead which can hurt performance on very small matrix sizes (like your example where m = 5).

In the spirit of full disclosure, I work on the Sun Performance Library. One advantage the Sun Performance Library has is that there are optimized versions that run on SPARC and X64 architectures(both intel and Amd). Perflib also runs on both Solaris and Linux. MKL runs very well on Intel processors, ACML runs very well on the AMD processors, ATLAS runs okay as long as it's built on the processor where it will run.

Each of the implementations has strengths and weaknesses. It's not impossible to find particular examples of where one library outperforms another. It depends on how your application calls the underlying library. If your data contains lots of zeros and you are using small matrices, it might very well mean that your best solution is the reference implementation.[/code]