     thrust::transform(forces.begin(), forces.end(), pts.begin(), forces.begin(), thrust::multiply<float>());
 
     // compute forces_i = forces_i - sum_j pij*qij*normalization*yj
-    alpha = -4.0f;
+    alpha = -1.0f;
     beta = 1.0f;
     cusparseSafeCall(cusparseScsrmm(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
-                            N, 2, N, N * L, &alpha, descr,
+                            N, 2, N, nnz, &alpha, descr,
                             thrust::raw_pointer_cast(forceProd.data()),
                             thrust::raw_pointer_cast(pijRowPtr.data()),
                             thrust::raw_pointer_cast(pijColInd.data()),

    Math::max_norm(d_knn_distances);

    // Compute the perplexity/pij of the KNN distribution

    // TODO: make sure pij calc works with sparse matrix (should work now?)

    // TODO: symmetrize pij so that it is stored in sparse csr format

    for (int step = 0; step < n_iter; step++) {
        // compute repulsive forces and normalization
        // TODO: allocate all these things
        BoundingBoxKernel<<<blocks * FACTOR1, THREADS1>>>(nnodes, nbodies, startl, childl, massl, posxl, posyl, maxxl, maxyl, minxl, minyl);
        ClearKernel1<<<blocks * 1, 1024>>>(nnodes, nbodies, childl);
        TreeBuildingKernel<<<blocks * FACTOR2, THREADS2>>>(nnodes, nbodies, errl, childl, posxl, posyl);
        ClearKernel2<<<blocks * 1, 1024>>>(nnodes, startl, massl);
        SummarizationKernel<<<blocks * FACTOR3, THREADS3>>>(nnodes, nbodies, countl, childl, massl, posxl, posyl);
        SortKernel<<<blocks * FACTOR4, THREADS4>>>(nnodes, nbodies, sortl, countl, startl, childl);
        ForceCalculationKernel<<<blocks * FACTOR5, THREADS5>>>(nnodes, nbodies, errl, dthf, itolsq, epssq, sortl, childl, massl, posxl, posyl, velxl, velyl, norml);
        // non-normalized xrep stored in velxl, yrep stored in velyl, norm_i stored in norml
        float norm = thrust::reduce(norml.begin(), norml.end(), 0.0f, thrust::plus<float>());

        // TODO: Write a kernel that does elementwise multiplication of qij * pij * (point_i - point_j)
        // This is O(N) because pij is sparse so if we iterate over pij, get i and j, compute qij etc., this could be written as a thrust transform?
        // Should result in an N x 2 force vector (possibly, we could do it inplace with velxl and velyl and just sum into those)
        
        // Finally, normalize the result (Need to check exactly where normalization happens... probably look at bh_tsne code)


                        thrust::device_vector<float> &pts,
                        thrust::device_vector<float> &forces)
{
    assert(nnz == 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(nnz, BLOCKSIZE);
    computePijxQij<<<NBLOCKS, BLOCKSIZE>>>(N, nnz,
                                            thrust::raw_pointer_cast(sparsePij.data()),
                                            thrust::raw_pointer_cast(pijRowPtr.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 = 1.0f;
    float beta = 0.0f;
    thrust::device_vector<float> ones(N*2);
    cusparseSafeCall(cusparseScsrmm(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                            N, 2, N, nnz, &alpha, descr,
                            thrust::raw_pointer_cast(forceProd.data()),
                            thrust::raw_pointer_cast(pijRowPtr.data()),
                            thrust::raw_pointer_cast(pijColInd.data()),
