       if (vely == NULL) {fprintf(stderr, "cannot allocate vely\n");  exit(-1);}
 
       if (cudaSuccess != cudaMalloc((void **)&errl, sizeof(int))) fprintf(stderr, "could not allocate errd\n");  CudaTest("couldn't allocate errd");
-      if (cudaSuccess != cudaMalloc((void **)&childl, sizeof(int) * (nnodes+1) * 8)) fprintf(stderr, "could not allocate childd\n");  CudaTest("couldn't allocate childd");
+      if (cudaSuccess != cudaMalloc((void **)&childl, sizeof(int) * (nnodes+1) * 4)) fprintf(stderr, "could not allocate childd\n");  CudaTest("couldn't allocate childd");
       if (cudaSuccess != cudaMalloc((void **)&massl, sizeof(float) * (nnodes+1))) fprintf(stderr, "could not allocate massd\n");  CudaTest("couldn't allocate massd");
       if (cudaSuccess != cudaMalloc((void **)&posxl, sizeof(float) * (nnodes+1))) fprintf(stderr, "could not allocate posxd\n");  CudaTest("couldn't allocate posxd");
       if (cudaSuccess != cudaMalloc((void **)&posyl, sizeof(float) * (nnodes+1))) fprintf(stderr, "could not allocate posyd\n");  CudaTest("couldn't allocate posyd");

                        thrust::device_vector<float> &pts,
                        thrust::device_vector<float> &forces)
{
    assert(N * L == sparsePij.size());
    assert(pijRowPtr.size() == N + 1);
    assert(pijColInd.size() == sparsePij.size());
    assert(forceProd.size() == sparsePij.size());

    const int BLOCKSIZE = 128;
    const int NBLOCKS = iDivUp(N * L, BLOCKSIZE);
    computePijxQij<<<NBLOCKS, BLOCKSIZE>>>(N, L,
                                            thrust::raw_pointer_cast(sparsePij.data()),
                                            thrust::raw_pointer_cast(pijColInd.data()),
                                            thrust::raw_pointer_cast(forceProd.data()),
                                            thrust::raw_pointer_cast(pts.data()));

    // compute forces_i = sum_j pij*qij*normalization*yi
    float alpha = 4.0f;
    float beta = 0.0f;
    thrust::device_vector<float> ones(N*2);
    cusparseSafeCall(cusparseScsrmm(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                            N, 2, N, N * L, &alpha, descr,
                            thrust::raw_pointer_cast(forceProd.data()),
                            thrust::raw_pointer_cast(pijRowPtr.data()),
                            thrust::raw_pointer_cast(pijColInd.data()),

    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,
                                            float* points, 
                                            unsigned int N_POINTS, 
                                            unsigned int N_DIMS, 
