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:
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.
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
withMKL_DYNAMIC=TRUE
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.