             double proposed_var_j = exp(proposed_theta[j + p]);
             // transform and reduction is over d_chi[j * n : (j+1) * n - 1]
             zsqr zsqrj(proposed_mu_j, proposed_var_j);
-            logdens_pop = thrust::transform_reduce(chi_iter_begin, chi_iter_begin+n,
+            logdens_pop = thrust::transform_reduce(chi_iter_begin, chi_iter_begin+ndata,
                                                 zsqrj, logdens_pop, thrust::plus<double>());
-            thrust::advance(chi_iter_begin, n);
+            thrust::advance(chi_iter_begin, ndata);
         }
 
         double logdens_old = thrust::reduce(d_logdens.begin(), d_logdens.end());

    thrust::device_vector<double> d_meas_unc = h_meas_unc;
    thrust::device_vector<double> d_chi = h_chi;
    thrust::device_vector<double> d_prop_cholfact = h_prop_cholfact;
    thrust::device_vector<double> d_logdens(n);  // log posteriors for an individual data point

    thrust::fill(d_logdens.begin(), d_logdens.end(), -1.0e300);



    // Cuda grid launch
    dim3 nThreads(256);
    dim3 nBlocks((ndata + nThreads.x-1) / nThreads.x);
    printf("nBlocks: %d\n", nBlocks.x);  // no more than 64k blocks!
    if (nBlocks.x > 65535)
    {
