ArrayFire_demo: New package

git-svn-id: svn://ultimatepp.org/upp/trunk@9359 f0d560ea-af0d-0410-9eb7-867de7ffcac7

bazaar/ArrayFire_demo/ArrayFire_demo.cpp

@@ -0,0 +1,36 @@
+#include <Core/Core.h>
+
+using namespace Upp;
+
+#include <arrayfire.h>
+
+using namespace af; 
+
+void MatMult_Bench();
+void Pi_Bench();
+void Vectorize_Bench();
+void Demo();
+
+CONSOLE_APP_MAIN
+{
+    try {
+		int device = 0;
+        setDevice(device);
+        info();
+
+		Pi_Bench();
+		MatMult_Bench();
+		Vectorize_Bench();
+		Demo();
+    } catch (af::exception& e) {
+        printf("\nError: %s\n", e.what());
+    } catch (...) {
+        printf("\nUnknown error\n");
+    }
+    printf("\nEnd");
+	#ifdef WIN32 // pause in Windows
+        printf(". Hit enter...");
+        fflush(stdout);
+        getchar();
+    #endif
+}

bazaar/ArrayFire_demo/ArrayFire_demo.upp

@@ -0,0 +1,23 @@
+uses
+	Core;
+
+library(AF_CPU) afCPU;
+
+library(AF_OPEN_CL) afOpenCL;
+
+library(AF_CUDA) afCUDA;
+
+options
+	;
+
+file
+	ArrayFire_demo.cpp,
+	matrix.cpp,
+	pi.cpp,
+	vectorize.cpp,
+	demo.cpp,
+	srcdoc.tpp;
+
+mainconfig
+	"" = "AF_CUDA";
+

bazaar/ArrayFire_demo/ArrayFire_demo.upp.bak

@@ -0,0 +1,20 @@
+uses
+	Core;
+
+library(AF_CPU) afCPU;
+
+library(AF_CL) afOpenCL;
+
+options
+	;
+
+file
+	ArrayFire_demo.cpp,
+	pi.cpp,
+	progress.h,
+	vectorize.cpp,
+	helloworld.cpp;
+
+mainconfig
+	"" = "";
+

bazaar/ArrayFire_demo/demo.cpp

@@ -0,0 +1,50 @@
+/*******************************************************
+ * 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
+ ********************************************************/
+
+#include <arrayfire.h> 
+
+using namespace af; 
+ 
+void Demo()  { 
+    printf("Create a 5-by-3 matrix of random floats on the GPU\n");
+    array A = randu(5,3, f32);
+    af_print(A);
+
+    printf("Element-wise arithmetic\n");
+    array B = sin(A) + 1.5;
+    af_print(B);
+
+    printf("Negate the first three elements of second column\n");
+    B(seq(0, 2), 1) = B(seq(0, 2), 1) * -1;
+    af_print(B);
+
+    printf("Fourier transform the result\n");
+    array C = fft(B);
+    af_print(C);
+
+    printf("Grab last row\n");
+    array c = C.row(end);
+    af_print(c);
+
+    printf("Create 2-by-3 matrix from host data\n");
+    float d[] = { 1, 2, 3, 4, 5, 6 };
+    array D(2, 3, d, af_source::afHost);
+    af_print(D);
+
+    printf("Copy last column onto first\n");
+    D.col(0) = D.col(end);
+    af_print(D);
+
+    // Sort A
+    printf("Sort A and print sorted array and corresponding indices\n");
+    array vals, inds;
+    sort(vals, inds, A);
+    af_print(vals);
+    af_print(inds);
+}

bazaar/ArrayFire_demo/init

@@ -0,0 +1,4 @@
+#ifndef _ArrayFire_demo_icpp_init_stub
+#define _ArrayFire_demo_icpp_init_stub
+#include "Core/init"
+#endif

bazaar/ArrayFire_demo/matrix.cpp

@@ -0,0 +1,33 @@
+#include <Core/Core.h>
+
+using namespace Upp;
+
+#include <arrayfire.h>
+
+using namespace af; 
+
+static array A, B; // populated before each timing
+
+static void MatMult() {
+    array C = matmul(A, B);  // matrix multiply
+    C.eval();                // ensure evaluated
+}
+
+void MatMult_Bench() {
+    printf("\nBenchmark N-by-N matrix multiply");
+    
+  	double peak = 0;
+    for (int n = 200; n <= 2000; n += 200) {
+        printf("\n%4d x %4d: ", n, n);
+        A = randu(n, n, f32);
+        B = randu(n, n, f32);
+        double time = timeit(MatMult); // time in seconds
+        double gflops = 2.0 * pow(n, 3) / (time * 1e9);
+        if (gflops > peak)
+            peak = gflops;
+
+        printf(" %f secs, %4.0f Gflops", time, gflops);
+        fflush(stdout);
+    }
+    printf("\nPeak %g GFLOPS\n", peak);
+}

