Are DAXPY, DCOPY, DSCAL overkills?

Question

I have implemented CG in FORTRAN by linking it to Intel MKL.

When there are statements like: (Refer Wikipedia)

 p=r; 
 x=x+alpha*p
 r=r-alpha*Ap;

or similar ones in QMR (in much greater quantity)

v_tld = r;
y = v_tld;
rho = norm( y );
w_tld = r;
z = w_tld;
xi = norm( z ); (and more)

Does it make sense to use BLAS Level 1 implementations such as DAXPY, DCOPY, DSCAL? The motivation for my question is:

1. I have 2 implementations of the algorithms. One wherein I have only linked Norms and MatVecs to MKL; copying, scaling and adding is done by Fortran's intrinsic functions and another where every possible subroutine is carried out by BLAS.

2. I was of the notion that nothing can get faster than BLAS. But, it turns out that my code using Fortran's intrinsic functions ran 100% faster than one with BLAS Level 1 subroutines (FWIW, This wasn't a small problem, it was solving a dense system of size 13k x 13k which filled up my 4 GB RAM). I was running both on 2 threads (on a 2 core machine) ifort QMR.f90 -mkl with MKL_DYNAMIC=TRUE

3. I had asked a question on SO regarding the extension of BLAS but as I tried to include BLAS Level 1 into my code, my code kept getting slower and slower.

Am I doing something wrong or is this expected?

Also, Does it make sense to try extend BLAS to do non-obvious operations like y = 2.89*x by DCOPY(n,2.89*x,1,y,1) or even DSCAL then DCOPY?

---

What is also interesting is, DDOT and DNRM2 improve performance. I attributed it to the fact that since they carry out double precision multiplications, putting them in parallel might help.

Supplementary Question : When do you decide whether a BLAS Level 1 operation is actually going to aid the performance?

Adding : Currently, I am running on a i3 2.13 GHz Laptop with 4 GB RAM and Debian 64 bit Proc info here. But, I get similar answers on an Intel Xeon 12 core Workstation with 24 GB RAM.

Answer 1

If your goal is really to squeeze as much performance out as possible, then it is important to remember:

1. The (BLAS) library might not have been tuned for your exact system/configuration.
2. Library developers make mistakes.

A vendor-tuned BLAS library should certainly be your default approach, but if you have taken the time to time individual kernels and noticed that some other implementation is faster, then, by all means, use the other implementation. Missing out on the usage of vector intrinsics could possibly lead to the large performance difference.

It is possible that your best bet for simple routines like daxpy and dscal is a hand-written loop which exploits vector intrinsics.

Answer 2

Given the state of optimizing compilers now a days, I don't think there's much voodoo in the linear BLAS routines, e.g. DAXPY, DCOPY, and DSCAL, that your compiler won't do already, e.g. SSE-vectorization and loop unrolling.

If the code is the same, the only difference between your routine and a call to MKL's BLAS is the overhead of the function call and whatever extra magic MKL might be trying to do in there. If this is the case, the difference between your code and MKL's code should be a constant, independent of the problem/vector size.

This question has interesting echoes of this question, which also uses DAXPY as an example.

Answer 3

The BLAS standard actually has several checks for correctness of the function arguments that are unnecessary in many situations. See this reference implementation of daxpy.f. Additionally, constants like INCX are usually known to you at compile-time, but may not be assumed by the implementation. BLAS calls cross compilation units, and I am not aware of any compilers that will be able to optimize these out without turning on whole program optimization.

• A funny sidenote to this is that the Intel Compiler now recognizes BLAS 3 matrix-matrix multiply loops, and will transform this code into the equivalent xgemm call, with enough optimization enabled.