use cublas*dgmm in linalg.dot about scikit-cuda

[untom]

FYI, I looked into this:

using dgmm is noticably faster: on my machine the following code takes 70ms with the current implementation, while an implementation based on *dgmm takes 20ms.

Test code:

from scikits.cuda import linalg as la
import pycuda.autoinit
import pycuda.gpuarray as gpu
import numpy as np
la.init()
n = 1024*4
x = np.asarray(np.random.rand(n, n), np.float32)
d = np.asarray(np.random.rand(n, 1), np.float32).reshape(-1)
xd = gpu.to_gpu(x)
dd = gpu.to_gpu(d)
%time la.dot_diag(dd, xd).get()

However, the current version allows for incompatible types between A and D. E.g. it is possible to multiply a np.complex128 vector with a np.float64 diagonal matrix. *dgmm requires that the diagonal matrix as well as the dense matrix be of the same type.

The current behaviour required in the linalg.pinv implementation. However pinv could certainly implemented using CULA's gelss, which is probably a faster solution anyway. Then, that requirement would go away.

Alternatively, one could use *dgemm if both arguments are of the same type, and use the old code if they are not (personally, I'd prefer the "cleaner" version of requiring D and A be of the same type).

As a sidenote, I think the cublas API is a bit cleaner than the scikit-cuda one: For me, it's more natural to specify whether D is multiplied from the Left/Right (= how cublas does it) than to specify whether A is transposed or not.

With that said, I can send a PR for my current code which replaces diag_dot with calls to *dgmm if you're interested (But I guess I'd have to change pinv to use gelss or switch to the old implementation if A and D have different types).

[lebedov]

A dgmm-based dot_diag() would be welcome. Requiring that all of the arguments have the same type seems reasonable since maintaining a single numerical storage type throughout a CUDA-based program is usually a given.

I would prefer to keep the current transpose-based interface to dot() and dot_diag() for the time being because some other code that my adviser's research group has written depends on it.

It would be interesting to see whether a gelss-based pinv() is be noticeably faster than the gesvd-based one irrespective of the proposed dot_diag() improvement (although I would have to investigate whether there are meaningful numerical differences between the two algorithms when applied to poorly conditioned matrices that might make one preferable to the other in certain circumstances).