# Running Monte Carlo Simulations on Multiple GPUs

Today I'd like to introduce James Lebak. James is a developer who works on GPU support in the Parallel Computing Toolbox.

### Contents

#### Basic Option Pricing on the GPU

One common use of multi-GPU systems is to perform Monte Carlo simulations. The use of more GPUs can increase the number of samples that can be simulated, resulting in a more accurate simulation.

Let's begin by looking at the basic option pricing simulation from the Parallel Computing Toolbox example "Exotic Option Pricing on a GPU using a Monte-Carlo Method". In this example, we run many concurrent simulations of the evolution of the stock price. The mean and distribution of the simulation outputs gives us a sense of the ultimate value of the stock.

The function `simulateStockPrice` that describes the stock price is a discretization of a stochastic differential equation. The equation assumes that prices evolve according to a log-normal distribution related to the risk-free interest rate, the dividend yield (if any), and the volatility in the market.

```
type simulateStockPrice;
```

function finalStockPrice = simulateStockPrice(stockPrice, rate,... dividend, volatility,... Tfinal, dT) % Discrete simulation of stock price t = 0; while t < Tfinal t = t + dT; dr = (rate - dividend - volatility*volatility/2)*dT; perturbation = volatility*sqrt(dT)*randn(); stockPrice = stockPrice*exp(dr + perturbation); end finalStockPrice = stockPrice; end

Achieving good performance on a GPU requires executing many operations in parallel. The GPU has thousands of individual computational units, and you get best performance from them when you execute hundreds of thousands or millions of parallel operations. This is similar to the advice we often give to vectorize functions in MATLAB for better performance. The function `simulateStockPrice` produces a single, scalar output, so if we run it as it is, it won't achieve good performance.

We can execute the simulation many times on the GPU using ARRAYFUN. The `arrayfun` function on the GPU takes a function handle and a `gpuArray` input, and executes the function on each element of the input `gpuArray`. This is an alternative to vectorizing the function, which would require changing more code. It returns one output for each element in the input. This is exactly what we want to do in a Monte Carlo simulation: we want to produce many independent outputs.

The `runSimulationOnOneGPU` function uses `arrayfun` to perform many simulations in parallel and returns the mean price.

```
type runSimulationOnOneGPU.m;
```

function mfp = runSimulationOnOneGPU(Nsamples) % Run a single stock simulation on the GPU and return the % mean file price on the CPU. stockPrice = 100; % Stock price starts at $100. dividend = 0.01; % 1% annual dividend yield. riskFreeRate = 0.005; % 0.5 percent. timeToExpiry = 2; % Lifetime of the option in years. sampleRate = 1/250; % Assume 250 working days per year. volatility = 0.20; % 20% volatility. % Create the input data. Any scalar inputs are expanded to the size of the % array inputs. In this case, the starting stock prices is a vector whose % length is the number of simulations to perform on the GPU at once, and % all of the other inputs are scalars. startPrices = stockPrice*gpuArray.ones(Nsamples, 1); % Run all Nsamples simulations in parallel with one call to arrayfun. finalPrices = arrayfun( @simulateStockPrice, ... startPrices, riskFreeRate, dividend, volatility, ... timeToExpiry, sampleRate ); mfp = gather(mean(finalPrices)); end

We execute this function with a relatively small number of samples and display the computed mean.

```
nSamples = 1e6;
meanFinalPrice = runSimulationOnOneGPU(nSamples);
disp(['Calculated mean final price of ', num2str(meanFinalPrice)]);
```

Calculated mean final price of 98.96

#### Using Multiple GPUs on One Machine

Assume that we want to run this simulation on a desktop machine with two GPUs. This will allow us to have more samples in the same amount of time, which by the law of large numbers should give us a better approximation of the mean final stock price.

With Parallel Computing Toolbox, MATLAB can perform calculations on multiple GPUs. In order to do so, we need to open a MATLAB pool with one worker for each GPU. One MATLAB worker is needed to communicate with each GPU.

The workers in the pool will do all the calculations, so the client doesn't need to use a `gpuDevice`. Deselect the device on the client, so that all the memory on all the devices is completely available to the workers.

gpuDevice([]);

Determine how many GPUs we have on this machine and open a local MATLAB pool of that size. This gives us one worker for each GPU on this machine.

```
nGPUs = gpuDeviceCount();
matlabpool('local', nGPUs);
```

Starting matlabpool using the 'local' profile ... connected to 2 workers.

