         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);
-        //mb_tmp = (-1.0/gdiff)*(mb_tmp - x);
-        //real g_proj_norm = Norm(mb_tmp);
+        
+        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;
 }

}

uint ChSolverAPGD::SolveAPGD(custom_vector<real> &x, const custom_vector<real> &b, const uint max_iter) {
    //real gdiff = 0.000001;
    real theta_k=1.0;
    custom_vector<real> ms, mg_tmp2, mb_tmp, x01(x.size()), d01;
    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;
    custom_vector<real>  x_candidate = 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);
