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 -mklwith
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,
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.