variants on dp question

dp-m2.chpl

// shows program modified to avoid array views in inner loops

class AbstractKernel {
  // assumes x[xfirst, ..] and y[yfirst, ..] have shape 1..p
  proc kernel(x, xfirst: int, y: [?D] ?T, yfirst: int, p: int): T
  {
    return x[xfirst, 0]; // out of bounds (intentionally?)
  }
}

class DotProduct: AbstractKernel {
  override proc kernel(x, xfirst: int, y: [?D] ?T, yfirst: int, p: int): T
  {
    var ret: T = 0: T;
    for i in 1..p {
      ret += x[xfirst, i] * y[yfirst, i];
    }
    return ret;
  }
}

/***************************************************************************/
use DynamicIters;
proc calculateKernelMatrix(K: AbstractKernel, data: [?D] ?T) /* : [?E] T */
{
  var n = D.dim(0).last;
  var p = D.dim(1).last;
  var E: domain(2) = {D.dim(0), D.dim(0)};
  var mat: [E] T;
  forall j in guided(1..n by -1) {
    for i in j..n {
      mat[i, j] = K.kernel(data, i, data, j, p);
      mat[j, i] = mat[i, j];
    }
  }
  return mat;
}

use Random;
use Time;

const n = 10000; const p = 784;
var x: [1..n, 1..p] real(32);
fillRandom(x);

/**
  To compile:
  chpl kernel.chpl --fast && ./kernel
*/

var K: AbstractKernel = new DotProduct();
var t = new Timer();
t.start();
var mat = calculateKernelMatrix(K, x);
t.stop();
writeln("Time taken (s): ", t.elapsed(TimeUnits.microseconds)/1000000);

dp-m3.chpl

// shows alternative strategy of using row pointers
// note speedup over dp-m2 is probably caused by removing dynamic dispatch, not row pointers

use CPtr;

// one day, compiler will be able to check these
/*
interface AbstractKernel {
  // assumes xrow and yrow have shape 0..#p
  proc kernel(xrow:c_ptr, yrow:c_ptr(?T), p: int): T
  {
    return xrow[0];
  }
}*/

record DotProduct /* implements AbstractKernel */ {
  proc kernel(xrow:c_ptr, yrow:c_ptr(?T), p: int): T
  {
    var ret: T = 0: T;
    for i in 0..#p {
      ret += xrow[i] * yrow[i];
    }
    return ret;
  }
}

/***************************************************************************/
use DynamicIters;
proc calculateKernelMatrix(K, data: [?D] ?T) /* : [?E] T */
{
  var n = D.dim(0).last;
  var p = D.dim(1).last;
  var E: domain(2) = {D.dim(0), D.dim(0)};
  var mat: [E] T;
  // code below assumes data starts at 1,1
  var rowPointers: [1..n] c_ptr(T) =
    forall i in 1..n do c_ptrTo(data[i, 1]);

  forall j in guided(1..n by -1) {
    for i in j..n {
      mat[i, j] = K.kernel(rowPointers[i], rowPointers[j], p);
      mat[j, i] = mat[i, j];
    }
  }
  return mat;
}

use Random;
use Time;

const n = 10000; const p = 784;
var x: [1..n, 1..p] real(32);
fillRandom(x);

/**
  To compile:
  chpl kernel.chpl --fast && ./kernel
*/

var K = new DotProduct();
var t = new Timer();
t.start();
var mat = calculateKernelMatrix(K, x);
t.stop();
writeln("Time taken (s): ", t.elapsed(TimeUnits.microseconds)/1000000);