       fm = ksn2f[j];
 
       m = j - vstart;
-      index = i*nvirt + m;
 
       d1p = omega_re - (em-en); d1pp = omega_im;
       d2p = omega_re + (em-en); d2pp = omega_im;
       
       alpha = d1p/(d1p*d1p + d1pp*d1pp) - d2p/(d2p*d2p + d2pp*d2pp);
       beta = -d1pp/(d1p*d1p + d1pp*d1pp) + d2pp/(d2p*d2p + d2pp*d2pp);
-      //printf("i = %d, j = %d, alpha = %f, beta = %f\n", i, j, alpha, beta);
-      old_re = nm2v_re[index];
-      old_im = nm2v_im[index];
-
-      nm2v_re[index] = (fn - fm)*(old_re*alpha - old_im*beta);
-      nm2v_im[index] = (fn - fm)*(old_re*beta + old_im*alpha);
+      old_re = nm2v_re[i*nvirt + m];
+      old_im = nm2v_im[i*nvirt + m];
+
+      nm2v_re[i*nvirt + m] = (fn - fm)*(old_re*alpha - old_im*beta);
+      nm2v_im[i*nvirt + m] = (fn - fm)*(old_re*beta + old_im*alpha);
+      //printf("i = %d, j = %d, m = %d, alpha = %f, beta = %f, old_re = %f, old_im = %f, nm2v_re = %f, nm2v_im = %f\n", 
+      //    i, j, m, alpha, beta, old_re, old_im, nm2v_re[index], nm2v_im[index]);
 
       //nm2v = nm2v * (fn-fm) * ( 1.0 / (comega - (em - en)) - 1.0 /(comega + (em - en)) );
     }

  v_dab_d = init_sparse_matrix_csr_gpu_float(v_dab_vals, v_dab_rowPtr, 
                  v_dab_col_ind, v_dab_shape[0], v_dab_shape[1], v_dab_nnz, v_dab_indptr_size);

  checkCudaErrors(cusparseCreate(&handle_cuparse));
  checkCudaErrors(cublasCreate(&handle_cublas));


  int j = blockIdx.y * blockDim.y + threadIdx.y; //nvirt
  float en=0.0, fn=0.0, em=0.0, fm=0.0, old_re, old_im;
  double d1p, d1pp, d2p, d2pp, alpha, beta;
  int m;

  if (i < nfermi)
  {
