function [x_out, iter] = qp_m_rwm4(alg_data, bt, max_iter)
% global wsxn
% global wsyn
% global wsvn
% global wsqn
% global wstmp_var_n2n
% global wstmp_var_pn
% global wstmp_var_p2n
% global wstmp_var_rn
% global wstmp_var_r2n
% global wsxx
% global wsyx
% global wsvx
% global wsqx
% global wstmp_var_n2x
% global wstmp_var_px
% global wstmp_var_p2x
% global wstmp_var_rx
% global wstmp_var_r2x
% % wsxn = 1e8*ones(81,1);
% % wsyn = 1e8*ones(81,1);
% % wsvn = 1e8*ones(81,1);
% % wsqn = 1e8*ones(81,1);
% % wstmp_var_n2n = 1e8*ones(81,1);
% % wstmp_var_pn = 1e8*ones(81,1);
% % wstmp_var_p2n = 1e8*ones(81,1);
% % wstmp_var_rn = 1e8*ones(81,1);
% % wstmp_var_r2n = 1e8*ones(81,1);
% % wsxx = -1e8*ones(81,1);
% % wsyx = -1e8*ones(81,1);
% % wsvx = -1e8*ones(81,1);
% % wsqx = -1e8*ones(81,1);
% % wstmp_var_n2x = -1e8*ones(81,1);
% % wstmp_var_px = -1e8*ones(81,1);
% % wstmp_var_p2x = -1e8*ones(81,1);
% % wstmp_var_rx = -1e8*ones(81,1);
% % wstmp_var_r2x = -1e8*ones(81,1);


ws.tmp_var_n_orig = (zeros(alg_data.vdim,1));
ws.tmp_var_n2_orig = (zeros(alg_data.vdim,1));
ws.v = (zeros(alg_data.vdim,1));
ws.y = (zeros(alg_data.vdim,1));
ws.y1 = (zeros(alg_data.vdim,1));
ws.q1 = (zeros(alg_data.vdim,1));
ws.q2 = (zeros(alg_data.vdim,1));
ws.q = (zeros(alg_data.vdim,1));
ws.r = (zeros(alg_data.vdim,1));
ws.tmp_var_n = (zeros(alg_data.vdim,1));
ws.tmp_var_n2 = (zeros(alg_data.vdim,1));
ws.x = (zeros(alg_data.vdim,1));
ws.y_old = (zeros(alg_data.vdim,1));
ws.tmp_var_p = (zeros(alg_data.vdim,1));
ws.tmp_var_p2 = (zeros(alg_data.vdim,1));
ws.tmp_var_r = (zeros(alg_data.vdim,1));
ws.tmp_var_r2 = (zeros(alg_data.vdim,1));
ws.v_old = (zeros(alg_data.vdim,1));

% theta = (1);
% theta_old = (1);
jj = (0);
ix = (0);
test = (0);
test0 = (0);
test1 = (0);

ws.q1 = alg_data.Gg;
ws.q2 = alg_data.Gb * bt;
ws.q = ws.q1 + ws.q2;
ws.r = ws.tmp_var_n_orig + ws.tmp_var_n2_orig;

while ((jj < max_iter))%0) && (cond < 0))
  jj = jj+1;
  
  ws.tmp_var_n = alg_data.H * ws.v;
  ws.tmp_var_n2 = ws.tmp_var_n + ws.q;
  ws.x = ws.v - ws.tmp_var_n2;
  ws.y_old = ws.y;
  % clip(ws->y,ws->l,ws->u,81);      - inlined
%  ws.y1 = max(ws.x, alg_data.h.Lb);
  for ix = 1:alg_data.vdim 
    if alg_data.h.Lb(ix) > ws.x(ix) 
      ws.y1(ix) = alg_data.h.Lb(ix); 
    else 
      ws.y1(ix) = ws.x(ix);
    end%if
  end%for ix