bazaar/ArrayFire_demo/pi.cpp

@@ -0,0 +1,44 @@
+#include <Core/Core.h>
+
+using namespace Upp;
+
+#include <arrayfire.h>
+
+using namespace af; 
+
+// generate millions of random samples
+static int samples = 40000000;
+
+/* 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()) / RAND_MAX;
+        float y = float(rand()) / 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(); }
+
+void Pi_Bench() {
+	double t_device = timeit(device_wrapper);
+	double t_host = timeit(host_wrapper);
+	printf("\nPI number benchmark\n");
+  	printf("device:  %.5f seconds to estimate  pi = %.8f\n", t_device, pi_device());
+   	printf("  host:  %.5f seconds to estimate  pi = %.8f\n", t_host, pi_host());
+   	printf("GPU/accelerated device is %.2f times faster than CPU host\n", t_host/t_device);
+}

bazaar/ArrayFire_demo/srcdoc.tpp/ArrayFire$en-us.tpp

@@ -0,0 +1,25 @@
+topic "ArrayFire";
+[ $$0,0#00000000000000000000000000000000:Default]
+[a83;*R6 $$1,0#31310162474203024125188417583966:caption]
+[{_}%EN-US 
+[s1; [+184 ArrayFire]&]
+[s0; [2 ArrayFire package includes simple ][^http`:`/`/arrayfire`.com`/^2 ArrayFire][2  
+benchmarks and demos.]&]
+[s0;2 &]
+[s0; [^http`:`/`/arrayfire`.com`/^2 ArrayFire][2  is a library for fast 
+GPU computing, supporting both Nvidia CUDA and OpenCL devices, 
+and it is open source (BSD 3`-clause). ]&]
+[s0;2 &]
+[s0; [2 It includes matrix algebra, algorithms, image processing, etc. 
+classes and functions.]&]
+[s0;2 &]
+[s0; [2 You can download the sources from ][^https`:`/`/github`.com`/arrayfire`/arrayfire^2 G
+itHub][2  or the binaries from the ][^http`:`/`/arrayfire`.com`/login`/`?redirect`_to`=http`%3A`%2F`%2Farrayfire`.com`%2Fdownload^2 A
+rrayFire][2  page.]&]
+[s0;2 &]
+[s0; [2 Benchmarks and demos are:]&]
+[s0;i150;O0; [2 Pi-|-|Pi number benchmark]&]
+[s0;i150;O0; [2 Matrix-|NxN matrix product benchmark]&]
+[s0;i150;O0; [2 Vectorize-|Different strategies to do operations between 
+vectors]&]
+[s0;i150;O0; [2 Demo-|Basic matrix algebra demos]]]

bazaar/ArrayFire_demo/vectorize.cpp

@@ -0,0 +1,139 @@
+#include <Core/Core.h>
+
+using namespace Upp;
+
+#include <arrayfire.h>
+
+using namespace af; 
+
+
+static array A, B;
+
+static array dist_naive(array a, array b) {
+    array dist_mat = constant(0, a.dims(1), b.dims(1));
+
+    for (int ii = 0; ii < a.dims(1); ii++) 			// Iterate through columns a
+        for (int jj = 0; jj < b.dims(1); jj++) 		// Iterate through columns of b
+            for (int kk = 0; kk < a.dims(0); kk++) 	// Get the sum of absolute differences
+                dist_mat(ii, jj) += abs(a(kk, ii) - b(kk, jj));
+
+    return dist_mat;
+}
+
+static array dist_vec(array a, array b) {
+    array dist_mat = constant(0, a.dims(1), b.dims(1));
+
+    for (int ii = 0; ii < a.dims(1); ii++) {		// Iterate through columns a
+        array avec = a(span, ii);
+
+        for (int jj = 0; jj < b.dims(1); jj++) {	// Iterate through columns of b
+            array bvec = b(span, jj);
+
+            dist_mat(ii, jj) = sum(abs(avec - bvec));// Get SAD using sum on the vector
+        }
+    }
+    return dist_mat;
+}
+
+static array dist_gfor1(array a, array b) {
+    array dist_mat = constant(0, a.dims(1), b.dims(1));
+
+    gfor (seq ii, a.dims(1)) {						// GFOR along columns of a
+        array avec = a(span, ii);
+
+        for (int jj = 0; jj < b.dims(1); jj++) {	// Itere through columns of b
+            array bvec = b(span, jj);
+
+            dist_mat(ii, jj) = sum(abs(avec - bvec));// Get SAD using sum on the vector
+        }
+    }
+    return dist_mat;
+}
+
+static array dist_gfor2(array a, array b) {
+    array dist_mat = constant(0, a.dims(1), b.dims(1));
+
+    gfor (seq jj, b.dims(1)) {						// GFOR along columns of b
+        array bvec = b(span, jj);
+
+        for (int ii = 0; ii < a.dims(1); ii++) {	// Iterate through columns of A
+            array avec = a(span, ii);
+
+            dist_mat(ii, jj) = sum(abs(avec - bvec));// Get SAD using sum on the vector
+        }
+    }
+    return dist_mat;
+}
+
+static array dist_tile1(array a, array b) {
+    int alen = a.dims(1);
+    int blen = b.dims(1);
+
+    array dist_mat = constant(0, alen, blen);
+
+    for (int jj = 0; jj < blen; jj++) {    			// Iterate through columns of b
+        											// Get the column vector of b
+        array bvec = b(span, jj);					// shape of bvec is (feat_len, 1)
+        											// Tile avec to be same size as a
+        array bvec_tiled = tile(bvec, 1, alen);		// shape of bvec_tiled is (feat_len, alen)
+
+        array sad = sum(abs(bvec_tiled - a));		// Get the sum of absolute differences
+
+        											// sad is row vector, dist_mat needs column vector
+        dist_mat(span, jj) = sad.T();				// transpose sad and fill in dist_mat
+    }
+    return dist_mat;
+}
+
+static array dist_tile2(array a, array b) {
+    int feat_len = a.dims(0);
+    int alen = a.dims(1);
+    int blen = b.dims(1);
+    
+    array a_mod = a;								// Shape of a is (feat_len, alen, 1)
+
+    array b_mod = moddims(b, feat_len, 1, blen);	// Reshape b from (feat_len, blen) to (feat_len, 1, blen)
+    
+    array a_tiled = tile(a_mod, 1, 1, blen);		// Tile both matrices to be (feat_len, alen, blen)
+    array b_tiled = tile(b_mod, 1, alen, 1);
+    												// Do The sum operation along first dimension
+    array dist_mod = sum(abs(a_tiled - b_tiled));	// Output is of shape (1, alen, blen)
+
+    array dist_mat = moddims(dist_mod, alen, blen);	// Reshape dist_mat from (1, alen, blen) to (alen, blen)
+    return dist_mat;
+}
+
+static void bench_naive()	{dist_naive(A, B);}
+static void bench_vec()		{dist_vec(A, B);}
+static void bench_gfor1()	{dist_gfor1(A, B);}
+static void bench_gfor2()	{dist_gfor2(A, B);}
+static void bench_tile1()	{dist_tile1(A, B);}
+static void bench_tile2()	{dist_tile2(A, B);}
+
+void Vectorize_Bench() {
+	printf("\nVectorize demo:");
+    // Do not increase the sizes
+    // dist_naive and dist_vec get too slow at large sizes
+    A = randu(3, 200);
+    B = randu(3, 300);
+
+    array d1 = dist_naive(A, B);
+    array d2 = dist_vec  (A, B);
+    array d3 = dist_gfor1(A, B);
+    array d4 = dist_gfor2(A, B);
+    array d5 = dist_tile1(A, B);
+    array d6 = dist_tile2(A, B);
+
+    printf("\nMax. Error for dist_vec  : %f", max<float>(abs(d1 - d2)));
+    printf("\nMax. Error for dist_gfor1: %f", max<float>(abs(d1 - d3)));
+    printf("\nMax. Error for dist_gfor2: %f", max<float>(abs(d1 - d4)));
+    printf("\nMax. Error for dist_tile1: %f", max<float>(abs(d1 - d5)));
+    printf("\nMax. Error for dist_tile2: %f", max<float>(abs(d1 - d6)));
+    printf("\n");
+    printf("\nTime for dist_naive: %2.2fms", 1000 * timeit(bench_naive));
+    printf("\nTime for dist_vec  : %2.2fms", 1000 * timeit(bench_vec  ));
+    printf("\nTime for dist_gfor1: %2.2fms", 1000 * timeit(bench_gfor1));
+    printf("\nTime for dist_gfor2: %2.2fms", 1000 * timeit(bench_gfor2));
+    printf("\nTime for dist_tile1: %2.2fms", 1000 * timeit(bench_tile1));
+    printf("\nTime for dist_tile2: %2.2fms", 1000 * timeit(bench_tile2));
+}