Mixing PyCULA - Expert

Within the mixed context environment created by mixed_init(), PyCULA allows you freely mix CULA LAPACK calls with your own calls to cuBLAS, cuFFT, (cuSP is coming soon!), and any other custom PyCUDA kernels you can cook up. The magic all occurs behind the scenes, but we only claim it will work with CULA if you use mixed_init() and cula_gpuarrays. Below are two examples using just this technique.

As a note: Right now we don’t support multiple devices or other black magic, but neither does anyone else, so you aren’t missing much. As other packages move forward, we will add such features.

CUDA Libraries: cuBLAS, cuFFT,cuSP

In the dawn of CUDA time, we began to wrap various libraries, including skeletons of cuBLAS and cudart. It turns out that a fellow python~cuda developer, Lev Givon, has provided a python package called scikits.cuda that includes wrappers for cuBLAS, and cuFFT. There are even some wrappers for CULA, but they are rather bare, and we recommend you use PyCULA for this aspect. The scikits.cuda cuBLAS wrappers are very nice, and here is a very simple example of how to use them:

# import PyCULA module
from PyCULA.cula import *

# import scki
import scikits.cuda.cublas as cublas

#initialize special mixed environment;

#make a numpy array; you may use float32 or float64 dtypes
a = numpy.array([[1,1],[0,1]], dtype=numpy.float32)
print a

#make a cula_Fpitched_gpuarray on gpu device like a
a_ = cula_Fpitched_gpuarray_like(a)

#note that a_ is transposed now, and has a pitch of 48.
print a_

#run PyCULA gpu_devdev* routine, this keeps results on device; print results
#PyCULA checks for pitch and performs all the necessary work behind the scenes.
golden_ = gpu_devdevsvd(a_)
print golden_

#Note that golden_ is a cula_gpuarray instance, but can be used exactly in place of PyCUDA.gpuarray

## run scikits.cuda.cublas's  cublasSnrm2 routine to find norm:
norm = cublas.cublasSnrm2(golden_.size, golden_.gpudata, 1)
print norm

Similarly, you can mix with other CUDA runtime libraries using scikits.cuda wrappers to your hearts content. Should you wish to use a different runtime libs wrapper, or your own wrappers... if you can make the code work with PyCUDA.GPUArrays then it should work just as well with cula_gpuarrays allowing you to mix with calls to CULA LAPACK by the scheme above.

PyCUDA Custom Kernels

Thats all fine and dandy, but you have your own kernels that you wish to run? Well, we can use PyCUDA to run your custom kernels. This is exactly how the PyCULA project started... we had some graph generators and wished to find the spectrums of such graphs quickly. Noticing that we were spending a lot of time making (and moving) large matrixes we decided it would be faster to ship a much smaller edge-list and use the GPU to expand it into a matrix on device in parallel.

Writing a Kernel Source Module

First we need to write our custom kernel. More information about this can be found in the PyCUDA documentation, but it basically works like this:

  1. Write some C code for your kernel
  2. Use PyCUDA.compiler.SourceModule to create a SourceModule
  3. Use ‘get_function’ to create a python callable function
  4. Call your CUDA-C from within python

This example is slightly complicated, but it is really useful! It shows how to use PyCULA to improve a real world computing problem. Within the code you can also see places where parameters can be tuned for performance. This example (and all the others) is included in the PyCULA/examples folder.

import numpy
from collections import defaultdict

from pycuda.compiler import SourceModule
from PyCULA.cula import *

# verbose printing for debugging etc
verbose = None
tuning = None
def vprint(s):
    global verbose
    if verbose:
        print s
def tprint(s):
    global tuning
    if tuning:
    print s

# cpu_expand_kernel
# Input edgelist as numpy.array, datatype
# For testing/comparing to host cpu
# Output discrete Laplacian as numpy.array in host memory

def cpu_expand_kernel(eList_nx2):
    #find size of nxn matrix
    n =  max(eList_nx2.flatten())
    vprint('.... expanding edgelist on cpu....')
    for row in eList_nx2:
        [i,j] = row
        #fix indexing to start at 1
        [i,j] = [i-1,j-1]
        M[i,j] = 1
        M[j,i] = 1
        if i==j:
            M[i,i] = M[i,i] - 1
            M[i,i] = M[i,i] - 1
            M[j,j] = M[j,j] - 1
    return M

# gpu_expand_kernel
# Input edgelist as numpy.array, datatype
# For fast GPU expansion
# Output discrete Laplacian as gpuarray in device memory

