function [x_out, iter] = qp_h_rwm4(bt, max_iter)

alg_data = coder.load('alg_data4.mat');

unroll_pe_o = 9;
unroll_p_im = 2; 
unroll_e_im = 81; 
unroll_ml_o = 9; 
unroll_ml_im = 81;

% theta = (1);
% theta_old = (1);
%jj = (0);
ix = (0);
test = (0);
test0 = (0);
test1 = (0);
btdim = int16(numel(bt));
vdim = int16(numel(alg_data.Gg));

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

last64 = uint16(vdim/2);
imax64 = uint16(bitshift(vdim,-1));
last32 = uint16(last64/2);
imax32 = uint16(bitshift(last64,-1));
last16 = uint16(last32/2);
imax16 = uint16(bitshift(last32,-1));
last8 = uint16(last16/2);
imax8 = uint16(bitshift(last16,-1));
last4 = uint16(last8/2);
imax4 = uint16(bitshift(last8,-1));
last2 = uint16(last4/2);
imax2 = uint16(bitshift(last4,-1));
last1 = uint16(last2/2);
imax1 = uint16(bitshift(last2,-1));

tmpprodmat = coder.nullcopy(zeros(vdim,vdim));
vt64 = coder.nullcopy(zeros(vdim,last64));
vt32 = coder.nullcopy(zeros(vdim,last32));
vt16 = coder.nullcopy(zeros(vdim,last16));
vt8 = coder.nullcopy(zeros(vdim,last8));
vt4 = coder.nullcopy(zeros(vdim,last4));
vt2 = coder.nullcopy(zeros(vdim,last2));
vt1 = coder.nullcopy(zeros(vdim,last1));

tmpprodmatm = coder.nullcopy(zeros(vdim,vdim));
vt64m = coder.nullcopy(zeros(vdim,last64));
vt32m = coder.nullcopy(zeros(vdim,last32));
vt16m = coder.nullcopy(zeros(vdim,last16));
vt8m = coder.nullcopy(zeros(vdim,last8));
vt4m = coder.nullcopy(zeros(vdim,last4));
vt2m = coder.nullcopy(zeros(vdim,last2));
vt1m = coder.nullcopy(zeros(vdim,last1));
%% prologue

ws.q1 = alg_data.Gg;
%coder.hdl.loopspec('stream')
ws.q2 = alg_data.Gb * bt;
ws.q = ws.q1 + ws.q2;
ws.r = ws.tmp_var_n_orig + ws.tmp_var_n2_orig;


%% main iteration loop

% while ((jj < max_iter))%0) && (cond < 0))
%   jj = jj+1;
%coder.hdl.loopspec('stream')
%for jj = 1:max_iter
for jj = 1:4 %%%%%%%%%%%%%%%%%%%%%%%%%%%20     
  
  %ws.tmp_var_n = alg_data.H * ws.v;
  
  % firstly compute the products - the whole matrix
  %coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
  coder.hdl.loopspec('unroll',unroll_ml_o)           % full unroll 81
  for ir = 1:vdim
    %coder.hdl.loopspec('unroll',vdim)           % full unroll 81
    coder.hdl.loopspec('unroll',unroll_ml_im)           % full unroll 81
    for ic = 1:vdim
      tmpprodmatm(ir,ic) = alg_data.H(ir,ic)*ws.v(ic);
    end%for ic
  end%for ir
  
  % ...then tree-sum the products   !!!hard-coded for up to 128 elements
  vt64m(:,last64) = tmpprodmatm(:,last64);
  %coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
  coder.hdl.loopspec('unroll',unroll_ml_o)           % full unroll 81
  for ir = 1:vdim
    %coder.hdl.loopspec('unroll', imax64)
    coder.hdl.loopspec('unroll', 40)
    for i = 1:imax64
      j = uint16(last64 + i);
      vt64m(ir,i) = tmpprodmatm(ir,i) + tmpprodmatm(ir,j);
    end
  end%for ir
  
  vt32m(:,last32) = vt64m(:,last32);
  %coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
  coder.hdl.loopspec('unroll',unroll_ml_o)           % full unroll 81
  for ir = 1:vdim
    %coder.hdl.loopspec('unroll', imax32)
    coder.hdl.loopspec('unroll', 20)
    for i = 1:imax32
      j = uint16(last32 + i);
      vt32m(ir,i) = vt64m(ir,i) + vt64m(ir,j);
    end
  end%for ir
  
  vt16m(:,last16) = vt32m(:,last16);
  %coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
  coder.hdl.loopspec('unroll',unroll_ml_o)           % full unroll 81
  for ir = 1:vdim
    %coder.hdl.loopspec('unroll', imax16)
    coder.hdl.loopspec('unroll', 10)
    for i = 1:imax16
      j = uint16(last16 + i);
      vt16m(ir,i) = vt32m(ir,i) + vt32m(ir,j);
    end
  end%for ir
  
  vt8m(:,last8) = vt16m(:,last8);
  %coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
  coder.hdl.loopspec('unroll',unroll_ml_o)           % full unroll 81
  for ir = 1:vdim
    %coder.hdl.loopspec('unroll', imax8)
    coder.hdl.loopspec('unroll', 5)
    for i = 1:imax8
      j = uint16(last8 + i);
      vt8m(ir,i) = vt16m(ir,i) + vt16m(ir,j);
    end
  end%for ir
  
  vt4m(:,last4) = vt8m(:,last4);
  %coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
  coder.hdl.loopspec('unroll',unroll_ml_o)           % full unroll 81
  for ir = 1:vdim
    %coder.hdl.loopspec('unroll', imax4)
    coder.hdl.loopspec('unroll', 3)
    for i = 1:imax4
      j = uint16(last4 + i);
      vt4m(ir,i) = vt8m(ir,i) + vt8m(ir,j);
    end
  end%for ir
  
  vt2m(:,last2) = vt4m(:,last2);
  %coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
  coder.hdl.loopspec('unroll',unroll_ml_o)           % full unroll 81
  for ir = 1:vdim
    %coder.hdl.loopspec('unroll', imax2)    %! with vdim=81, imax2=1, no loop
    for i = 1:imax2
      j = uint16(last2 + i);
      vt2m(ir,i) = vt4m(ir,i) + vt4m(ir,j);
    end
  end%for ir
  
  vt1m(:,last1) = vt2m(:,last1);
  %coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
  coder.hdl.loopspec('unroll',unroll_ml_o)           % full unroll 81
  for ir = 1:vdim
    %for i = 1:imax1
    i = imax1;                        % imax1 should always be 1, no loop
    j = uint16(last1 + i);
    vt1m(ir,i) = vt2m(ir,i) + vt2m(ir,j);
    %end
  end%for ir
  
  ws.tmp_var_n = vt1m;
  %end of 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:vdim 
  %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:vdim 
  %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;
  %coder.hdl.loopspec('stream')
  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; 
  %coder.hdl.loopspec('stream')
  for ix = 1:vdim
  %for ix = 1:alg_data.vdim
    %coder.hdl.loopspec('stream')
    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

