                         thrust::device_vector<float> &pts,
                         thrust::device_vector<float> &forces)
 {
-    assert(N * L == sparsePij.size());
+    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(N * L, BLOCKSIZE);
-    computePijxQij<<<NBLOCKS, BLOCKSIZE>>>(N, L,
+    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 = 4.0f;
+    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, 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()),

      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 **)&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");

/******************************************************************************/
__global__
void computePijxQij(int N, 
                    int nnz, 
                    volatile float * __restrict pij,
                    volatile int   * __restrict pijRowPtr,
                    volatile int   * __restrict pijColInd,
                    volatile float * __restrict forceProd,
                    volatile float * __restrict pts)
{
    register int TID, i, j, start, end;
    register float ix, iy, jx, jy, dx, dy, tmp;
    TID = threadIdx.x + blockIdx.x * blockDim.x;
    if (TID > nnz) return;
    start = 0; end = N + 1;
    i = (N + 1) >> 1;
    while (!(pijRowPtr[i] <= TID && pijRowPtr[i+1] > TID)) {
      if (pijRowPtr[i] > TID) {
        end = i;
      } else {
        start = i;
      }
      i = (start + end) >> 1;
    }
    
    j = pijColInd[TID - i];
    
    ix = pts[i]; iy = pts[N + i];
    jx = pts[j]; jy = pts[N + j];
