         return 2;
     }
 
+    /*
+    Set up parallel random number generations on the GPU
+    */
     curandState* devStates;  // Create state object for random number generators on the GPU
     // Allocate memory on GPU for RNG states
     CUDA_CALL(cudaMalloc((void **)&devStates, nThreads.x * nBlocks.x *

        double* p_meas = thrust::raw_pointer_cast(&d_meas[0]);
        double* p_meas_unc = thrust::raw_pointer_cast(&d_meas_unc[0]);
        double* p_chi = thrust::raw_pointer_cast(&d_chi[0]);
        double* p_prop_cholfact = thrust::raw_pointer_cast(&d_prop_cholfact[0]);
        double* p_logdens = thrust::raw_pointer_cast(&d_logdens[0]);
        int current_iter = i + 1;
        update_chi<<<nBlocks,nThreads>>>(p_chi, p_meas, p_meas_unc, ndata, p_logdens, p_prop_cholfact, devStates, current_iter);

    h_theta[4] = log(sig_beta * sig_beta);
    h_theta[5] = log(sig_temp * sig_temp);

    // Allocate memory for arrays of chi, measurement, and measurement uncertainties on host
    thrust::host_vector<double> h_meas(ndata * mfeat);  // The measurements, mfeat values for each of ndata data points
	thrust::host_vector<double> h_meas_unc(ndata * mfeat);  // The measurement uncertainties
    thrust::host_vector<double> h_chi(ndata * p);  // Unknown characteristics, p values for each of ndata data points
    thrust::host_vector<double> h_true_chi(ndata * p);

    // Generate characteristics and data
    thrust::host_vector<double> h_true_chi(ndata * p); // true values of the chi-values
    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

    std::cout << "Generating some data..." << std::endl;
    generate_data(ndata, h_theta, h_true_chi, h_meas, h_meas_unc);

    // Allocate memory for cholesky factors on host.
    //
    // Cholesky factor of chi proposal distribution for each data point, used by the Metropolis algorithm. The proposal
    // covariance matrix of chi for each data point is Covar_i = PropChol_i * transpose(PropChol_i).
    thrust::host_vector<double> h_chi_cholfactor(ndata * dim_cholfactor);
    thrust::fill(h_chi_cholfactor.begin(), h_chi_cholfactor.end(), 0.0);

    // Cholesky factor for the covariance matrix of the theta proposal distribution
    int dim_cholfactor_theta = dim_theta * dim_theta - ((dim_theta - 1) * dim_theta) / 2;
    thrust::host_vector<double> h_theta_cholfact(dim_cholfactor_theta);

    // Initialize the chi values, as well as the cholesky factors for the theta and chi proposals
    initialize_parameters(ndata, h_theta, h_chi, h_theta_cholfact, h_chi_cholfactor, h_true_chi);


    // Allocate memory for arrays on device and copy the values from the host
    thrust::device_vector<double> d_meas = h_meas;
    thrust::device_vector<double> d_meas_unc = h_meas_unc;
    thrust::device_vector<double> d_chi = h_chi;
    thrust::device_vector<double> d_chi_cholfactor = h_chi_cholfactor;
    thrust::device_vector<double> d_logdens(ndata);  // log posteriors for an individual data point

    // Initialize the conditional log-posterior of chi|theta
    thrust::fill(d_logdens.begin(), d_logdens.end(), -1.0e300);

	// Load the initial guess for theta into constant memory
    double* p_theta = &h_theta[0];
	CUDA_CALL(cudaMemcpyToSymbol(c_theta, p_theta, h_theta.size() * sizeof(*p_theta)));

