Black scholes mgl

% This is a multi-GPU version of black_scholes_example.m
function black_scholes_mgl(varargin)
  disp '** Black-Scholes Multi-GPU Example **'
 
  %%-- setup --%%
  % pricing parameters
  risk = 0.02; % risk-free interest
  volatility = 0.30;
  ngpu = getfield(ginfo, 'gpu_count');
  % determine precision
  use_double = determine_precision(varargin);
  if use_double   [cpu_cls gpu_cls] = deal('double', @gdouble);
  else            [cpu_cls gpu_cls] = deal('single', @gsingle); end
 
  %%-- iterate over progressively larger option pools --%%
  fprintf(' %8s     %9s', 'options', 'cpu rate')
  for i = 1:ngpu, fprintf('            %d gpu', i), end
  fprintf('\n')
  fprintf(' --------     ---------')
  for i = 1:ngpu, fprintf('    ---------------'), end
  fprintf('\n')
 
  current_device = gselect;
  flops = getfield(ginfo('all'), 'flops');
 
  for options = 100e3*(1:20)
    % generate random input at this size
    stock_price   = randu(5,    30, [options 1], cpu_cls); % ~ U[5,30]
    strike_price  = randu(1,   100, [options 1], cpu_cls); % ~ U[1,100]
    maturity_time = randu(0.25, 10, [options 1], cpu_cls); % ~ U[0.25,10]
 
    % evaluate on CPU
    time = timeit(@blackscholes, stock_price, strike_price, risk, volatility, maturity_time);
    rate_cpu = options / time;
    fprintf(' %8d     %-9.3e', options, rate_cpu);
 
    % try spreading across more and more GPUs
    for n = 1:ngpu
      % transfer same input out to each GPU
      ratio = flops(1:n);
      splits = [0, ceil(cumsum(options * (ratio / sum(ratio))))];
      for i = 1:n
        range = [(splits(i)+1):splits(i+1)];
        stock{i}    = gpu_cls(stock_price(range, :));
        strike{i}   = gpu_cls(strike_price(range, :));
        maturity{i} = gpu_cls(maturity_time(range, :));
      end
 
      % execute
      time = timeit(@blackscholes_multi, n, stock, strike, risk, volatility, maturity,'all');
      rate_gpu = options / time;
      speedup = rate_gpu / rate_cpu;
 
      fprintf('    %-9.3e (%2.0fx)', rate_gpu, speedup);
    end
    fprintf('\n')
  end
 
  gselect(current_device);
end
 
 
function use_double = determine_precision(args)
  is_double = getfield(ginfo, 'compute') >= 1.3;
  if numel(args)
    switch args{1}
     case 'double',  use_double = true;
     case 'single',  use_double = false;
     otherwise       error('Input arguments are incorrect.')
    end
  else
    use_double = is_double;
  end
  if use_double
    disp '(using double-precision)'
  else
    disp '(using single-precision)'
  end
end
 
 
% uniform random distribution ~ U[min,max]
function X = randu(min,max, varargin)
  X = min + (max-min) * rand(varargin{:});
end
 
 
function blackscholes_multi(ngpu, stock, strike, risk, volatility, maturity)
  for i = 1:ngpu
    gselect(i)
    [c p] = blackscholes(stock{i}, strike{i}, risk, volatility, maturity{i});
    geval(c,p)
  end
end
 
 
% calculate call/put option prices based on Black-Scholes European model
% from Black, Scholes, "Pricing of Options and Corporate Liabilities" (1973)
function [Vcall, Vput] = blackscholes(S,X,R,V,T)
  % S = Underlying stock price
  % X = Strike Price
  % R = Risk free rate of interest
  % V = Volatility
  % T = Time to maturity
 
  d1 = (log(S ./ X) + (R + (V.*V)*0.5) .* T) ./ (V.*sqrt(T));
  d2 = d1 - (V.*sqrt(T));
 
  % cumulative normal distribution (probability that a variable will be less than x std)
  eps = 1e-16;
  cnd_d1 = 0.5 + sign(d1+eps) .* erf(sign(d1).*d1/sqrt(2))/2;
  cnd_d2 = 0.5 + sign(d2+eps) .* erf(sign(d2).*d2/sqrt(2))/2;
 
  Vcall = S .* cnd_d1  - X .* exp(-R.*T) .* cnd_d2;
  Vput  = X .* exp(-R.*T) .* (1-cnd_d2) - (S .* (1-cnd_d1));
end