My MATLAB GPU blog

Here are the results of my experiments with MATLAB's GPU capabilities -- these experiments all require the Parallel Computing Toolbox.

All timing results are on my HP Z420 PC which has a 4-core Intel Xeon E5-1620 CPU running at 3.6GHz, and an NVIDIA GeForce GTX 670 with 1344 cores running at 1GHz. Note that the double precision performance on the NVIDIA GeForce cards is much less than the single precision performance.

25/02/13:

My first experiment is a really simple Monte Carlo testcase: mc.m.

The single precision is about 13 times faster on the GPU, while the double precision is about 6 times faster. I think this comparison is with a CPU execution using 3 CPU threads.

The code experiments with the vector length, and also demonstrates the built-in function "gpuDevice" and the output from "whos".

My second experiment (kernel_test.m, kernel_test.cu) is based on this Mathworks webpage.

It demonstrates the ability of MATLAB to execute a CUDA kernel very simply, without the complexity of using the standard MEX interface. Further details about kernel parameters which can be set are given on this Mathworks webpage.

MATLAB has a built-in function called "paralleldemo_gpu_backslash" which assesses the performance of A\b on the GPU compared to the CPU.

On my system the single precision is up to 4 times faster on the GPU, but the double precision is only 45% faster.

More information on this testcase is available on this Mathworks webpage

This Mathworks webpage has this zip file with an application which benchmarks GPU performance for matrix multiplication, backslash (solving Ax=b) and FFTs.

26/02/13:

I generally find that Mathworks documentation is good once you find the thing you are looking for, but it can be frustrating initially when looking for information without a good idea of what's available and what the functions are called.

A good place to start looking is this Mathworks page: GPU computing

Another built-in demo function is called "paralleldemo_gpu_optionpricing" which does lots of Monte Carlo path calculations for option pricing.

This is a particularly good example of the use of "arrayfun" to create vectorised code -- I think I'll come back to this in the future to construct my own test codes.

Running it on my system, it immediately aborts with a message that "This example cannot run if kernel execution on the GPU can time-out." This is because my GPU also drives my display, and so there may be a timeout if a kernel runs for more than 5 seconds. (This is definitely a problem with Windows systems -- I'm not sure it's a problem with Linux systems.) Anyway, this code (my_paralleldemo_gpu_optionpricing.m) is the same code with the safety check commented out and the function renamed. The kernels ran for less than 2 seconds on my system, so this wasn't a problem.

More information on this testcase is available on this Mathworks webpage

More GPU examples are available on this Mathworks webpage -- click on "GPU Computing with MATLAB"

The first FFT example doesn't seem to be available as a file, so I've cut-and-pasted it here (fft_example.m).

One thing to notice in the code is the inheritance of array type; "timeVec" is explicitly declared to be a gpuArray (i.e.~an array that lives on the GPU) and so therefore "signal", "transformedSignal" and "powerSpectrum" are all gpuArrays. I've added a "whos" at the end of the code to check this.

However, "randn(size(timeVec))" performs the random number generation on the CPU -- to do it on the GPU it would need to be "gpuArray.randn(size(timeVec))"

To construct simple matrices on the GPU, such as all zeros or all ones, there are special "static methods" which are documented here. Note in particular gpuArray.colon which is the GPU equivalent of A:D:B. There is also information there on random number generation on the GPU using gpuArray.rand.

There is currently no support for sparse matrices on GPUs -- I don't know whether this might come in the future.

27/02/13:

Here's my new Monte Carlo testcase: mc2.m.

This uses arrayfun to implement a Monte Carlo option pricing calculation with 2D paths. The speedup factor for this is approximately 30, but this higher speedup factor is partly because the CPU implementation (which uses exactly the same function code) seems to use only 1 thread.

This suggests that the vector operations in my first MC code (mc.m) are multi-threaded, but not the arrayfun in my new code. That strikes me as odd, but maybe I'm missing something. Perhaps it's connected with the issues involved in parallelising random number generation and getting the same results regardless of the number of threads used.

Another point to note is that the arguments in the function called by arrayfun must all be scalars. In this 2D example it means the starting values for each dimension have to be passed in separately. In a high dimensional case (e.g. dimension 50?) this could be a nuisance. In the CPU version, I think this can be handled by using cell arrays, with each cell array element holding a vector, but this is not an option for the GPU version and I don't see an obvious alternative.

Here's an updated version of the same testcase: mc2c.m.

