function [x_f, y_f, z_f, s_f, fo, k] = pcQP_gen_2(G, g, A, b, C, d, x, y, z, s)
% Resuelve el problema de Programación Cuadrática
% min 0.5x’Gx + g’x
% s. a    Ax =  b         A=nA x n
%         Cx >= d         C=nC x n
% mediante el método predictor corrector de punto interior. Parte de (x,y,z,s)
%
 [nA,n]=size(A); [nC,~]=size(C); e=ones(nC,1);         % ¡OJO! n==ncc
 k=0 ; maxk=200; eps_L=1e-9; bc=1+max(norm(b),norm(d));
 if nargin<7, s=ones(nC,1); y=ones(nA,1); z=ones(nC,1); x=sqrt(n)*ones(n,1); end
 % Residuos
 rL  = G*x + g - A'*y - C'*z;
 rA  = -A*x + b;
 rC  = -C*x + s + d;
 rsz = s.*z;
 mu  = mean(rsz);               residuo=norm([rL;rA;rC;rsz])/bc;
 fprintf('\n              No fac.           Error         Brecha     Error');
 fprintf('\n              primal           Lagrang.        dual      relati.\n');
 fprintf('  Iter    A*x-b, C*x-s-d   G*x+g-A''*y-C''*z     z''*s      total\n');
 fprintf('  --------------------------------------------------------------\n');
 while k<=maxk && residuo>eps_L
  % Se resuelve el sistema para determinar dirección de Newton pura: predictor
    Mat = [G          -A'          -C'                 ; ...
          -A          zeros(nA,nA) zeros(nA,nC)        ; ...
          -C          zeros(nC,nA) sparse(-diag(s./z))];
    rhs = [-rL; -rA; -rC+rsz./z];
    [L,D,P] = ldl(Mat);
    dxyz_a = P*(L'\(D\(L\(P'*rhs))));
    dx_a   = dxyz_a(1:length(x));
    dy_a   = dxyz_a(length(x)+1:length(x)+length(y));
    dz_a   = dxyz_a(length(x)+length(y)+1:length(x)+length(y)+length(z));
    ds_a   = -((rsz+s.*dz_a)./z);
  % Cálculo de alpha_af
    alpha_a = 1;
    idx_z = find(dz_a<0);
    if isempty(idx_z)==0, alpha_a = min(alpha_a,min(-z(idx_z)./dz_a(idx_z))); end
    idx_s = find(ds_a<0);
    if isempty(idx_s)==0, alpha_a = min(alpha_a,min(-s(idx_s)./ds_a(idx_s))); end
    mu_a = ((z+alpha_a*dz_a)'*(s+alpha_a*ds_a))/nC;    % Affine duality gap
    sigma = (mu_a/mu)^3;                               % Centering parameter
  % Se resuelve dirección corrector
    rsz = rsz + ds_a.*dz_a - sigma*mu*e;
    rhs = [-rL; -rA; -rC+rsz./z];
    dxyz = P*(L'\(D\(L\(P'*rhs))));
    dx = dxyz(1:length(x));
    dy = dxyz(length(x)+1:length(x)+length(y));
    dz = dxyz(length(x)+length(y)+1:length(x)+length(y)+length(z));
    ds = -((rsz+s.*dz)./z);
  % Cálculo de alpha
    alpha = 1;
    idx_z = find(dz<0);
    if isempty(idx_z)==0, alpha = min(alpha,min(-z(idx_z)./dz(idx_z))); end
    idx_s = find(ds<0);
    if isempty(idx_s)==0, alpha = min(alpha,min(-s(idx_s)./ds(idx_s))); end
  % Obtención nuevo punto
    tau=max(.9995,1-mu)*alpha;        % Con factor eta=0.9995
    x = x+tau*dx; y = y+tau*dy; z = z+tau*dz; s = s+tau*ds; k = k+1;
  % Recalcular residuos en nuevo punto, brecha dual e imprimir iteración
    rL  = G*x + g - A'*y - C'*z;
    rA  = -A*x + b;
    rC  = -C*x + s + d;             nrC=norm([rA;rC]);
    rsz = s.*z;
    mu  = mean(rsz);                residuo=norm([rL;rA;rC;rsz])/bc;
    fprintf('%5i %15.2e %17.2e %13.2e %10.2e\n',k,nrC,norm(rL),mu,residuo);
 end
 x_f = x; y_f = y; z_f = z; s_f = s; fo = 0.5*x'*G*x+g'*x;
end