Python matrix provide with numpy.dot()

Question

During my acquaintance with CUDA in Python (numba lib), I implemented matrix provide methods:

• Just with numpy.dot()
• Strassen algorithm with numpy.dot()
• Blocks method on GPU
• Strassen algorithm on GPU

So I tested it on 2 types of data:

• numpy.random.randint(0, 5, (N, N)) # with int32 elements
• numpy.random.random((N, N)) # with float64 elements

For int32 i obtained expected result, where my GPU algroithms performed better than CPU with numpy:

However, on float64 type, numpy.dot() outperformed all my GPU methods:

So, question is: Why is numpy.dot() so fast with float64 arrays, and does numpy use the GPU?

Answer 1

A typical installation of numpy will be dynamically linked against a BLAS library, which provides routines for matrix-matrix and matrix-vector multiplication. For example, when you use np.dot() on a pair of float64 arrays, numpy will call the BLAS dgemm routine in the background. Although these library functions run on the CPU rather than the GPU, they are often multithreaded, and are very finely tuned for performance. A good BLAS implementation, such as MKL or OpenBLAS, will probably be hard to beat in terms of performance, even on the GPU*.

However, BLAS only supports floating point types. If you call np.dot() on integer arrays, numpy will fall back on using a very simple internal C++ implementation, which is single-threaded and much slower than a BLAS dot on two floating point arrays.

Without knowing more about how you conducted those benchmarks, I would bet that a plain call to numpy.dot would also comfortably beat your other 3 methods for float32, complex64 and complex128 arrays, which are the other 3 types supported by BLAS.

---

* One possible way to beat standard BLAS would be to use cuBLAS, which is a BLAS implementation that will run on an NVIDIA GPU. The scikit-cuda library seems to provide Python bindings for it, although I've never used it myself.