 }
 
 uint ChSolverAPGD::SolveAPGD(custom_vector<real> &x, const custom_vector<real> &b, const uint max_iter) {
-   /* real gdiff = 0.000001;
-    custom_vector<real> ms, mg_tmp2, mb_tmp;
+    //real gdiff = 0.000001;
+    real theta_k=1.0;
+    custom_vector<real> ms, mg_tmp2, mb_tmp, x01, d01;
     custom_vector<real> mg = ShurProduct(x) - b;
-    real L_k = Norm(ShurProduct(x - 1)) / Norm(x - 1);
+    Thrust_Fill(x01, 1.0);
+    d01=x-x01;
+    real L_k = Norm(ShurProduct(d01)) / Norm(d01);
     real t_k = 1 / L_k;
     custom_vector<real>  my = x;
     custom_vector<real>  mx = x;

        }
        real theta_k1 = (-pow(theta_k, 2) + theta_k * sqrt(pow(theta_k, 2) + 4)) / 2.0;
        real  beta_k1 = theta_k * (1.0 - theta_k) / (pow(theta_k, 2) + theta_k1);
        my = mx + beta_k1 * (mx - x);
        ms = mx - x;
        if (Dot(mg, ms) > 0) {
            my = mx;
            theta_k1 = 1.0;

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