{"id":9,"date":"2020-06-11T21:32:03","date_gmt":"2020-06-11T20:32:03","guid":{"rendered":"http:\/\/retrofocus.com\/?p=9"},"modified":"2021-01-09T22:29:47","modified_gmt":"2021-01-09T21:29:47","slug":"estimating-pi-using-the-monte-carlo-method-with-nvidia-gpu","status":"publish","type":"post","link":"http:\/\/retrofocus.com\/?p=9","title":{"rendered":"Estimating Pi using the Monte Carlo Method with Nvidia GPU"},"content":{"rendered":"\n

What does it feel to command 1280 parallel cores? Having earlier used GPU for machine learning with Keras, I wanted to have a look at directly programming the GPU with C\/C++. This is thus my first trial with GPU programming based on Nvidia CUDA framework. I’m using Nvidia 1060 GPU but I guess this would work on others as well. <\/p>\n\n\n\n

Here is my code: <\/p>\n\n\n\n

#include<stdio.h>\n#include <cuda.h>\n#include <curand_kernel.h>\n__global__ void mc_pi(unsigned long long *countti,long iter,unsigned long long *palluu){\n        int ii;\n        curandState_t rng;\n        int tid = threadIdx.x + blockIdx.x * blockDim.x;\n        \/\/int blockid=blockIdx.x*blockDim.x;\n        countti[tid]=0;\n        \/\/countti[threadIdx.x]=0;\n        curand_init(clock64(), tid, 0, &rng);\n        for(ii=0;ii<iter;ii++){\n        float x = curand_uniform(&rng); \n        float y = curand_uniform(&rng); \n        if ((x*x+y*y) <=1 )\n                {countti[tid] += 1 ;    \n                }\n        }\n        __syncthreads();        \n        if (!(tid&1)) {countti[tid]=countti[tid]+countti[tid+1];}\n        __syncthreads();\n        if (!(tid&3)) {countti[tid]=countti[tid]+countti[tid+2];}\n        __syncthreads();\n        if (!(tid&7)){countti[tid]=countti[tid]+countti[tid+4];}\n        __syncthreads();\n        if (!(tid&15)) {countti[tid]=countti[tid]+countti[tid+8];}\n        __syncthreads();\n        if (!(tid&31)) {countti[tid]=countti[tid]+countti[tid+16];}\n        __syncthreads();\n        if (!(tid&63)) {countti[tid]=countti[tid]+countti[tid+32];}\n        __syncthreads();\n        if (!(tid&127)) {countti[tid]=countti[tid]+countti[tid+64];}\n        __syncthreads();\n        if (threadIdx.x ==0) {countti[tid]=countti[tid]+countti[tid+128];\n        palluu[blockIdx.x]=countti[tid];}       \n        \/\/palluu[blockIdx.x]=blockIdx.x;\n        __syncthreads();\n        \n}\n\n__global__ void reduce_blocks(unsigned long long *palluu){\nint tid = threadIdx.x;\/\/ + blockIdx.x * blockDim.x;\n\nint i=0;\nint oso= 0;\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif ((oso&1) ==0 && oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+1];}}\n__syncthreads();\n\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif (oso%4 ==0 && oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+2];}}\n        __syncthreads();\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif (oso%8 ==0 && oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+4];}}\n        __syncthreads();\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif (oso%16 ==0 && oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+8];}}\n        __syncthreads();\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif (oso%32 ==0 && oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+16];}}\n        __syncthreads();\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif (oso%64 ==0 &&oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+32];}}\n        __syncthreads();\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif (oso%128 ==0 && oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+64];}}\n        __syncthreads();\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif (oso%256 ==0&&oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+128];}}\n        __syncthreads();\n\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif (oso%512 ==0&&oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+256];}}\n __syncthreads();\n\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif (oso%1024 ==0&&oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+512];}}\n __syncthreads();\n\n\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif (oso%2048 ==0&&oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+1024];}}\n __syncthreads();\n\nfor(i=0;i<6;i++) {\noso=tid+(i*blockDim.x);\nif (oso%4096 ==0&&oso<3906) {palluu[oso]=palluu[oso]+palluu[oso+2048];}}\n __syncthreads();\n\n\n}\n\n\n\nint main() {\n        int N=999936;\n        int blocks=(N+255)\/256;\n        int threadsperblock=256;\n        int iterations=1000000;\n        int  i=0;\n        long total=0;\n        long long totalos=0;\n  unsigned long long *countti;\n  unsigned long long *resu;\n  unsigned long long *hostresu,hostresu2;\n        printf(\"blocks:%d, threadsperblock:%d\\n\",blocks,threadsperblock);\n  hostresu=(unsigned long long*)malloc(blocks*sizeof(long long));\n  cudaMalloc(&countti, N*sizeof(unsigned long long));\n  cudaMalloc(&resu,blocks * sizeof(unsigned long long));\n   mc_pi<<<blocks,threadsperblock>>>(countti,iterations,resu); \n   \/\/reduce_blocks<<<1,651>>>(resu);\n  reduce_blocks<<<1,651>>>(resu);\n\n  cudaMemcpy((long long*)&hostresu2,resu,sizeof(long long),cudaMemcpyDeviceToHost);\n  printf(\"Estimating Pi using the Monte Carlo Method in Nvidia 1060 parallel computation\\n: %.15f, based on %llu iterations\\n\",(double)hostresu2\/threadsperblock\/blocks\/iterations*4,threadsperblock*blocks*iterations);\n\nreturn 0;\n}\n<\/code><\/pre>\n\n\n\n
I’m using two kernel functions: mc_pi is used for actual pi monte carlo computations as well as for reducing the results on cuda block level. Second kernel reduce_blocks is used for reduction over 3906 blocks. The block size and number of blocks were quickly picked to come close to one million threads in total. <\/p>\n\n\n\n
Compilation, when you have Nvidia cuda framework installed:<\/p>\n\n\n\n
nvcc -o pi4cu pi4.cu <\/code><\/pre>\n\n\n\n
Performance: <\/p>\n\n\n\n
time .\/pi4cu \nblocks:3906, threadsperblock:256\nEstimating Pi using the Monte Carlo Method in Nvidia 1060 parallel computation\n: 3.141592344990079, based on 3503587328 iterations\n\nreal    0m27.299s\nuser    0m18.580s\nsys     0m8.704s\n<\/code><\/pre>\n\n\n\n
So 27 seconds it takes for 3.5 billion iterations maxing GPU to 100%. When compared again my single threaded python version of pi monte carlo function below, the performance improvement is 87x for GPU version. I’m almost sure this is still far from optimal. I read e.g. that GPU isn’t fast in the modulo operations in the reducers, thus in the mc_pi I already changed the modulos to bitwise operators but didn’t yet do it for the block level reducer. <\/p>\n\n\n\n
For reducers generally there would be different approaches to choose between. You can read more about those here: https:\/\/developer.download.nvidia.com\/assets\/cuda\/files\/reduction.pdf<\/a><\/p>\n\n\n\n
<\/p>\n\n\n\n
import random\nimport math\n\nsis=0\nlukum=40000000\n\nfor i in range(lukum):\n        a=random.random()\n        b=random.random()\n        if (a**2+b**2) <= 1:\n                sis+=1\nprint(\"4*sis\/lukum:\", 4*sis\/lukum)\n<\/code><\/pre>\n\n\n\n
Accuracy? <\/p>\n\n\n\n
My estimate is 3.141592344990079. Let's ask google:\n\npi-3.141592344990079<\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
What next with cuda – I would like to try some financial monte carlo simulations or perhaps deoloying graph algorithms on GPU. <\/p>\n","protected":false},"excerpt":{"rendered":"
What does it feel to command 1280 parallel cores? Having earlier used GPU for machine learning with Keras, I wanted to have a look at…<\/p>\n
Continue readingEstimating Pi using the Monte Carlo Method with Nvidia GPU<\/span><\/a><\/div>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":[],"categories":[3,4],"tags":[12],"_links":{"self":[{"href":"http:\/\/retrofocus.com\/index.php?rest_route=\/wp\/v2\/posts\/9"}],"collection":[{"href":"http:\/\/retrofocus.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"http:\/\/retrofocus.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"http:\/\/retrofocus.com\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"http:\/\/retrofocus.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=9"}],"version-history":[{"count":8,"href":"http:\/\/retrofocus.com\/index.php?rest_route=\/wp\/v2\/posts\/9\/revisions"}],"predecessor-version":[{"id":40,"href":"http:\/\/retrofocus.com\/index.php?rest_route=\/wp\/v2\/posts\/9\/revisions\/40"}],"wp:attachment":[{"href":"http:\/\/retrofocus.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=9"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"http:\/\/retrofocus.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=9"},{"taxonomy":"post_tag","embeddable":true,"href":"http:\/\/retrofocus.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=9"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}