getting_started/vectorize.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:
* https://arrayfire.com/licenses/BSD-3-Clause
********************************************************/
#include <arrayfire.h>
#include <stdio.h>
#include <af/util.h>
using namespace af;
array A, B;
static array dist_naive(array a, array b) {
array dist_mat = constant(0, a.dims(1), (int)b.dims(1));
// Iterate through columns a
for (int ii = 0; ii < (int)a.dims(1); ii++) {
// Iterate through columns of b
for (int jj = 0; jj < (int)b.dims(1); jj++) {
// Get the sum of absolute differences
for (int kk = 0; kk < (int)a.dims(0); kk++) {
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, (int)a.dims(1), (int)b.dims(1));
// Iterate through columns a
for (int ii = 0; ii < (int)a.dims(1); ii++) {
array avec = a(span, ii);
// Iterate through columns of b
for (int jj = 0; jj < (int)b.dims(1); jj++) {
array bvec = b(span, jj);
// get SAD using sum on the vector
dist_mat(ii, jj) = sum(abs(avec - bvec));
}
}
return dist_mat;
}
static array dist_gfor1(array a, array b) {
array dist_mat = constant(0, (int)a.dims(1), (int)b.dims(1));
// GFOR along columns of a
gfor(seq ii, (int)a.dims(1)) {
array avec = a(span, ii);
// Itere through columns of b
for (int jj = 0; jj < (int)b.dims(1); jj++) {
array bvec = b(span, jj);
// get SAD using sum on the vector
dist_mat(ii, jj) = sum(abs(avec - bvec));
}
}
return dist_mat;
}
static array dist_gfor2(array a, array b) {
array dist_mat = constant(0, (int)a.dims(1), (int)b.dims(1));
// GFOR along columns of b
gfor(seq jj, (int)b.dims(1)) {
array bvec = b(span, jj);
// Iterate through columns of A
for (int ii = 0; ii < (int)a.dims(1); ii++) {
array avec = a(span, ii);
// get SAD using sum on the vector
dist_mat(ii, jj) = sum(abs(avec - bvec));
}
}
return dist_mat;
}
static array dist_tile1(array a, array b) {
// int feat_len = (int)a.dims(0); // Same as (int)b.dims(0);
int alen = (int)a.dims(1);
int blen = (int)b.dims(1);
array dist_mat = constant(0, alen, blen);
// Iterate through columns of b
for (int jj = 0; jj < blen; jj++) {
// Get the column vector of b
// shape of bvec is (feat_len, 1)
array bvec = b(span, jj);
// Tile avec to be same size as a
// shape of bvec_tiled is (feat_len, alen)
array bvec_tiled = tile(bvec, 1, alen);
// Get the sum of absolute differences
array sad = sum(abs(bvec_tiled - a));
// sad is row vector, dist_mat needs column vector
// transpose sad and fill in dist_mat
dist_mat(span, jj) = sad.T();
}
return dist_mat;
}
static array dist_tile2(array a, array b) {
int feat_len = (int)a.dims(0);
int alen = (int)a.dims(1);
int blen = (int)b.dims(1);
// Shape of a is (feat_len, alen, 1)
array a_mod = a;
// Reshape b from (feat_len, blen) to (feat_len, 1, blen)
array b_mod = moddims(b, feat_len, 1, blen);
// Tile both matrices to be (feat_len, alen, blen)
array a_tiled = tile(a_mod, 1, 1, blen);
array b_tiled = tile(b_mod, 1, alen, 1);
// Do The sum operation along first dimension
// Output is of shape (1, alen, blen)
array dist_mod = sum(abs(a_tiled - b_tiled));
// Reshape dist_mat from (1, alen, blen) to (alen, blen)
array dist_mat = moddims(dist_mod, 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); }
int main(int, char **) {
try {
// 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("Max. Error for dist_vec : %f\n", max<float>(abs(d1 - d2)));
printf("Max. Error for dist_gfor1: %f\n", max<float>(abs(d1 - d3)));
printf("Max. Error for dist_gfor2: %f\n", max<float>(abs(d1 - d4)));
printf("Max. Error for dist_tile1: %f\n", max<float>(abs(d1 - d5)));
printf("Max. Error for dist_tile2: %f\n", max<float>(abs(d1 - d6)));
printf("\n");
printf("Time for dist_naive: %2.2fms\n", 1000 * timeit(bench_naive));
printf("Time for dist_vec : %2.2fms\n", 1000 * timeit(bench_vec));
printf("Time for dist_gfor1: %2.2fms\n", 1000 * timeit(bench_gfor1));
printf("Time for dist_gfor2: %2.2fms\n", 1000 * timeit(bench_gfor2));
printf("Time for dist_tile1: %2.2fms\n", 1000 * timeit(bench_tile1));
printf("Time for dist_tile2: %2.2fms\n", 1000 * timeit(bench_tile2));
} catch (const af::exception &ex) {
fprintf(stderr, "%s\n", ex.what());
throw;
}
return 0;
}