 	int cell = blockIdx.x;
 	for(int j=0; j<cell; j++)
 	{
-		sample_covariance[cell*N]
+		sample_covariance[cell*N] = 0.0;
 	}
 }
 

	int N = 20; /* Number of Neurons */
	int samps_per_block = 1024;
	int nblocks = 1024

	/* Compute sizes of various data structures */
	size_t mean_size = N*sizeof(float);
	size_t cov_size = N*(N-1)/2*sizeof(float);
	size_t samples_size = N*samps_per_block*nblocks*sizeof(float);

	/* allocate result memory */
	float *alpha_res = (float*)malloc(mean_size);
	float *beta_res = (float*)malloc(cov_size);

	/* Allocate device memory */
	float *sample_mean;

	int i, j;
	int cell;
	float s;
	float df, bdif;
	for(samp=0; samp < nsamps; samp++)
	{
		int sampstart = blk*N*nsamps + samp*N;
		/* need to include nsweeps */
		for(cell=0; cell < N; cell++)
		{

			for(j=0; j<=cell-1; j++)
			{
				bdif += beta[j*N-j*(j+1)/2 + cell];
			}
			for(j=cell+1; j<=N-1; j++)
			{
				bdif += beta[cell*N - cell*(cell+1)/2 +j];
			}
			
			df = -1.0*alpha[cell] - bdif;
			p1 =expf(df)/(1+expf(df));
			if(curand_uniform(&state[blk]) < p1)
			{
