Calculation of Pi

After my initial encountering with OpenCL, I wanted to explore more. So I wrote a small program to calculate Pi .
This uses a very simple algorithm to calculate Pi

  1. Inscribe a circle in a square
  2. Randomly generate points in the square
  3. Determine the number of points in the square that are also in the circle
  4. Let r be the number of points in the circle divided by the number of points in the square
  5. PI ~ 4 r
  6. Note that the more points generated, the better the approximation

This gave the value of Pi as 3.141172 . Pretty decent, I suppose 🙂
 

The code

 

import pyopencl as cl 
import numpy
import random

class circle(object):
    def __init__(self, r=100):
        self.r = r
    def inside_the_circle(self, x, y):
        if x**2+y**2 <= r**2: 
           return 1
        else:
           return 0

class CL:
    def __init__(self, size=1000):
        self.size = size
        self.ctx = cl.create_some_context()
        self.queue = cl.CommandQueue(self.ctx)

    def load_program(self):
        f = """
        __kernel void picalc(__global float* a, __global float* b, __global float* c)
        {
         unsigned int i = get_global_id(0);

        if (a[i]*a[i]+ b[i]*b[i] < 100*100)
           {
            c[i] = 1;
           }
        else
           {
            c[i]=0;
           }
        }
        """
        self.program = cl.Program(self.ctx, f).build()

    def popcorn(self):
        mf = cl.mem_flags
        x = [random.uniform(-100,100) for i in range(self.size)]

        y = [random.uniform(-100,100) for i in range(self.size)]

        self.a = numpy.array(x, dtype=numpy.float32)
        self.b = numpy.array(y, dtype=numpy.float32)

        self.a_buf = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.a)
        self.b_buf = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.b)
        self.dest_buf = cl.Buffer(self.ctx, mf.WRITE_ONLY, self.b.nbytes)

    def execute(self):
        self.program.picalc(self.queue, self.a.shape, None, self.a_buf, 
                            self.b_buf, self.dest_buf)
        self.c = numpy.empty_like(self.a)
        cl.enqueue_read_buffer(self.queue, self.dest_buf, self.c).wait()
        
    def calculate_pi(self):
        number_in_circle = 0
        for i in self.c:
            number_in_circle = number_in_circle + i
        pi = number_in_circle*4 / self.size
        print 'pi = ', pi



if __name__ == '__main__':
    ex = CL(1000000)
    ex.load_program()
    ex.popcorn()
    ex.execute()
    ex.calculate_pi()

 

Advertisements

One thought on “Calculation of Pi

  1. Thanks your code has helped me a lot to start learning OpenCL. I will compare it to my current Parallel Python code for performance. Thank you.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s