nVidia, for example, has CUBLAS, which promises 7-14x speedup. Naively, this is nowhere near the theoretical throughput of any of nVidia's GPU cards. What are the challenges in speeding up linear algebra on GPUs, and are there faster linear algebra routings already available?
I can't answer the second half of your question as far as other implementations out there but I can provide some insight as to the challenges. For reference, I personally used ViennaCL on a nVidia GTX 560 Ti with 2GB of memory for my benchmarks.
Over serial code on a mid-range i5, I saw speed-ups for dense matrix multiplications of approximately 40x. For operations such as a vector-scalar multiply I saw as high as 1000x speed-ups. The 800 pound gorilla in the room, however, is memory bandwidth. For most commercial GPUs, you will be using something like PCIe which limits you to about 6GB/s of throughput. In my case, while the computation was 40x faster, the three matrix copies (two to the GPU, and one back) each took about as much time as just doing the computation on the CPU.
The problem then with any general library for GPU linear algebra is going to be that they can't really re-use objects on the GPU, because they don't know what you are going to do with them. So every call to a compute kernel will likely require copying to the GPU, then copying the result back. This will eat up a large portion of the gains.
If you can reuse objects such as matrices, then you could write the higher level algorithms to avoid as much memory management as possible, but a library would be hard pressed to do this efficiently.
I hope that this helps, and I am sure there are other people here who are much more experienced in this, but these are the experiences and impressions I got during my short foray into GPU computing.
Let me focus only on CUDA and BLAS.
Speedup over an host BLAS implementation is not a good metric to assess throughput, since it depends on too many factors, although I agree that speedup is usually what one cares about.
If you look at the benchmarks published by NVIDIA and take into account that the Tesla M2090 has 1331 Gigaflops (single precision) and 665 Gigaflops (double prec.) peak performance, you will see that for SGEMM and DGEMM we have a measured throughput nearly at 60% of the theoretical one, which is pretty good.
But how do you define the flops performance? Flop count / elapsed time, where flop count is $2\,mnk$ ($m\times k$ and $k\times n$ are the matrix dimensions), and elapsed time can include or not the transfer time from host to GPU memory and back. (As Godric correctly points out, this make a big difference!)
As what regards sustained floating point throughput, I think that flops should be computed without taking into account data and result transfer times, and this makes speedup comparisons difficult. Furthermore you have to take into account the matrix size, since best performance is for big matrices.
Bottom line: speedup of a real life application can be very different from peak measured performance on linear algebra routines, since you have to take into account GPU initialization, data transfer times, etc. etc.
So I won't answer your question about the fastest library, since the question makes no sense unless a precise metric and problem is defined. All this said, i think that cuBLAS and MAGMA are a very good starting point.