%  ws.y = min(ws.y1, alg_data.h.Ub);
  for ix = 1:alg_data.vdim 
    if alg_data.h.Ub(ix) < ws.y1(ix) 
      ws.y(ix) = alg_data.h.Ub(ix); 
    else 
      ws.y(ix) = ws.y1(ix);
    end%if
  end%for ix
  % clip end
  ws.tmp_var_p = ws.y - ws.y_old;
  
  % theta_old = theta;
  % theta = (1+sqrt(1+4*theta_old^2))/2;
  % ws.tmp_var_p2 = (theta_old-1)/theta * ws.tmp_var_p;
  ws.tmp_var_p2 = alg_data.beta_vec(jj) * ws.tmp_var_p;
  ws.v_old = ws.v;
  ws.v = ws.tmp_var_p2 + ws.y;
  % restart(ws->y,ws->y_old,ws->v,ws->v_old,ws->tmp_var_r,ws->tmp_var_r2,81);    - inlined
  ws.tmp_var_r = ws.v_old - ws.y;
  ws.tmp_var_r2 = ws.y - ws.y_old;
%  test = ws.tmp_var_r' * ws.tmp_var_r2;
  test = 0; 
  for ix = 1:alg_data.vdim
    test0 = ws.tmp_var_r(ix) * ws.tmp_var_r2(ix);
    test1 = test0 * alg_data.Einv2(ix,ix);
    test = test + test1;
  end%for ix
  if test > 0
    ws.v = ws.y_old;
    ws.y = ws.y_old;
    % theta = 1; % reset theta?       % would need separate beta_vec index!
  end
  % restart end
end

ws.tmp_var_n_orig = alg_data.R * ws.y;
x_out = ws.tmp_var_n_orig + ws.r;
iter = jj;

% for ix = 1:alg_data.vdim
%   if wsxn(ix)>ws.x(ix), wsxn(ix)=ws.x(ix); end
%   if wsyn(ix)>ws.y(ix), wsyn(ix)=ws.y(ix); end
%   if wsvn(ix)>ws.v(ix), wsvn(ix)=ws.v(ix); end
%   if wsqn(ix)>ws.q(ix), wsqn(ix)=ws.q(ix); end
%   if wstmp_var_n2n(ix)>ws.tmp_var_n2(ix), wstmp_var_n2n(ix)=ws.tmp_var_n2(ix); end
%   if wstmp_var_pn(ix)>ws.tmp_var_p(ix), wstmp_var_pn(ix)=ws.tmp_var_p(ix); end
%   if wstmp_var_p2n(ix)>ws.tmp_var_p2(ix), wstmp_var_p2n(ix)=ws.tmp_var_p2(ix); end
%   if wstmp_var_rn(ix)>ws.tmp_var_r(ix), wstmp_var_rn(ix)=ws.tmp_var_r(ix); end
%   if wstmp_var_r2n(ix)>ws.tmp_var_r2(ix), wstmp_var_r2n(ix)=ws.tmp_var_r2(ix); end
%   if wsxx(ix)<ws.x(ix), wsxx(ix)=ws.x(ix); end
%   if wsyx(ix)<ws.y(ix), wsyx(ix)=ws.y(ix); end
%   if wsvx(ix)<ws.v(ix), wsvx(ix)=ws.v(ix); end
%   if wsqx(ix)<ws.q(ix), wsqx(ix)=ws.q(ix); end
%   if wstmp_var_n2x(ix)<ws.tmp_var_n2(ix), wstmp_var_n2x(ix)=ws.tmp_var_n2(ix); end
%   if wstmp_var_px(ix)<ws.tmp_var_p(ix), wstmp_var_px(ix)=ws.tmp_var_p(ix); end
%   if wstmp_var_p2x(ix)<ws.tmp_var_p2(ix), wstmp_var_p2x(ix)=ws.tmp_var_p2(ix); end
%   if wstmp_var_rx(ix)<ws.tmp_var_r(ix), wstmp_var_rx(ix)=ws.tmp_var_r(ix); end
%   if wstmp_var_r2x(ix)<ws.tmp_var_r2(ix), wstmp_var_r2x(ix)=ws.tmp_var_r2(ix); end
% end%for


