 		}
 		#endif
 		//pAp=p'*Ap
-		temp = p[blockIdx.x*blockDim.x + threadIdx.x]
+		temp = sharedx[threadIdx.x]
 				*ap[blockIdx.x*blockDim.x + threadIdx.x];
 		temp = blockReduceSum(temp);
 		if(threadIdx.x == 0){

	float temp = 0;
	for(int i = 0; i < f; i++)
		//this is math correct and coalesced because A is symmetric
		temp += A[blockIdx.x*f*f + f*i + threadIdx.x]*x[blockIdx.x*f + i];
	r[blockIdx.x*blockDim.x + threadIdx.x] = b[blockIdx.x*blockDim.x + threadIdx.x] - temp;
	//p=r;
	p[blockIdx.x*blockDim.x + threadIdx.x] = r[blockIdx.x*blockDim.x + threadIdx.x];


	for(int iter = 0; iter < cgIter; iter++){
		//ap=A*p;
		//sharedx <-- p
		sharedx[threadIdx.x] = p[blockIdx.x*blockDim.x + threadIdx.x];
		__syncthreads();
		//WARN: set temp to zero since the next operation is +=!
		temp = 0;
		for(int i = 0; i < f; i++)
			//this is math correct and coalesced because A is symmetric
			temp += A[blockIdx.x*f*f + f*i + threadIdx.x]*sharedx[i];
		ap[blockIdx.x*blockDim.x + threadIdx.x] = temp;
		__syncthreads();
		#ifdef DEBUG
