 
 __global__
 void
-cudaDocontractFieldFieldScalar_kernel(double * d_left, double * d_right,
+cudaDocontractFieldFieldScalar_kernel(const double * const __restrict__ d_left, const double * const __restrict__ d_right,
 double * d_out,
 int numCells,
 int numLeftFields,

				cuda_hostRight(cl,rbf, qp) = in_c_r_p(cl,rbf,qp);
				omp_hostRight(cl,rbf,qp) = in_c_r_p(cl,rbf,qp);

				cudaRight[cl * p * r + rbf * p + qp] = in_c_r_p(cl,rbf,qp);
			}
			for(int lbf = 0; lbf < l; ++lbf) {
				cuda_hostLeft(cl, lbf, qp) = in_c_l_p(cl,lbf, qp);
				omp_hostLeft(cl,lbf, qp) = in_c_l_p(cl,lbf,qp);

				cudaLeft[cl * p * l + lbf * p + qp] = in_c_l_p(cl,lbf,qp);
			}
			//cudaRightColMajor[cl + c*qp] = in_r_c_p(cl,qp);
			//cudaLeftColMajor[cl + c*qp] = in_l_c_p(cl,qp);



int main(int argc, char* argv[]) {
	int c=10000, p=10, l=10, r=10;

	FieldContainer<double> in_c_l_p(c, l, p);
	FieldContainer<double> in_c_r_p(c, r, p);
