             int idx = blockIdx.x * blockDim.x + threadIdx.x;
             if(idx < N)
             {
-                atomicExch(&sf[idx], s*sf[idx]);
+                sf[idx] *= s;
             }
         }
 		void scale(scalarfield & sf, scalar s)

            cudaDeviceSynchronize();
		}

        __global__ void cu_sum(const scalar *sf, scalar *out, size_t N)
        {
            int idx = blockIdx.x * blockDim.x + threadIdx.x;
            if(idx < N)
            {
                atomicAdd(out, sf[idx]);
            }
        }
		scalar sum(const scalarfield & sf)
		{
            scalarfield ret(1, 0);
            int n = sf.size();
            cu_sum<<<(n+1023)/1024, 1024>>>(sf.data(), ret.data(), n);
            cudaDeviceSynchronize();
            return ret[0];
		}

		scalar mean(const scalarfield & sf)
		{
            scalarfield ret(1, 0);
            int n = sf.size();
            cu_sum<<<(n+1023)/1024, 1024>>>(sf.data(), ret.data(), n);
            cudaDeviceSynchronize();
            ret[0] = ret[0]/n;
            return ret[0];
		}


            cudaDeviceSynchronize();
        }

        Vector3 sum(const vectorfield & vf)
		{
            int N = vf.size();
            int threads = 512;
            int blocks = min((N + threads - 1) / threads, 1024);

            vectorfield ret(1, {0,0,0});
            cu_sum<<<blocks, threads>>>(vf.data(), ret.data(), N);
            cudaDeviceSynchronize();
            return ret[0];
		}

		Vector3 mean(const vectorfield & vf)
		{
            int N = vf.size();
            int threads = 512;
            int blocks = min((N + threads - 1) / threads, 1024);

            vectorfield ret(1, {0,0,0});
            cu_sum<<<blocks, threads>>>(vf.data(), ret.data(), N);
            cudaDeviceSynchronize();
            ret[0] = ret[0]/N;
            return ret[0];
		}


