 }
 
 uint ChSolverAPGD::SolveAPGD(custom_vector<real> &x, const custom_vector<real> &b, const uint max_iter) {
-    //real gdiff = 0.000001;
+    real gdiff = 0.000001;
+    real lastgoodres=10e30;
     real theta_k=1.0;
-    custom_vector<real> ms, mg_tmp2, mb_tmp, x01(x.size()), d01;
+    custom_vector<real> ms, mg_tmp2, mb_tmp(x.size()), x01(x.size()), d01(x.size());
     custom_vector<real> mg = ShurProduct(x) - b;
 
     Thrust_Fill(x01, 1.0);

    real t_k = 1 / L_k;
    custom_vector<real>  my = x;
    custom_vector<real>  mx = x;
    for (current_iteration = 0; current_iteration < max_iter; current_iteration++) {
        mg = ShurProduct(my) - b;
        SEAXPY(-t_k, mg, my, mx);   // mx = my + mg*(-t_k);

        t_k = 1 / L_k;
        x = mx;
        theta_k = theta_k1;
        
        SEAXPY(-gdiff, mg_tmp2, x, mb_tmp); //mb_tmp=x+mg_tmp2*(-gdiff)
        Project(mb_tmp);
        d01=mb_tmp-x;
        mb_tmp = (-1.0/gdiff)*d01;
        real g_proj_norm = Norm(mb_tmp);
        
        if(g_proj_norm < lastgoodres) {
        	lastgoodres = g_proj_norm;
        	x_candidate = x;
        }
        
        residual=lastgoodres;
        if (residual < tolerance) {
            break;
        }
    }
    x=x_candidate;
    return current_iteration;
}
