title: Use lambdas to get polymorphism, recursion, and true OO style in CUDA?
body: <p>I'm trying to migrate a big C++11 library to CUDA. The library builds a hierarchical tree of "Nodes", sharing the same virtual interface but implementing different algorithms, and then calculates the value of the root node for a number of trials. In each trial, the library changes some setting in the leaf nodes of the hierarchy, and then proceeds to recursively recalculate all node values upwards towards the root.</p>
  
  <p>The library has the trial change on the outermost loop, and heavily uses polymorphic virtual classes. It is not feasible to redesign this aspect of it. What I could redesign are the value calculation functions, in order to vectorize them and make them run in kernel space.</p>
  
  <p>Here's a cut down of the original, serial library:</p>
  
  <pre><code>// virtual ancestor
  class Node {
  public:
      // Return value of this node for the current trial.
      // This method will recursively call the value() method
      // of the direct children of this Node.
      virtual float value() = 0;
  };
  
  // Implement y = a1*x1 + a2*x2 + ... + an*xn,
  // where x* are the values  of the direct children of this node
  // and a* are constants
  class LinearCombination: virtual public Node {
  public:
      // skip: constructor
  
      float value() {
          float acc = 0;
          for (size_t i = 0; i &lt; children.size(); i++) {
              // Recursively get value of child
              acc += children[i]-&gt;value() * weights[i];
          }
          return acc;
      }
  
  protected:
      std::vector&lt;Node *&gt; children;
      std::vector&lt;float&gt; weights;
  };
  
  
  int main() {
      Node * root;
      // skip: initialise tree of nodes
      for (auto trial: trials) {
          // skip: setup global trial settings. This will change
          // the result of the leaf Nodes.
          auto value = root-&gt;value();
          // skip: dump to disk
      }
  }
  </code></pre>
  
  <p>My first conversion attempt changes the value() method to calculate the value of each node in a vectorized way, that is for all trials at once. All flow control remains in plain C++, and only the actual maths is moved to the graphics card.
  This is necessary because virtual polymorphism doesn't work inside a kernel, and it's very difficult to implement tree data structures (Node ** children can't be used).</p>
  
  <pre><code>class Node {
  public:
      // Return value for all scenarios
      // result is a buffer in the device memory
      virtual void value(float * result, int n) = 0;
  };
  
  class LinearCombination: virtual public Node {
  public:
      // skip: constructor
  
      void value(float * result, int n) {
          hemi::parallel_for(0, n, [=] HEMI_LAMBDA (int i) {
              result[i] = 0.0;
          });
  
          // Allocate temporary buffer to store the value of the underlyings
          float * scratch;
          cudaMalloc((void **)&amp;scratch, n * sizeof(float));
  
          for (size_t child_id = 0; child_id &lt; children.size(); child_id++) {
              // Recursively get value of child
              children[child_id]-&gt;value(scratch, n);
              auto weight = weights[child_id];
  
              hemi::parallel_for(0, n, [=] HEMI_LAMBDA (int i) {
                  result[i] += scratch[i] * weight;
              });
          }
  
          cudaFree(scratch);
      }
  
  protected:
      std::vector&lt;Node *&gt; children;
      std::vector&lt;float&gt; weights;
  };
  
  
  int main() {
      Node * root;
      // skip: allocate nodes
  
      // Create buffer for output value of the root nodes
      float * value;
      cudaMalloc((void **)&amp;value, n * sizeof(float));
  
      // skip: initialise vectors of data for the leaf nodes
      root-&gt;value(value, n);
      // skip: dump to disk
      cudaFree(value);
  }
  </code></pre>
  
  <p>The above technically works, but its performance is poor - on a GeForce GTX 970, running on 500,000 trials in parallel, it goes only 10x as fast as the serial algorithm goes on a single CPU - put it on a 16-core computer, and the GPU is slower.</p>
  
  <p>This is unsurprising, as in the linear combination example above the value() function performs 3n+1 memory accesses (where n is the number of children), which could be completely avoided if the whole computation were done inside a single kernel.</p>
  
  <p>So I came up with the idea of using the new C++11 Lambda support in CUDA 7.5:</p>
  
  <pre><code>class Node {
  public:
      // Return __device__ lambda which returns the value
      // of the node for a single trial
      virtual std::function&lt;float (int)&gt; valueFunc() = 0;
  };
  
  
  class LinearCombination: virtual public Node {
  public:
      // skip: constructor
  
      std::function&lt;float (int)&gt; valueFunc() {
          auto func = [=] HEMI_LAMBDA (int i) {
              return 0.0;
          };
  
          for (size_t child_id = 0; child_id &lt; children.size(); child_id++) {
              auto childFunc = children[child_id]-&gt;valueFunc();
              auto weight = weights[child_id];
  
              func = [=] HEMI_LAMBDA (int i) {
                  return func(i) + childFunc(i) * weight;
              };
          }
          return func;
      }
  
  protected:
      std::vector&lt;Node *&gt; children;
      std::vector&lt;float&gt; weights;
  };
  
  
  int main() {
      Node * root;
      // skip: allocate nodes
  
      // Create buffer for output value of the root node
      float * scratch;
      cudaMalloc((void **)&amp;scratch, n * sizeof(float));
  
      // skip: initialise vectors of data for the leaf nodes
      auto valueFunc = root-&gt;valueFunc();
      hemi::parallel_for(0, n, [=] HEMI_LAMBDA (int i) {
          scratch[i] = valueFunc(i);
      });
  
      // skip: dump to disk
      cudaFree(scratch);
  }
  </code></pre>
  
  <p>The idea of the above is that there is one big kernel that processes the whole tree, assembled at runtime as a recursion of scalar lambdas, so the whole tree calculation performs ONE memory write, plus whatever input vector the leaf nodes need to read.</p>
  
  <p>However, it doesn't compile, and I can't understand if it's just a matter of syntax or if what I'm trying to do is outright impossible.</p>
  
  <p>If the above can't be fixed, are there any alternative solutions to the problem? As mentioned earlier, refactoring the whole library to be less recursive, less object-oriented, or less based on virtual polymorphism is not an option.</p>