benchmarks/pi.cpp

/*******************************************************
 * Copyright (c) 2014, ArrayFire
 * All rights reserved.
 *
 * This file is distributed under 3-clause BSD license.
 * The complete license agreement can be obtained at:
 * http://arrayfire.com/licenses/BSD-3-Clause
 ********************************************************/
 
/*
   monte-carlo estimation of PI
 
   algorithm:
   - generate random (x,y) samples uniformly
   - count what percent fell inside (top quarter) of unit circle
*/
 
#include <arrayfire.h>
#include <math.h>
#include <stdio.h>
#include <cstdlib>
using namespace af;
 
// generate millions of random samples
static int samples = 20e6;
 
/* Self-contained code to run host and device estimates of PI.  Note that
   each is generating its own random values, so the estimates of PI
   will differ. */
static double pi_device() {
    array x = randu(samples, f32), y = randu(samples, f32);
    return 4.0 * sum<float>(sqrt(x * x + y * y) < 1) / samples;
}
 
static double pi_host() {
    int count = 0;
    for (int i = 0; i < samples; ++i) {
        float x = float(rand()) / float(RAND_MAX);
        float y = float(rand()) / float(RAND_MAX);
        if (sqrt(x * x + y * y) < 1) count++;
    }
    return 4.0 * count / samples;
}
 
// void wrappers for timeit()
static void device_wrapper() { pi_device(); }
static void host_wrapper() { pi_host(); }
 
int main(int argc, char** argv) {
    try {
        int device = argc > 1 ? atoi(argv[1]) : 0;
        setDevice(device);
        info();
 
        printf("device:  %.5f seconds to estimate  pi = %.5f\n",
               timeit(device_wrapper), pi_device());
        printf("  host:  %.5f seconds to estimate  pi = %.5f\n",
               timeit(host_wrapper), pi_host());
    } catch (exception& e) {
        fprintf(stderr, "%s\n", e.what());
        throw;
    }
 
    return 0;
}