     while ((nnodes & (WARPSIZE-1)) != 0) nnodes++;
     nnodes--;
 
+    thrust::device_vector<int> errl(1);
     thrust::device_vector<int> startl(nnodes + 1);
     thrust::device_vector<int> childl((nnodes + 1) * 4);
-    thrust::device_vector<float> massl(nnodes + 1); // TODO: probably don't need massl
+    thrust::device_vector<float> massl(nnodes + 1, 1); // TODO: probably don't need massl
     thrust::device_vector<int> countl(nnodes + 1);
     thrust::device_vector<int> startl(nnodes + 1);
     thrust::device_vector<int> sortl(nnodes + 1);
     thrust::device_vector<float> norml(nnodes + 1);
+    float eta = 10.0f;
+    float norm;
+    // These variables currently govern the tolerance (whether it recurses on a cell)
+    float epssq = 0.05 * 0.05;
+    float itolsq = 1.0f / (0.5 * 0.5);
 
     for (int step = 0; step < n_iter; step++) {
         // compute attractive forces

    // divide column by resulting vector
    Broadcast::broadcast_matrix_vector(pij, sums, K, N, thrust::divides<float>(), 1, 1.0f);
}

thrust::device_vector<float> BHTSNE::tsne(cublasHandle_t &dense_handle, 
                                          cusparseHandle_t &sparse_handle,

      CudaTest("kernel 4 launch failed");

      cudaEventRecord(start, 0);
      ForceCalculationKernel<<<blocks * FACTOR5, THREADS5>>>(nnodes, nbodies, errl, itolsq, epssq, sortl, childl, massl, posxl, posyl, velxl, velyl, norml);
      cudaEventRecord(stop, 0);  cudaEventSynchronize(stop);  cudaEventElapsedTime(&time, start, stop);
      timing[5] += time;
      CudaTest("kernel 5 launch failed");
