     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);

        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);
    }
    return current_iteration;
}

}

uint ChSolverAPGD::SolveAPGD(custom_vector<real> &x, const custom_vector<real> &b, const uint max_iter) {
    real gdiff = 0.000001;
    real lastgoodres=10e30;
    real theta_k=1.0;
    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);
