 		d_right, d_out, numCells,numPoints);
 
 	cudaMemcpy(h_out, d_out, sizeof(double) * numCells, cudaMemcpyDeviceToHost);
+
+  cudaFree(d_right);
+  cudaFree(d_left);
+  cudaFree(d_out);
 
 }
 

			_outputFields(elementIndex) = tmpVal;
		}
};




template<class DeviceType, class LeftViewType, class RightViewType, class OutputViewType>
struct contractDataDataScalarFunctor1D {
	typedef DeviceType device_type;
	LeftViewType _leftFields;
	RightViewType _rightFields;
	OutputViewType _outputFields;
	int _numLeftFields;
	int _numRightFields;
	int _numPoints;
	int _dim1Tensor;
	int _dim2Tensor;
	int _numCells;

	contractDataDataScalarFunctor1D(LeftViewType leftFields,
			RightViewType rightFields,
			OutputViewType outputFields,
			int numLeftFields,
			int numRightFields,
			int numPoints,
			int dim1Tensor,
			int dim2Tensor,
			int numCells) :
		_leftFields(leftFields),
		_rightFields(rightFields),
		_outputFields(outputFields),
		_numLeftFields(numLeftFields),
		_numRightFields(numRightFields),
		_numPoints(numPoints),
		_dim1Tensor(dim1Tensor),
		_dim2Tensor(dim2Tensor),
		_numCells(numCells)
	{
		// Nothing to do
	}

	KOKKOS_INLINE_FUNCTION
		void operator()(const unsigned int elementIndex) const {

			for (int lbf = 0; lbf < _numLeftFields; lbf++) {
				for (int rbf = 0; rbf < _numRightFields; rbf++) {
					double tmpVal = 0;
					for (int qp = 0; qp < _numPoints; qp++) {
						for (int iTens1 = 0; iTens1 < _dim1Tensor; iTens1++) {
							for (int iTens2 = 0; iTens2 < _dim2Tensor; iTens2++) {
								tmpVal +=
									_leftFields(lbf*_numPoints*_dim1Tensor*_dim2Tensor*_numCells +
											qp*_dim1Tensor*_dim2Tensor*_numCells +
											iTens1*_dim2Tensor*_numCells + iTens2*_numCells + elementIndex)
									*_rightFields(rbf*_numPoints*_dim1Tensor*_dim2Tensor*_numCells +
											qp*_dim1Tensor*_dim2Tensor*_numCells +
											iTens1*_dim2Tensor*_numCells + iTens2*_numCells + elementIndex);
							} // D2-loop
						} // D1-loop
					} // P-loop
					_outputFields(lbf*_numRightFields*_numCells +
							rbf*_numCells + elementIndex) = tmpVal;
				} // R-loop
			} // L-loop
		}
};





// Serial contractDataDataScalar.  Contracts FieldContainers of doubles.

	clock_gettime(CLOCK_MONOTONIC, &tic);

	//repeat the calculation 5 times so we can average out some randomness
	for(int i = 0; i < numRepeats; ++i){
		contractDataDataScalarSerial(out2_c, inl_c_p, inr_c_p);
	}

