 
 		double temp = 0;
 		for (int qp = 0; qp < numPoints; qp++) {
-			temp += d_left[myMatrix*numPoints*numLeftFields + qp*numLeftFields + matrixRow] *
+			temp += d_left[myMatrix*numPoints*numLeftFields + numPoints*matrixRow + qp] *
 							d_right[myMatrix*numPoints*numRightFields + qp*numRightFields + matrixCol];
 		}
 		d_out[myID]= temp;


	dim3 blockSize(1024);
	dim3 gridSize((numCells * numLeftFields * numRightFields / 1024) + 1);

	cudaDocontractFieldFieldScalar_kernel<<<gridSize, blockSize>>>(d_left,
			d_right, d_out, numCells, numLeftFields, numRightFields, numPoints);

	cudaMemcpy(h_out, d_out, sizeof(double) * numCells * numLeftFields * numRightFields, cudaMemcpyDeviceToHost);

}

		int numCells,
		int numLeftFields,
		int numRightFields,
		int numPoints,
		timespec * tic,
		timespec * toc) {

	double * d_right;
	double * d_left;