def gpu_expand_kernel(edges):
    def degrees(edges):
        h = defaultdict(int)
        for e in edges:
            i, j = e
            h[i] += 1
            h[j] += 1
        return [d for (v,d) in sorted(h.items())]
    D = degrees(edges)
    n =  len(D)
    m = len(edges)
    nsquare = int(n*n)
    vprint('.... expanding edgelist on gpu....')
    # host to device
    # Allocate n^2 zero array on GPU device

    M_gpu = cula_zeros(nsquare,dt)

    mod_expand = SourceModule("""
        __global__ void offDiagonal( float *M_gpu, int n, float *edges, int q, int blocksize1)
            int edge = threadIdx.x + (blocksize1*blockIdx.x) + q;
            int i = (int) edges[2*edge] - 1;
            int j = (int) edges[2*edge + 1] - 1;
            M_gpu[i * n + j] = 1.0;
            M_gpu[j * n + i] = 1.0;

        __global__ void onDiagonal(float *M_gpu, int n, float *degrees, int p,int blocksize2)
            int vertex = threadIdx.x + (blocksize2*blockIdx.x) + p;
            M_gpu[vertex*n + vertex] = degrees[vertex];

    mod_expand_64 = SourceModule("""
        __global__ void offDiagonal( double *M_gpu, int n, double *edges, int q, int blocksize1)
            int edge = threadIdx.x + (blocksize1*blockIdx.x) + q;
            int i = (int) edges[2*edge] - 1;
            int j = (int) edges[2*edge + 1] - 1;
            M_gpu[i * n + j] = 1.0;
            M_gpu[j * n + i] = 1.0;

        __global__ void onDiagonal(double *M_gpu, int n, double *degrees, int p,int blocksize2)
            int vertex = threadIdx.x + (blocksize2*blockIdx.x) + p;
            M_gpu[vertex*n + vertex] = degrees[vertex];

    #Define GPU C Kernels in pyCUDA
    if dt==numpy.float32:
        offDiagonal = mod_expand.get_function("offDiagonal")
        onDiagonal = mod_expand.get_function("onDiagonal")
    elif dt==numpy.float64:
        offDiagonal = mod_expand_64.get_function("offDiagonal")
        onDiagonal = mod_expand_64.get_function("onDiagonal")

    #Manual Block index tuning for Kernel1
    if m <= blocksize1:
        #offdiagonal kernel call small
        #offdiagonal kernel call full blocks
    if d !=0:
        #offdiagonal kernel call partial blocks

    del edges_gpu

    #Manual Block Index tuning for Kernel2
    if n <= blocksize2:
        # ondiagonal kernel call small
        # ondiagonal kernel call, full blocks
        if b != 0:
        # ondiagonal kernel call, partial blocks

    del D_gpu

    # tprint for tuning
    tprint('length of edgelist')
    tprint('number of blocks kernel1')
    tprint('partial block length kernel1')
    tprint('length of "degree" list')
    tprint('number of whole blocks kernel2')
    tprint('partial block length kernel2')

    #reshape M_gpu to nxn
    M_gpu.shape = n,n

    return M_gpu

pheww, that was long. Your functions need not be this complicated, but I wanted to provide an example that really improved performance over the best known alternatives(read below). The next step is to call this module then follow with some CULA linear algebra, and fortunately this is quite easy from here on out!

Calling Kernel and Mixing With CULA

So here the steps are simpler:

  1. We need to import our tools, including our kernel_modules created above
  2. Run mixed_init() to create a proper mixed enviornment
  3. Make an input data set (provided)
  4. Call our Custom Kernel
  5. Call CULA syev to find the eigenvalues of the resulting matrix

As code it looks like:

from PyCULA.cula import *
import numpy as np
import kernel_modules                       #import custom kernel module code


# This is an edgelist, it represents a reduced intermediary form of a larger sparse matrix.
E=np.array([[1 , 2 ],
[ 1 , 16 ],
[ 2 , 3],
[3 , 4 ],
[ 4 , 5 ],
[ 5 , 6 ],
[ 6 , 7 ],
[ 7 , 8 ],
[ 8 , 9 ],
[ 8 , 5 ],
[ 9 , 10 ],
[ 9 , 4 ],
[ 10 , 11 ],
[ 11 , 12 ],
[ 12 , 13 ],
[ 13 , 14 ],
[ 13 , 10 ],
[ 14 , 15 ],
[ 14 , 9 ],
[ 15 , 16 ],
[ 15 , 4 ],
[ 16 , 17 ],
[ 16 , 3 ],
[ 17 , 18 ],
[ 18 , 19 ],
[ 19 , 20 ],
[ 19 , 16 ],
[ 20 , 1 ]], dtype=np.float64)

print E

# expand reduced form into full matrix using GPU kernel
# holds array M on GPU
M = kernel_modules.gpu_expand_kernel(E)

print M

eigvals = gpu_devsyev(M)            # run CULA syev routine on GPU array M
print eigvals                               # print solutions returned by CULA as numpy array on host

When I claimed this is an improvement, I really meant it. To chug along a similar problem on a 10kx10k or larger graph in MATLAB may take anywhere from 10-40minutes using a fast machine. Just using the GPU based CULA for matrix computations yielded times in the 5-10 minute range on large problems with good hardware. Using this reduced intermediary form and GPU expansion technique above, the some computations can be had in 1 or 2 minutes. That is a 5-40x speedup depending on precision and problem size.

Best of all, is that since you are developing in Python, its easy to input data, store your results, and share your science with the world.

Thanks Again!

Thanks for checking out our tutorial documentation. Should you have more questions, the docstrings from the source code follow in the next section. Here you will find more detailed information about what particular wrappers are expecting as input and will be returning as output. There are MANY functions, so use search to speed your query, and don’t forget that docstring magic allows you to get help interactively in the interpreter:

>>>from PyCULA.cula import *

    Help on function gpu_svd in module PyCULA.cula:

    gpu_svd(A, left_vectors=False, right_vectors=False)
        Takes numpy.array A and computes the SVD, returns Singular Values (optionally left/right vectors) as np.arrays.

        Keyword arguments:
        A -- input matrix
        left_vectors -- default == False; computes left singular vectors
        right_vectors -- default == False; computes right singular vectors

        Refer to CULA docs regarding gesvd for specific application.

~Cheers from the PyCULA team.