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

	int i, j;
	int cell;
	float s;
	float up, dp, df, p1;
	for(samp=0; samp < nsamps; samp++)
	{

		int sampstart = blk*N*nsamps + samp*N;
		/* need to include nsweeps */
		for(cell=0; cell < N; cell++)
		{
			samples[sampstart + cell] = 1.0;
			for(i = 0; i < N; i++)
			{
				up += -1.0*alpha[i]*samples[sampstart + i];
				for(j=0; j<i; j++)
				{
					up += -1.0*beta[(i + j)-1]*samples[sampstart+ i]*samples[sampstart + j];
				}
			}
			samples[sampstart + cell] = 0.0;
			for(i = 0; i < N; i++)
			{
				dp += -1.0*alpha[i]*samples[sampstart+ i];
				for(j=0; j<i; j++)
				{
					dp += -1.0*beta[(i + j)-1]*samples[sampstart + i]*samples[sampstart + j];
				}
			}
			df = up - dp;
			p1 =expf(df)/(1+expf(df));
			if(curand_uniform(&state[blk]) < p1)
			{

	curandState *d_state;
	cudaMalloc(&d_state, nblocks);

	/* generate random initial conditoins */
	int i;
	for(i=0; i<N; i++)
	{
		alpha_res[i] = 2*rand()-1;
	}
	for(i=0; i<N*(N-1)/2; i++)
	{
		beta_res[i]=2*rand()-1;
	}

	/* copy initial conditions over to device */
	cudaMemcpy(alpha, alpha_res, mean_size, cudaMemcpyHostToDevice);
	cudaMemcpy(beta, beta_res, cov_size, cudaMemcpyHostToDevice);

	/* sample */
	printf("Sampling...\n");
	init_rand<<<nblocks, 1>>>(d_state);
	gibbs_sample<<<nblocks, 1>>>(d_state, alpha, beta, samples, N, samps_per_block, nsweeps)
	printf("Finished.  nsamps=%d", nsamps_tot);
	printf("Computing Sample mean...");
	compute_sample_mean<<<N, 1>>>(samples, sample_mean, nsamps_tot, N);

	/*copy back sample mean*/
	cudaMemcpy(sample_mean_res, sample_mean, cudaMemcpyDeviceToHost);

	/* Display*/
	for(i=0; i<N; i++)
	{
		printf("%f\n", sample_mean_res[i]);
	}


	/* free memory */
	free(alpha_res);
	free(beta_res);
	free(sample_mean_res);

	cudaFree(sample_mean);
	cudaFree(sample_covariance);
	cudaFree(samples);
	cudaFree(alpha);
	cudaFree(beta);
	cudaFree(d_state);

	return 0;
}
