         scalar dot(const vectorfield & vf1, const vectorfield & vf2)
         {
             int n = vf1.size();
-            static scalarfield sf(n);
+            static scalarfield sf(n, 0);
+            Vectormath::fill(sf, 0);
             scalar ret;
 
             // Dot product

            int blocks = min((N + threads - 1) / threads, 1024);

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

            int blocks = min((N + threads - 1) / threads, 1024);

            static scalarfield ret(1, 0);
            Vectormath::fill(ret, 0);
            cu_sum<<<blocks, threads>>>(sf.data(), ret.data(), N);
            cudaDeviceSynchronize();
            ret[0] = ret[0]/N;
