# Are DAXPY, DCOPY, DSCAL overkills?

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.

• What hardware are you running on? – Pedro Apr 15 '12 at 14:11
• Don't assume there is nothing faster than BLAS/LAPACK. They are optimized for utility, not necessarily gold-medal-winning speed. When speed is your need, you might try this. – Mike Dunlavey Sep 4 '12 at 12:49
• DCOPY(n,2.89*x,1,y,1) is not going to do what you want. It's straight up wrong. The function you want is DAXPY. – Jeff Hammond Jan 25 '14 at 21:26
• MKL_DYNAMIC=TRUE is terrible for performance. I know of no scientific code that benefits from this. Turn it off and set the thread number via MKL_NUM_THREADS/OMP_NUM_THREADS and turn on OMP_SCHEDULE=STATIC. – Jeff Hammond Jan 25 '14 at 21:35

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.

• While I cannot logically refute (2), I don't believe that it is pertinent here. DCOPY, DSCAL and DAXPY are almost trivial to implement properly, hence I doubt people make mistakes. The issue is that their triviality derives from the fact that these functions hit the hardware limit very quickly and thus very few optimizations are effective. – Jeff Hammond Jan 25 '14 at 21:28

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.

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.
• Have you run the experiment where you eliminate the conditionals in DAXPY to see if that has any impact on performance? I seriously doubt it does. – Jeff Hammond Jan 25 '14 at 21:32
• No, but I have written pure assembly code and outperformed vendor-provided DAXPY in BLAS on a couple platforms :) – Aron Ahmadia Jan 27 '14 at 0:49

BLAS1 functions represent a set of kernels that are bandwidth limited because their compute intensity is low. In particular, these kernels do O(1) flops per memory access. This means that on modern hardware, they run at a small fraction of peak and there is essentially nothing you can do about. The best implementation of BLAS1 will check for alignment and modulo the FPU vector length and perform at bandwidth peak, which is likely to be 5-10% of compute peak.

When you write these operations explicitly in the source, a good compiler immediately recognizes them and inlines some optimal implementation that is equivalent to the BLAS1 one noted above. However, because the compiler knows more about the context, it can avoid certain branches (not that these matter that much) and function call overhead, as well as potentially perform higher-order transformations in the code that would be blocked by a function call to an opaque library.

There are a variety of experiments you can perform to determine what is actually affecting the performance of your code. They are pretty obvious so I won't list them here.