             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();
             return ret[0];

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

            static scalarfield out_min(blocks);
            static scalarfield out_max(blocks);
            static scalarfield temp(1);

            cu_MinMax<<<blocks, threads>>>(&vf[0][0], out_min.data(), out_max.data(), N);
            cu_MinMax<<<1, 1024>>>(out_min.data(), out_min.data(), temp.data(), blocks);

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

            // Dot product
