     kFeHat<<<CONFIG(host_settings.num_mpm_markers)>>>(pos.data_d, marker_Fe.data_d, grid_vel.data_d,
                                                       marker_Fe_hat.data_d);
 
-    kSVD<<<CONFIG(host_settings.num_mpm_markers)>>>(marker_Fe_hat.data_d, PolarR.data_d, PolarS.data_d);
+    // kSVD<<<CONFIG(host_settings.num_mpm_markers)>>>(marker_Fe_hat.data_d, PolarR.data_d, PolarS.data_d);
 
     kApplyForces<<<CONFIG(host_settings.num_mpm_markers)>>>(pos.data_d,            // input
                                                             marker_Fe_hat.data_d,  // input

        }
    }
}
CUDA_GLOBAL void kSVD(Mat33* marker_Fe_hat, Mat33* PolarR, Mat33* PolarS) {
    const int p = blockIdx.x * blockDim.x + threadIdx.x;
    if (p < device_settings.num_mpm_markers) {
        Mat33 U, V, R, S, W;
        real3 E;
        SVD(marker_Fe_hat[p], U, E, V);
        // Perform polar decomposition F = R*S
        R = MultTranspose(U, V);
        S = V * MultTranspose(Mat33(E), V);

        PolarR[p] = R;
        PolarS[p] = S;
    }
}
CUDA_GLOBAL void kMultiplyA(const real3* sorted_pos,  // input
                            const real* v_array,
                            const real* old_vel_node_mpm,

                                                            marker_volume.data_d,  // input
                                                            node_mass.data_d,      // input
                                                            PolarR.data_d,         // output
                                                            PolarS.data_d,         // output
                                                            grid_vel.data_d);      // output

    kRhs<<<CONFIG(host_settings.num_mpm_nodes)>>>(node_mass.data_d, grid_vel.data_d, rhs.data_d);
