     dy = iy - jy;
     tmp = 1 / (1 + dx*dx + dy*dy);
     forceProd[TID] = pij[TID] * tmp;
-
 }
 
 // computes unnormalized attractive forces
 void computeAttrForce(int N,
-                        int L,
+                        int nnz,
                         cusparseHandle_t &handle,
                         cusparseMatDescr_t &descr,
                         thrust::device_vector<float> &sparsePij,

/******************************************************************************/
__global__
void computePijxQij(int N, 
                    int L, 
                    volatile float * __restrict pij,
                    volatile int   * __restrict pijInds,
                    volatile float * __restrict forceProd,
                    volatile float * __restrict pts)
{
    register int TID, i, j;
    register float ix, iy, jx, jy, dx, dy, tmp;
    TID = threadIdx.x + blockIdx.x * blockDim.x;
    if (TID > N * L) return;
    

    i = TID / L;
    j = TID % L;
    j = pijInds[i*L + j];
    
    ix = pts[i]; iy = pts[N + i];
    jx = pts[j]; jy = pts[N + j];

    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, 
