 	dim3 gridSize((numCells * numLeftFields * numRightFields / 1024) + 1);
 
 	cudaDocontractFieldFieldScalar_kernel<<<gridSize, blockSize>>>(d_left,
-			d_right, d_out, numCells, numLeftCells, numRightCells, numPoints);
-
-	cudaMemcpy(h_out, d_out, sizeof(double) * numCells * numLeftCells,
-									numRightCells, cudaMemcpyDeviceToHost);
+			d_right, d_out, numCells, numLeftFields, numRightFields, numPoints);
+
+	cudaMemcpy(h_out, d_out, sizeof(double) * numCells * numLeftFields * numRightFields, cudaMemcpyDeviceToHost);
 
 }
 /*

		int matrixRow = matrixIndex / numRightFields;
		int matrixCol = matrixIndex % numRightFields;

		tmpVal += leftFields(cl, lbf, qp)*rightFields(cl, rbf, qp);
		double temp = 0;
		for (int qp = 0; qp < numPoints; qp++) {
			temp += d_left[myMatrix*numPoints*numLeftFields + matrixRow*numPoints + qp] *

	for (int cl = 0; cl < c; ++cl) {
		for(int lbf = 0; lbf < l; ++lbf) {
			for(int rbf = 0; rbf < r; ++rbf) {
				out1_c_l_r(cl,lbf,rbf) = cudaOut[cl * l * r + lbf * r + rbf];
			}
		}
