Earlier quoted context omitted.
Are you sure that BLAS uses different matrix multiplication algorithm ? I havent looked at the actual code in a while, but from what I have seen they use the standard vanilla matrix multiplication algorithm whose complexity is no less than O(N^3) (for square matrices). They achieve the speed by exploiting the deep cache structure of modern machines. Essentialy by blocking (sometimes unrolling) loops and consuming dat…
BLAS is an API. Just because the reference implementation does things in a particular way does not necessarily mean that that your vendor implementation cannot do things differently (e.g. Apple's numerics team has spent a considerable amount of time tuning BLAS for Apple hardware). For the particular problem discussed in the original article, there are at least two ways the multiplication A'A could potentially be mad…
I would also hazard a guess that the tuning that you mention does not involve a change in the algorithm but are essentially reordering the steps of the algorithm to obtain better caching. I was responding to the parent post which conjectured that BLAS implementations use different algorithms. If you look at your two suggestions, none of them actually change the complexity class of the number of floating operations, but wall clock time oh absolutely.
Although there are matrix multiplication algorithms that have a complexity less than O(N^3) the constants for these are so large enough that the sizes of matrices for which there will be any appreciable benefit are extremely rare to come by.