 
         // Finally, adapt the Cholesky factor of the proposal distribution
         // TODO: check that the ratio is finite before doing this step
+        double unit_norm = 0.0;
+        for (int j=0; j<c_p; j++) {
+            unit_norm += snorm_deviate[j] * snorm_deviate[j];
+        }
+        unit_norm = sqrt(unitnorm);
         double decay_sequence = 1.0 / pow(current_iter, c_decay_rate);
-        
+        scaled_coef = sqrt(decay_sequence * fabs(ratio - c_target_rate)) / unit_norm;
+        for (int j=0; j<c_p; j++) {
+            scaled_proposal[j] *= scaled_coef;
+        }
+        bool downdate = (ratio < c_target_rate);
+        CholUpdateR1(cholfactor, scaled_proposal, downdate);
+        
+        // copy cholesky factor for this data point back to global memory
+        cholfact_index = 0
+        for (int j=0; j<dim_cholfactor; j++) {
+            prop_cholfact[n * j + i] = cholfactor[j];
+        }
 	}
 }
 

    // 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

{
    // Random number generator and distribution needed to simulate some data
	// unsigned int seed = 98724732;
    unsigned int seed = 123455;
	static thrust::minstd_rand rng(seed);
	thrust::random::experimental::normal_distribution<double> snorm(0., 1.0);
    thrust::random::uniform_real_distribution<double> uniform(0.0, 1.0);
