Fix bug in linalg.diag

`linalg.diag` will corrupt memory (and cublas state) when extracting the diagonal of matrices that are larger than 32x32.  To see this, try:

        v = np.array(np.random.rand(64, 64), np.float32)
        v_gpu = gpuarray.to_gpu(v)
        d_gpu = linalg.diag(v_gpu)
        assert np.all(np.diag(v) == d_gpu.get())

The assert will fail, and all subsequent calls to PyCUDA stuff will also fail.  I'm not entirely sure why the original code worked, since -- at least to my understanding, the elementwise kernel used for the copy will start a thread for each element of the large matrix, and try to write an element in the diagonal. (So it has way to many writes that should fall out of bounds).

Anyhow, this PR should fix the bug. As additional benefit, it uses a cublas function so it doesn't need to go through a pycuda-compile cycle on its first call.