function [x,lambda]=reduktion(C,d,A,b)
% ======================================================================
% Eine Implementation einer Reduktionsmethode fuer 
% quadratische, linear restringierte Optimierungsaufgaben der Form
%    Minimiere f(x) = 0.5 x^T C x + d^T x
%    unter der Nebenbedingung
%              A x - b = 0
%    
%    loest auch Probleme ohne Nebenbedingung (A=[], b=[])
%
% Erstellt fuer die Veranstaltung
% Weiterbildenes Studium fuer Ingenieure WS 98/99
%
% Helge Baumann, 13.9.1998, geaendert S.F. 
% ----------------------------------------------------------------------
% Aurufparameter:
% C : n x n Matrix
% A : m x n Matrix
% d : n Vektor
% b : n Vektor
% ---------------------------------------------------------------------
% Rueckgabeparameter:
% x : Loesung (Spaltenvektor!)
% lambda : Lagrange - Parameter
% ---------------------------------------------------------------------
% Programmvariablen:
% Q,R : Bilden QR- Zerlegung der Matrix A.
% Y,Z : Matrizen mit der Eigenschaft: A*Y = I, A*Z = 0, [Y Z] invertierbar.
% Cred: Matrix aus der reduzierten quadratischen Zielfunktion 
% D : Cholesky - zerlegte Matrix Cred zum Test auf positive Definitheit.
% u : Loesung des reduzierten LGS.
% ======================================================================

x=[]; lambda = [];

if nargin<4, error('fehlerhafte Eingabe!');end;

[m, n] = size(A);

% ======================================================================
% Beginn der Reduktionsmethode
% ======================================================================
%% --- Bestimmung der Matrizen Y und Z
[Q,R] = qr(A') ;    % QR-Zerlegung von A
R = R([1:m],:);     % Eleminiere die Nullzeilen aus R. R ist nun m x m Matrix
% Q1 = Q(:,[1:m]);
% Q2 = Q(:,[m+1:n]);
% A' = [Q1 Q2]*[R;0]=Q1*R		%Geiger 9.1.6
% Bestimmung von Y := Q1 (R^-1)^T durch Loesen 
% der n LGS von  R Y^T = Q1^T.
Y = (R \ Q(:,[1:m])')';    
% Z = Q2
Z = Q(:,[m+1:n]);   % Die letzten n-m Spalten werden zu Z zusammengefasst.

if isempty(A), 	    % damit fuer das unrestringierte Problem, die 
    Z=1;	    % die richtige Loesung berechnet wird. 
    Y=eye(size(C));
    b=zeros(length(C),1);
end;

%% --- Loese das LGS (Z^T C Z)u = -Z^T (d + C Y b)
Cred = Z' * C * Z;
[D,p] = chol(Cred);   % Cholesky - Zerlegung fuer Test auf pos. Definitheit
if (p > 0) 
   disp('Matrix  Z^T C Z  nicht positiv definit!');
   disp('Es konnte kein Minimum bestimmt werden');
   return;
end;

u = Cred \ (-Z' * (d + C*Y*b));

% Berechne x
x = Y*b + Z*u;

% Berechne Lagrangeparameter
lambda = Y' * (C*x + d);
if isempty(A), lambda=[]; end; % keine Nebenbedingung => kein lambda

return;












