From 69f2911e04ffb1b19eef1fafb8c040af271f656e Mon Sep 17 00:00:00 2001
From: Tor Aamodt <aamodt@ece.ubc.ca>
Date: Thu, 15 Jul 2010 18:09:46 -0800
Subject: creating branch for adding support for CUDA 3.x and Fermi

[git-p4: depot-paths = "//depot/gpgpu_sim_research/fermi/distribution/": change = 6829]
---
 benchmarks/CUDA/LIB/libor.cu | 387 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 387 insertions(+)
 create mode 100644 benchmarks/CUDA/LIB/libor.cu

(limited to 'benchmarks/CUDA/LIB/libor.cu')

diff --git a/benchmarks/CUDA/LIB/libor.cu b/benchmarks/CUDA/LIB/libor.cu
new file mode 100644
index 0000000..a245862
--- /dev/null
+++ b/benchmarks/CUDA/LIB/libor.cu
@@ -0,0 +1,387 @@
+
+/* Program to compute swaption portfolio using NVIDIA CUDA */
+
+#include <stdio.h>
+#include <cutil.h>
+
+// parameters for nVidia device execution
+
+#define BLOCK_SIZE 64
+#define GRID_SIZE 64
+
+// parameters for LIBOR calculation
+
+#define NN 80
+#define NMAT 40
+#define L2_SIZE 3280 //NN*(NMAT+1)
+#define NOPT 15
+#define NPATH 4096
+
+// constant data for swaption portfolio: stored in device memory,
+// initialised by host and read by device threads
+
+__constant__ int    N, Nmat, Nopt, maturities[NOPT]; 
+__constant__ float  delta, swaprates[NOPT], lambda[NN];
+
+
+/* Monte Carlo LIBOR path calculation */
+
+__device__ void path_calc(float *L, float *z)
+{
+  int   i, n;
+  float sqez, lam, con1, v, vrat;
+
+  for(n=0; n<Nmat; n++) {
+    sqez = sqrtf(delta)*z[n];
+    v = 0.0;
+
+    for (i=n+1; i<N; i++) {
+      lam  = lambda[i-n-1];
+      con1 = delta*lam;
+      v   += __fdividef(con1*L[i],1.0+delta*L[i]);
+      vrat = __expf(con1*v + lam*(sqez-0.5*con1));
+      L[i] = L[i]*vrat;
+    }
+  }
+}
+
+
+/* forward path calculation storing data
+   for subsequent reverse path calculation */
+
+__device__ void path_calc_b1(float *L, float *z, float *L2)
+{
+  int   i, n;
+  float sqez, lam, con1, v, vrat;
+
+  for (i=0; i<N; i++) L2[i] = L[i];
+   
+  for(n=0; n<Nmat; n++) {
+    sqez = sqrt(delta)*z[n];
+    v = 0.0;
+
+    for (i=n+1; i<N; i++) {
+      lam  = lambda[i-n-1];
+      con1 = delta*lam;
+      v   += __fdividef(con1*L[i],1.0+delta*L[i]);
+      vrat = __expf(con1*v + lam*(sqez-0.5*con1));
+      L[i] = L[i]*vrat;
+
+      // store these values for reverse path //
+      L2[i+(n+1)*N] = L[i];
+    }
+  }
+}
+
+
+/* reverse path calculation of deltas using stored data */
+
+__device__ void path_calc_b2(float *L_b, float *z, float *L2)
+{
+  int   i, n;
+  float faci, v1;
+
+  for (n=Nmat-1; n>=0; n--) {
+    v1 = 0.0;
+    for (i=N-1; i>n; i--) {
+      v1    += lambda[i-n-1]*L2[i+(n+1)*N]*L_b[i];
+      faci   = __fdividef(delta,1.0+delta*L2[i+n*N]);
+      L_b[i] = L_b[i]*__fdividef(L2[i+(n+1)*N],L2[i+n*N])
+              + v1*lambda[i-n-1]*faci*faci;
+ 
+    }
+  }
+}
+
+/* calculate the portfolio value v, and its sensitivity to L */
+/* hand-coded reverse mode sensitivity */
+
+__device__ float portfolio_b(float *L, float *L_b) 
+{
+  int   m, n;
+  float b, s, swapval,v;
+  float B[NMAT], S[NMAT], B_b[NMAT], S_b[NMAT];
+
+  b = 1.0;
+  s = 0.0;
+  for (m=0; m<N-Nmat; m++) {
+    n    = m + Nmat;
+    b    = __fdividef(b,1.0+delta*L[n]);
+    s    = s + delta*b;
+    B[m] = b;
+    S[m] = s;
+  }
+
+  v = 0.0;
+
+  for (m=0; m<N-Nmat; m++) {
+    B_b[m] = 0;
+    S_b[m] = 0;
+  }
+
+  for (n=0; n<Nopt; n++){
+    m = maturities[n] - 1;
+    swapval = B[m] + swaprates[n]*S[m] - 1.0;
+    if (swapval<0) {
+      v     += -100*swapval;
+      S_b[m] += -100*swaprates[n];
+      B_b[m] += -100;
+    }
+  }
+
+  for (m=N-Nmat-1; m>=0; m--) {
+    n = m + Nmat;
+    B_b[m] += delta*S_b[m];
+    L_b[n]  = -B_b[m]*B[m]*__fdividef(delta,1.0+delta*L[n]);
+    if (m>0) {
+      S_b[m-1] += S_b[m];
+      B_b[m-1] += __fdividef(B_b[m],1.+delta*L[n]);
+    }
+  }
+
+  // apply discount //
+
+  b = 1.0;
+  for (n=0; n<Nmat; n++) b = b/(1.0+delta*L[n]);
+
+  v = b*v;
+
+  for (n=0; n<Nmat; n++){
+    L_b[n] = -v*delta/(1.0+delta*L[n]);
+  }
+
+  for (n=Nmat; n<N; n++){
+    L_b[n] = b*L_b[n];
+  }
+
+  return v;
+}
+
+
+/* calculate the portfolio value v */
+
+__device__ float portfolio(float *L)
+{
+  int   n, m, i;
+  float v, b, s, swapval, B[40], S[40];
+	
+  b = 1.0;
+  s = 0.0;
+
+  for(n=Nmat; n<N; n++) {
+    b = b/(1.0+delta*L[n]);
+    s = s + delta*b;
+    B[n-Nmat] = b;
+    S[n-Nmat] = s;
+  }
+
+  v = 0.0;
+
+  for(i=0; i<Nopt; i++){
+    m = maturities[i] -1;
+    swapval = B[m] + swaprates[i]*S[m] - 1.0;
+    if(swapval<0)
+      v += -100.0*swapval;
+  }
+
+  // apply discount //
+
+  b = 1.0;
+  for (n=0; n<Nmat; n++) b = b/(1.0+delta*L[n]);
+
+  v = b*v;
+
+  return v;
+}
+
+
+__global__ void Pathcalc_Portfolio_KernelGPU(float *d_v, float *d_Lb)
+{
+  const int     tid = blockDim.x * blockIdx.x + threadIdx.x;
+  const int threadN = blockDim.x * gridDim.x;
+
+  int   i,path;
+  float L[NN], L2[L2_SIZE], z[NN];
+  float *L_b = L;
+  
+  /* Monte Carlo LIBOR path calculation*/
+
+  for(path = tid; path < NPATH; path += threadN){
+    // initialise the data for current thread
+    for (i=0; i<N; i++) {
+      // for real application, z should be randomly generated
+      z[i] = 0.3;
+      L[i] = 0.05;
+    }
+    path_calc_b1(L, z, L2);
+    d_v[path] = portfolio_b(L,L_b);
+    path_calc_b2(L_b, z, L2);
+    d_Lb[path] = L_b[NN-1];
+  }
+}
+
+
+__global__ void Pathcalc_Portfolio_KernelGPU2(float *d_v)
+{
+  const int     tid = blockDim.x * blockIdx.x + threadIdx.x;
+  const int threadN = blockDim.x * gridDim.x;
+
+  int   i, path;
+  float L[NN], z[NN];
+  
+  /* Monte Carlo LIBOR path calculation*/
+
+  for(path = tid; path < NPATH; path += threadN){
+    // initialise the data for current thread
+    for (i=0; i<N; i++) {
+      // for real application, z should be randomly generated
+      z[i] = 0.3;
+      L[i] = 0.05;
+    }	   
+    path_calc(L, z);
+    d_v[path] = portfolio(L);
+  }
+}
+
+
+////////////////////////////////////////////////////////////////////////
+// Main program
+////////////////////////////////////////////////////////////////////////
+
+int main(int argc, char **argv){
+    
+  // 'h_' prefix - CPU (host) memory space
+
+  float  *h_v, *h_Lb, h_lambda[NN], h_delta=0.25;
+  int     h_N=NN, h_Nmat=NMAT, h_Nopt=NOPT, i;
+  int     h_maturities[] = {4,4,4,8,8,8,20,20,20,28,28,28,40,40,40};
+  float   h_swaprates[]  = {.045,.05,.055,.045,.05,.055,.045,.05,
+                            .055,.045,.05,.055,.045,.05,.055 };
+  double  v, Lb; 
+
+  unsigned int hTimer;
+  double  gpuTime;
+
+  // 'd_' prefix - GPU (device) memory space
+
+  float  *d_v,*d_Lb;
+
+  // initialise card and timer
+  int deviceCount;                                                         
+  CUDA_SAFE_CALL_NO_SYNC(cudaGetDeviceCount(&deviceCount));                
+  if (deviceCount == 0) {                                                  
+      fprintf(stderr, "There is no device.\n");                            
+      exit(EXIT_FAILURE);                                                  
+  }                                                                        
+  int dev;                                                                 
+  for (dev = 0; dev < deviceCount; ++dev) {                                
+      cudaDeviceProp deviceProp;                                           
+      CUDA_SAFE_CALL_NO_SYNC(cudaGetDeviceProperties(&deviceProp, dev));   
+      if (deviceProp.major >= 1)                                           
+          break;                                                           
+  }                                                                        
+  if (dev == deviceCount) {                                                
+      fprintf(stderr, "There is no device supporting CUDA.\n");            
+      exit(EXIT_FAILURE);                                                  
+  }                                                                        
+  else                                                                     
+      CUDA_SAFE_CALL(cudaSetDevice(dev));  
+
+  CUT_SAFE_CALL( cutCreateTimer(&hTimer) );
+
+  for (i=0; i<NN; i++) h_lambda[i] = 0.2;
+
+  // Copy all constants into constant memory
+
+  cudaMemcpyToSymbol(N, &h_N, sizeof(h_N));
+  cudaMemcpyToSymbol(Nmat, &h_Nmat, sizeof(h_Nmat));
+  cudaMemcpyToSymbol(Nopt, &h_Nopt, sizeof(h_Nopt));
+  cudaMemcpyToSymbol(delta, &h_delta, sizeof(h_delta));
+  cudaMemcpyToSymbol(maturities, &h_maturities, sizeof(h_maturities));
+  cudaMemcpyToSymbol(swaprates, &h_swaprates, sizeof(h_swaprates));
+  cudaMemcpyToSymbol(lambda, &h_lambda, sizeof(h_lambda));        
+
+  // Allocate memory on host and device
+
+  h_v      = (float *)malloc(sizeof(float)*NPATH);
+  CUDA_SAFE_CALL( cudaMalloc((void **)&d_v, sizeof(float)*NPATH) );
+  h_Lb     = (float *)malloc(sizeof(float)*NPATH);
+  CUDA_SAFE_CALL( cudaMalloc((void **)&d_Lb, sizeof(float)*NPATH) );
+
+  // Execute GPU kernel -- no Greeks
+
+  CUDA_SAFE_CALL( cudaThreadSynchronize() );
+  CUT_SAFE_CALL( cutResetTimer(hTimer) );
+  CUT_SAFE_CALL( cutStartTimer(hTimer) );
+    
+  // Set up the execution configuration
+
+  dim3 dimBlock(BLOCK_SIZE);
+  dim3 dimGrid(GRID_SIZE);
+
+  // Launch the device computation threads
+
+  Pathcalc_Portfolio_KernelGPU2<<<dimGrid, dimBlock>>>(d_v);
+  CUT_CHECK_ERROR("Pathcalc_Portfolio_kernelGPU2() execution failed\n");
+  CUDA_SAFE_CALL( cudaThreadSynchronize() );
+
+  // Read back GPU results and compute average
+
+  CUDA_SAFE_CALL( cudaMemcpy(h_v, d_v, sizeof(float)*NPATH,
+                  cudaMemcpyDeviceToHost) );
+  CUT_SAFE_CALL( cutStopTimer(hTimer) );
+  gpuTime = cutGetTimerValue(hTimer);
+
+  v = 0.0;
+  for (i=0; i<NPATH; i++) v += h_v[i];
+  v = v / NPATH;
+
+  printf("v  = %15.8f\n", v);
+  printf("Time(No Greeks) : %f msec\n", gpuTime);
+
+  // Execute GPU kernel -- Greeks
+
+  CUDA_SAFE_CALL( cudaThreadSynchronize() );
+  CUT_SAFE_CALL( cutResetTimer(hTimer) );
+  CUT_SAFE_CALL( cutStartTimer(hTimer) ); 
+
+  // Launch the device computation threads
+
+  Pathcalc_Portfolio_KernelGPU<<<dimGrid, dimBlock>>>(d_v,d_Lb);
+  CUT_CHECK_ERROR("Pathcalc_Portfolio_kernelGPU() execution failed\n");
+  CUDA_SAFE_CALL( cudaThreadSynchronize() );
+
+  // Read back GPU results and compute average
+
+  CUDA_SAFE_CALL( cudaMemcpy(h_v, d_v, sizeof(float)*NPATH,
+                  cudaMemcpyDeviceToHost) );
+  CUDA_SAFE_CALL( cudaMemcpy(h_Lb, d_Lb, sizeof(float)*NPATH,
+                  cudaMemcpyDeviceToHost) );
+  CUT_SAFE_CALL( cutStopTimer(hTimer) );
+  gpuTime = cutGetTimerValue(hTimer);
+    
+  v = 0.0;
+  for (i=0; i<NPATH; i++) v += h_v[i];
+  v = v / NPATH;
+
+  Lb = 0.0;
+  for (i=0; i<NPATH; i++) Lb += h_Lb[i];
+  Lb = Lb / NPATH;
+
+  printf("v  = %15.8f\n", v);
+  printf("Lb = %15.8f\n", Lb);
+  printf("Time (Greeks)   : %f msec\n", gpuTime);
+
+  // Release GPU memory
+
+  CUDA_SAFE_CALL( cudaFree(d_v));
+  CUDA_SAFE_CALL( cudaFree(d_Lb));
+       
+  // Release CPU memory
+
+  free(h_v);
+  free(h_Lb);
+
+  CUT_SAFE_CALL( cutDeleteTimer(hTimer) );
+  CUT_EXIT(argc, argv);
+}
-- 
cgit v1.3