 
 // Perform the MHA update on the n parameters, done in parallel on the GPU. This is what the kernel does.
 __global__
-void update_chi(double* chi, double* meas, double* meas_unc, int n, double* logdens, double* jump_sigma,
+void update_chi(double* chi, double* meas, double* meas_unc, int n, double* logdens, double* prop_cholfact,
                 curandState* state, int current_iter)
 {
 	int i = blockDim.x * blockIdx.x + threadIdx.x;

    
	for (int i=0; i<n; i++) {
        // Loop over the data indices
		for (int j=0; j<p; j++) {
            // Loop over the chi indices to generate true value of the characteristics
			h_true_chi[i + n * j] = h_theta[j] + sqrt(h_theta[p + j]) * snorm(rng);
            // Initialize the scale of the chi proposal distributions to just be some small constant value
            h_jump_sigma[i + n * j] = 0.01;
        }
        for (int k=0; k<m; k++) {
            // Loop over the feature indices to generate measurements, given the characteristics

    
    // Cholesky factor of proposal distribution for each data point, used by the Metropolis algorithm. The proposal
    // covariance matrix for each data point is Covar_i = PropChol_i * transpose(PropChol_i).
    dim_cholfactor = p * p - ((p - 1) * p) / 2  // cholesky factor is a lower-diagonal matrix
    thrust::host_vector<double> h_prop_cholfact(n * dim_chol_factor);
    thrust::fill(h_prop_cholfact.begin(), h_prop_cholfact.end(), 0.0);
    
	// Create simulated characteristics and data
    thrust::host_vector<double> h_true_chi(n * p);
    double sigma_msmt[5] = {2.2e-4, 3.3e-4, 5.2e-4, 3.7e-4, 2.2e-4};  // Standard deviation for the measurement errors
    double lambda[5] = {70.0, 170.0, 250.0, 350.0, 500.0};  // frequency bands correspond to Herschel PACS and SPIRES wavelengths
    double frequencies[5];
    for (int k=0; k<m; k++) {
        frequencies[k] = 1e6 / lambda[i];