%%epilogue

%ws.tmp_var_n_orig = alg_data.R * ws.y;      % treesum version, with expanded indermediate variables, separate strictly simple nested loops

% firstly compute the products - the whole matrix 
%coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
coder.hdl.loopspec('unroll',unroll_pe_o)           % full unroll 81
for ir = 1:vdim
%for ir = 1:alg_data.vdim
  %coder.hdl.loopspec('unroll',vdim)           % full unroll 81
  coder.hdl.loopspec('unroll',unroll_e_im)           % full unroll 81
  for ic = 1:vdim
  %for ic = 1:alg_data.vdim
    tmpprodmat(ir,ic) = alg_data.R(ir,ic)*ws.y(ic); 
  end%for ic
end%for ir
  
% ...then tree-sum the products   !!!hard-coded for up to 128 elements
vt64(:,last64) = tmpprodmat(:,last64);
%coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
coder.hdl.loopspec('unroll',unroll_pe_o)           % full unroll 81
for ir = 1:vdim
  %coder.hdl.loopspec('unroll', imax64)
  coder.hdl.loopspec('unroll', 40)
  for i = 1:imax64
    j = uint16(last64 + i);
    vt64(ir,i) = tmpprodmat(ir,i) + tmpprodmat(ir,j);
  end
end%for ir

vt32(:,last32) = vt64(:,last32);
%coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
coder.hdl.loopspec('unroll',unroll_pe_o)           % full unroll 81
for ir = 1:vdim
  %coder.hdl.loopspec('unroll', imax32)
  coder.hdl.loopspec('unroll', 20)
  for i = 1:imax32
    j = uint16(last32 + i);
    vt32(ir,i) = vt64(ir,i) + vt64(ir,j);
  end
end%for ir

vt16(:,last16) = vt32(:,last16);
%coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
coder.hdl.loopspec('unroll',unroll_pe_o)           % full unroll 81
for ir = 1:vdim
  %coder.hdl.loopspec('unroll', imax16)
  coder.hdl.loopspec('unroll', 10)
  for i = 1:imax16
    j = uint16(last16 + i);
    vt16(ir,i) = vt32(ir,i) + vt32(ir,j);
  end
end%for ir

vt8(:,last8) = vt16(:,last8);
%coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
coder.hdl.loopspec('unroll',unroll_pe_o)           % full unroll 81
for ir = 1:vdim
  %coder.hdl.loopspec('unroll', imax8)
  coder.hdl.loopspec('unroll', 5)
  for i = 1:imax8
    j = uint16(last8 + i);
    vt8(ir,i) = vt16(ir,i) + vt16(ir,j);
  end
end%for ir

vt4(:,last4) = vt8(:,last4);
%coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
coder.hdl.loopspec('unroll',unroll_pe_o)           % full unroll 81
for ir = 1:vdim
  %coder.hdl.loopspec('unroll', imax4)
  coder.hdl.loopspec('unroll', 3)
  for i = 1:imax4
    j = uint16(last4 + i);
    vt4(ir,i) = vt8(ir,i) + vt8(ir,j);
  end
end%for ir

vt2(:,last2) = vt4(:,last2);
%coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
coder.hdl.loopspec('unroll',unroll_pe_o)           % full unroll 81
for ir = 1:vdim
  %coder.hdl.loopspec('unroll', imax2)    %! with vdim=81, imax2=1, no loop
  for i = 1:imax2
    j = uint16(last2 + i);
    vt2(ir,i) = vt4(ir,i) + vt4(ir,j);
  end
end%for ir

vt1(:,last1) = vt2(:,last1);
%coder.hdl.loopspec('stream', vdim)           % full stream 81 (can't stream outer loop less)
coder.hdl.loopspec('unroll',unroll_pe_o)           % full unroll 81
for ir = 1:vdim
  %for i = 1:imax1
  i = imax1;                        % imax1 should always be 1, no loop  
  j = uint16(last1 + i);
  vt1(ir,i) = vt2(ir,i) + vt2(ir,j);
  %end
end%for ir

ws.tmp_var_n_orig = vt1; 
%end of ws.tmp_var_n_orig = alg_data.R * ws.y;

x_out = ws.tmp_var_n_orig + ws.r;
iter = max_iter;                      % with fixed iterations, redundant output 