The `runSimulationOnManyGPUs` function contains a `parfor` loop that executes `nIter` times. Each iteration of the `parfor` loop executes `nSamples` independent Monte Carlo simulations on one GPU and returns the mean final price over all those simulations. The output of the overall simulation is the mean of the individual means, because each `parfor` iteration executes the same number of independent simulations. When the loop finishes we have executed `nIter*nSamples` simulations.

```
type runSimulationOnManyGPUs.m;
```

function [tout, mfp] = runSimulationOnManyGPUs(nSamples, nIter) % Executes a total of nSamples*nIter simulations, distributing the % simulations over the workers in a matlabpool in groups of nSamples each. tic; parfor ix = 1:nIter meanFinalPrice(ix) = runSimulationOnOneGPU(nSamples); end mfp = mean(meanFinalPrice); tout = toc;

To run `nSamples` iterations on each GPU, we call `runSimulationOnManyGPUs` with `nIter` equal to the number of GPUs.

[tout, meanFinalPrice] = runSimulationOnManyGPUs(nSamples, nGPUs); disp(['Performing simulations on the GPU took ', num2str(tout), ' s']); disp(['Calculated mean final price of ', num2str(meanFinalPrice)]);

Performing simulations on the GPU took 8.2226 s Calculated mean final price of 98.9855

It is important that the results are returned as regular MATLAB arrays rather than as `gpuArrays`. If the result is returned as a `gpuArray` then the client will share a GPU with one of the workers, which would unnecessarily use more memory on the device and transfer data over the PCI bus.

#### Multi-GPU Execution Details

The recommended setup for MATLAB when using multiple GPUs is, as we discussed above, to open one worker for each GPU. Let's dig in a little more to understand the details of how the workers interact with GPUs.

When workers share a single machine with multiple GPUs, MATLAB automatically assigns each worker in a parallel pool to use a different GPU by default. You can see this by using the spmd command to examine the index of the device used by each worker.

spmd gd = gpuDevice; idx = gd.Index; disp(['Using GPU ',num2str(idx)]); end

Lab 1: Using GPU 1 Lab 2: Using GPU 2

NVIDIA GPUs can be operated in one of four compute modes: default, exclusive thread, exclusive process, or prohibited. The compute mode is shown in the 'ComputeMode' field of the structure returned by gpuDevice. If a GPU is in prohibited mode, no worker will be assigned to use that GPU. If a GPU is in exclusive process or exclusive thread mode, only one MATLAB worker will attempt to access that GPU.

It is possible for multiple workers to share the same GPU, if the GPU is in 'Default' compute mode. When this is done, the GPU driver serializes accesses to the device. You should be aware that there will be a performance penalty because of the serialization. There will also be less memory available on the GPU for each worker, because two or more workers will be sharing the same GPU memory space. In general it will be hard to achieve good performance when multiple workers share the same GPU.

You can customize the way that MATLAB assigns workers to GPUs to suit your own needs by overriding the MATLAB function `selectGPU`. See the help for `selectGPU` for more details.

MATLAB initializes each worker to use a different random number stream on the GPU by default. Doing so ensures that the values obtained by each worker are different. In this example, we have for simplicity elected to use the default random number generation stream on the GPU, the well-known MRG32K3A stream. While MRG32K3A is a highly regarded algorithm with good support for parallelism, there are other streams that you can select that may provide better performance. The documentation page Using GPUArray describes the options available for controlling random number generation on the GPU, and lists the different streams that you can select.

#### Using Multiple GPUs in a Cluster

With all of these details in place, it's time to consider how to put multiple GPUs to use when hundreds of millions or even billions of simulations are desired. Such a situation might arise, for example, when we want different values of key parameters such as volatility or the dividend in each simulation.

Let's assume we want to run 800 million simulations, and that we have a cluster with 16 NVIDIA C2050 GPUs available. The problem is that a single GPU has a limited amount of memory. A single NVIDIA C2050 compute card has about 3 GB of memory and can compute on the order of $5\times10^7$ points in `runSimulationOnOneGPU` without running out of memory. We can call `runSimulationOnOneGPU` 16 times to complete all the desired simulations. On a single machine, doing so takes around 10 minutes.

We measured the time that it takes to run all $8\times 10^8$ simulations with up to 16 GPUs. As you can see, we got nearly-perfect linear scaling in this experiment.

How many GPUs do you use in your application, and how do they interact? Let me know in the comments.