 	cudaDeviceSynchronize();
 	clock_gettime(CLOCK_MONOTONIC, toc);
 	cudaMemcpy(h_out, d_out, sizeof(double) * numCells * numLeftFields * numRightFields, cudaMemcpyDeviceToHost);
+
+  cudaFree(d_right);
+  cudaFree(d_left);
+  cudaFree(d_out);
 
 }
 

	clock_gettime(CLOCK_MONOTONIC, &toc);
	const double elapsedTime_serial = getElapsedTime(tic, toc);

	std::cout << "trying cuda" << std::endl;
	//Now try the cuda version, start with warmup
	cudaDocontractFieldFieldScalar(cudaOut,cudaLeft,cudaRight, c, l, r, p, t1, t2, &tic, &toc);


    int lbf = matrixIndex / numRightFields;
    int rbf = matrixIndex % numRightFields;
    int sub = dim1Tensor * dim2Tensor;
    int left = myCell * numLeftFields * numPoints * sub + lbf * numPoints * sub;
    int right = myCell * numPoints * sub * numRightFields;
    int rsub = sub * numRightFields;

    double temp = 0;
    for (int qp = 0; qp < numPoints; qp++) {
      for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) {
        for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) {
          temp += d_left[left + qp * sub +
                         iTens1 * dim2Tensor +
                         iTens2] *
                  d_right[right + qp * rsub +
                          iTens1 * dim2Tensor * numRightFields +
                          iTens2 * numRightFields +
                          rbf];