This compares the performance using arrayfun and the vectorisation used before in mc.m. The results indicate that for CPU execution it is much better to use the vector form (which is multithreaded) than the arrayfun form (which is not). For GPU execution, it is slightly better to use the arrayfun form, probably because it leads to a single CUDA kernel whereas the vector form may lead to a separate kernel for each pass through the timestep iteration loop -- this is why the execution is a bit faster in a second vector version which uses partial loop unrolling.

Here's a really simple 2D finite difference application: fd.m.

The speedup factor is just over 10 in single precision, and almost 8 in double precision, relative to the CPU execution which appears to use 3 threads. This is very good performance with really simple programming -- the core piece of code is the same for all 4 versions (CPU/GPU, single/double precision).

28/02/13:

In a number of my experiments above I'm checking on performance using both single and double precision arithmetic.

MATLAB users should be aware that if S and D are single and double precision variables, respectively, then S+D and S*D are single precision variables. This is documented on this Mathworks webpage. Although it doesn't say, I suspect that D is converted to single precision and then the operation is carried out, rather than S being converted to double precision, the operation being done in double precision, and then the result being converted to single precision.

This behaviour is different from C/C++ where the output would be double precision, with the value of S being "promoted" to double precision before the operation is carried out (see, for example, this webpage).

Using MATLAB on a CPU, I would almost always use double precision, which is the default. There is almost no performance difference between single and double precision, and so the only benefit of single precision is that it requires less memory for big arrays, and that's rarely a problem for me in my applications.

However, on the GPU there can be a significant performance difference between single precision and double precision arithmetic, and so that's why I try both in my tests.

A comment about timing: it's best to run an application twice to get good timing results. The first time, there can be some extra overhead (checking which GPU is available, or maybe generating the CUDA code for execution on the GPU?) which disappears in subsequent runs.

04/03/13:

My new Apple MacBook Pro has a 384-core NVIDIA GT 650M GPU. Compared to my desktop GTX 670, it has 2/7 of the cores, 1/2 of the graphics memory (1GB versus 2GB) and about 2/5 of the memory bandwidth (up to 80 GB/s versus 190 GB/s)
Installing MATLAB was easy as usual, but trying to run GPU code produced an error because it seems the Parallel Computing toolbox does not come with a built-in copy of the CUDA toolkit. Installing CUDA (again a simple task) fixed this.
One oddity: after installing CUDA I was able to use command "gpuDevice" to check the GPU I had, but running application code produced nonsense. Being used to the peculiarities of software installation, I rebooted the system, and now things seem to work correctly. This advice might have been included in the installation documentation, but of course I didn't read it.
Initial experiments give a disappointing performance, roughly a factor 10 slower than my GTX 670. I'll need to investigate this and seek help from others.

06/03/13:

Jos Martin, who heads the development of the Parallel Computing toolbox at Mathworks, came to Oxford today and very kindly gave us two seminars, one on Mathworks approach to software engineering, and a second on the development of the GPU capabilities in MATLAB. He's going to send us a copy of his slides and then I'll put a link here.

Following discussions with Jos, I've cleaned up a number of my codes above so that they pass most of the MATLAB lint tests, and I correctly initialise the GPU random number generator in mc.m. Also, I've put in a "wait(gpu)" instruction to make sure the GPU operations have completed before "toc" computes the elapsed time.

One other thing I discovered while showing Jos the demo codes on my MacBook is that the performance is now much better than on Monday. I hate it when that happens, although I hate it more when the performance gets worse. One of the things I did on Monday was to update the CUDA driver, but when I tested it on Monday the performance was no better. However, I may not have rebooted the system, and it's possible (maybe even likely) that was required. I suspect this is why it ran faster today.

The GPU performance now is roughly a factor 3.5 slower than the GTX 670, which is about what I would expect.

A tip from Jos: to find all of the built-in GPU demos, type "paralleldemo_gpu" in the Command Window and then Tab for auto-completion -- there are 7 of them, including the option pricing example discussed above:
- paralleldemo_gpu_arrayfun
- paralleldemo_gpu_gpu_backslash
- paralleldemo_gpu_devices
- paralleldemo_gpu_fft
- paralleldemo_gpu_fft2
- paralleldemo_gpu_mandelbrot
- paralleldemo_gpu_optionpricing

07/03/13:

Sobol quasi-random numbers/points are very useful for Monte Carlo simulation -- in the best cases, with N points the error is O(1/N), whereas with regular Monte Carlo it would be O(1/sqrt(N)). I was very pleased when Mathworks added sobolset, set and scramble to the Statistics toolbox, but unfortunately there are no GPU versions of these yet.