 }
 //
 
-void Multiply(gpu_vector<real>& input, gpu_vector<real>& output, gpu_vector<real>& r) {
-    //    int size = input.size();
-    //    kMultiplyA<<<CONFIG(size)>>>(sorted_pos,  // input
-    //                                 input.data_d, old_vel_node_mpm.data_d,
-    //                                 marker_Fe_hat.data_d,  // input
-    //                                 marker_Fe.data_d,      // input
-    //                                 marker_Fp.data_d,      // input
-    //                                 marker_volume.data_d,  // input
-    //                                 output.data_d);
-    //
-    //    kMultiplyB<<<CONFIG(size)>>>(input, old_vel_node_mpm.data_d, node_mass.data_d, r.data_d, result_array.data_d);
-}
+void Multiply(gpu_vector<real>& input, gpu_vector<real>& output) {
+    int size = input.size();
+    kMultiplyA<<<CONFIG(size)>>>(pos.data_d,  // input
+                                 input.data_d, old_vel_node_mpm.data_d,
+                                 marker_Fe_hat.data_d,  // input
+                                 marker_Fe.data_d,      // input
+                                 marker_Fp.data_d,      // input
+                                 marker_volume.data_d,  // input
+                                 output.data_d);
+
+    kMultiplyB<<<CONFIG(size)>>>(input.data_d, old_vel_node_mpm.data_d, node_mass.data_d, output.data_d);
+}
+
+CUDA_GLOBAL void kSubtract(int size, real* x, real* y) {
+    const int i = blockIdx.x * blockDim.x + threadIdx.x;
+    if (i < size) {
+        // printf("sg: %d %f\n", i, y[i]);
+        y[i] = y[i] - x[i];
+    }
+}
+
 template <bool inner>
 CUDA_GLOBAL void kResetGlobals(int size) {
     if (inner) {

        Mat33 Fe_hat_inv_transpose = InverseTranspose(marker_Fe_hat[p]);

        real dJ = J * InnerProduct(Fe_hat_inv_transpose, delta_F);
        Mat33 dF_inverse_transposed = -Fe_hat_inv_transpose * Transpose(delta_F) * Fe_hat_inv_transpose;
        Mat33 dJF_inverse_transposed = dJ * Fe_hat_inv_transpose + J * dF_inverse_transposed;
        Mat33 RD = Rotational_Derivative(marker_Fe_hat[p], delta_F);

    real data, block_sum;
    const int tid = block_start + threadIdx.x;
    if (tid < num_items) {
        data = mg[tid] * mg[tid];

        block_sum = BlockReduce(temp_storage).Reduce(data, cub::Sum(), num_valid);

        if (threadIdx.x == 0) {
            atomicAdd(&dot_g_proj_norm[0], block_sum);
        }
        // printf("resid [%f %f]\n", mg[tid], dot_g_proj_norm[0]);
    }
}
void MPM_BBSolver(gpu_vector<real>& r, gpu_vector<real>& delta_v) {
    const uint size = r.size();
    gpu_vector<real> dot_g_proj_norm(1);
    ml.resize(size);
    mg.resize(size);
