function [xk,yk,lmdk,fs,ck,k,L] = NLP_IP_nc(fcname,gname,hname,x0,y0,lmd0,tau0)
% Método de punto interior para min. f(x)
%                               s. a g(x) >= 0
%
 xk = x0(:); yk = y0(:); q = length(y0); lmdk = lmd0(:); 
 tau = tau0; n = length(x0); E2 = eye(n); ei = 1e-5; eo = 1e-6;
 q1 = q + 1; L = 0; K = 0; do = 1;
 % Bucle exterior
 while do >= eo,
    xk0 = xk; yk0 = yk; lmdk0 = lmdk;
    k = 0; di = 1;
  % Bucle interior
    while di >= ei,
       bt = 0;
       fck = feval(fcname,xk); Gk = feval(gname,xk);
       fk  = fck(1);    ck  = fck(2:q1);
       gk  = Gk(:,1);   Ak  = Gk(:,2:q1)'; 
       yki = 1./yk;     wwk = yki.*lmdk;
       Yki = diag(yki); Lk  = diag(lmdk);
       xlk = [n;q;xk;lmdk];
       Hk  = feval(hname,xlk); 
       ew  = min(eig(Hk));
       if ew <= 0
          eb = -1.2*ew + 1e-6;
          Hk = Hk + eb*E2;
       end
       Nk  = Hk + Ak'*(Lk*Yki)*Ak;
       R   = cholinc(sparse(Nk),'inf');
       Nki = R\(R'\eye(n));                            % Mat = [-Hk         zeros(n,q)  Ak'      ;...      
       gak = tau*yki - lmdk;                           %        zeros(q,n) -diag(wwk)  -eye(q)   ;...
       rk  = yk - ck;                                  %        Ak         -eye(q)      zeros(q)];
       dx = Nki*(-gk + tau*Ak'*yki + Ak'*(wwk.*rk));   % [L D P] = ldl(Mat);
       dL = wwk.*(rk - Ak*dx) + gak;                   % rhs = [gk-Ak'*lmdk; -gak; rk];
       dy = Ak*dx - rk;                                % [dx; dy; dL] = P*(L'\(D\(L\(P'*rhs))));
     % Comprobar si la dirección obtenida es de descenso en función de mérito
       ksi = gk - tau*Ak'*yki; 
       sk = -ksi'*(Nki*ksi) + tau*(yki'*rk) + ksi'*(Nki*Ak')*(wwk.*rk)...
            -bt*norm(rk)^2;
       if sk >= 0,
          bmin = (sk+bt*norm(rk)^2)/(norm(rk)^2);      % bmin = sk/(norm(rk)^2);
          bt   = 10*bmin; 
       end
     % Cálculo de line searchen x e y
       awk  = max([-dy./yk; -dL./lmdk]);
       ak   = 0.95/awk;
       aks  = lsearch_merit(fcname,xk,yk,dx,dy,ak,tau,bt);
       xk   = xk + aks*dx;
       yk   = yk + aks*dy;
       lmdk = lmdk + aks*dL;
       k  = k + 1;
       di = aks*(norm(dx) + norm(dy) + norm(dL));
    end
  % Recalcular parámetro de barrera tau.
    tsi = q*min(yk.*lmdk)/(yk'*lmdk);
    tau = 0.1*(min([0.05*(1-tsi)/tsi,2]))^3*(yk'*lmdk)/q;
    K  = k + K; L = L + 1;
    do = norm(xk-xk0) + norm(yk-yk0) + norm(lmdk-lmdk0);
 end
 fck = feval(fcname,xk); fs  = fck(1); k = K;
end

function [a psm] = lsearch_merit(fcname,xk,yk,dx,dy,ak,tau,bt)
 aa = (0:0.01:1)*ak;
 ps = zeros(101,1);
 for i = 1:101
    ai = aa(i);
    xki = xk + ai*dx;
    yki = yk + ai*dy;
    fci = feval(fcname,xki);
    fi = fci(1);
    cki = fci(2:end);
    yci = yki - cki;
    ps(i) = fi - tau*sum(log(yki)) + 0.5*bt*(yci'*yci);  % Fun. mérito
 end
 [psm,ind] = min(ps); a = min(1,aa(ind));
end



function z = f_ex15_6(x)
 x1 = x(1); x2 = x(2);
 ck = [-x1^2+x2; -x1+x2^2];
 fk = (x1-2)^2+(x2-1)^2;
 z = [fk; ck];
end

function z = g_ex15_6(x)
 x1 = x(1); x2 = x(2);
 Ak = [-2*x1, 1 ; -1, 2*x2];
 gk = [2*(x1-2); 2*(x2-1)];
 z = [gk,  Ak'];
end

function z = h_ex15_6(xlk)
 lmdk = xlk(5:6);
 lmd1 = lmdk(1); lmd2 = lmdk(2);
 z = [2+2*lmd1,  0; 0,  2-2*lmd2];
end

function z = f_ex15_7(x)
 x1 = x(1); x2 = x(2);
 ck = [-x1-x2+2; -x1^2+x2];
 fk = (x1-2)^2+(x2-1)^2;
 z = [fk; ck];
end

function z = g_ex15_7(x)
 x1 = x(1); x2 = x(2);
 Ak = [-1, -1 ; -2*x1, 1];
 gk = [2*(x1-2); 2*(x2-1)];
 z = [gk,  Ak'];
end

function z = h_ex15_7(xlk)
 lmdk = xlk(5:6);
 lmd1 = lmdk(1); lmd2 = lmdk(2);
 z = [2+2*lmd2, 0; 0, 2];
end

function z = f_ex15_8(x)
 x1 = x(1); x2 = x(2);
 ck = [-0.25*x1^2-x2^2+1; x1-2*x2+1; -x1+2*x2-1];
 fk = (x1-2)^2+(x2-1)^2;
 z = [fk; ck];
end


function z = g_ex15_8(x)
 x1 = x(1); x2 = x(2);
 Ak = [-0.5*x1, -2*x2; 1, -2; -1, 2];
 gk = [2*(x1-2); 2*(x2-1)];
 z = [gk,  Ak'];
end

function z = h_ex15_8(xlk)
 lmdk = xlk(5:7);
 lmd1 = lmdk(1); lmd2 = lmdk(2); lmd3 = lmdk(3);
 z = [2+0.5*lmd1, 0; 0, 2+2*lmd1];
end