%
%             Actividad practica del 8 de febrero, 2012  
%
%    Programacion cuadratica. Metodos basados en el espacio nulo y en el
%    complemento de Schur. JL Morales Departamento de Matematicas. I T A M
%
%     PROBLEMA:
%                 minimizar  1/2 p'*W*p + g'*p                       (1)
%                sujeta a   A*p + c =0
%   
%     en donde  W es simetrica y A es de rango completo
%
%--------------------------------------------------------------------------
%
randn('state',16);     n = 10; m = 5;
%
W = randn(n,n);        A = randn(m,n);   
W = W'*W;              % .... W es simetrica
g = randn(n,1);        c = randn(m,1);
%
% ... calculo del paso normal y del paso tangencial via una factorizacion
%     QR (completa) de A'.
%
[Q,T] = qr(A');
Q1    = Q(:,1:m);      %  base para calcular el paso normal
Q2    = Q(:,m+1:n);    %  base del espacio nulo para calcular el 
                       %  paso tangencial
R     = T(1:m,1:m);
p_N   = -Q1*( R'\c );  %  paso normal 
P     = Q2'*W*Q2;            
me    = min(eig(P));   %  verificar que la proyeccion de W en el espacio nulo de
                       %  A es positiva definida, y que el rango de A es m
rA = rank(A);                              
fprintf(' menor valor propio de W proyectada ........   %8.2e \n', me);
fprintf(' rango de A ................................ %3i   \n', rA);
if  rA < m || me < 0  
    fprintf(' la matriz A es de rango incompleto o la Hessiana proyectada no es positiva definida');
    break
end
%
% ... sustituyendo p = p_N + Q2*d_z en la funcion objetivo (1) tenemos el 
%     PROBLEMA
%               minimizar  1/2 d_z'*Q2'*W*Q2*d_z + ( Wp_N + g )'*Q2*d_z
%   
d_z = -P\( Q2'*( W*p_N + g ) );
p1 = p_N + Q2*d_z;
%
%      Solucion via un metodo basado en el complemento de Schur
%
mW = min(eig(W));      % verificar que W es spd
fprintf(' menor valor propio de W ........   %8.2e \n', mW);
if mW <= 0 
    fprintf(' la Hessiana no es positiva definida');
    break
end
%
s = (A*inv(W)*A')\( -c + A*inv(W)*g );
p2 = inv(W)*( -g + A'*s );
%
%     Solucion via la matriz de KKT
%
K = [  W       A'      ;
       A   zeros(m,m) ];   b = [ -g ; -c ];
%   
z = K\b;
p3 = z(1:n);  pnorm = norm(p1);
fprintf( ' errores %8.2e   %8.2e  \n', norm(p1-p2)/pnorm, norm(p1-p3)/pnorm );