function [x,l]=aktivemenge(C,d,A,b,x,mt)
% ======================================================================
% Eine Implementation einer Strategie der aktiven Menge fuer 
% quadratische, linear restringierte Optimierungsaufgaben der Form
%    Minimiere f(x) = 0.5 x^T C x + d^T x
%    unter der Nebenbedingung
%              a_i^T x >= b_i , i=1,...,mt 	(Ungleichungen)
%              a_i^T x  = b_i , i=mt+1,...,m    (Gleichungen)
% 
% Erstellt fuer die Veranstaltung
% Weiterbildenes Studium fuer Ingenieure WS 98/99
%
% S.F.
% ----------------------------------------------------------------------
% Aurufparameter:
% C : n x n Matrix
% A : m x n Matrix, mit den Zeilen a_i, i=1,...,m
% d : n Vektor
% b : n Vektor
% x : zulaessiger Punkt
% mt: Zahl der Ungleichungen
% ---------------------------------------------------------------------
% Rueckgabeparameter:
% x : Loesung (Spaltenvektor!)
% l : Lagrange - Parameter
% ---------------------------------------------------------------------
% Programmvariablen:
% ======================================================================

if nargin < 4, error('fehlerhafte Eingabe!'); end;
if nargin < 5, mt=0; end;

[m, n] = size(A);
% ------------------------------------------------------------------
% Ueberpruefung, ob zulaessiger Punkt
% ----------------------------------------------------------------
erg=A*x-b;
if ~isempty(find(erg(1:mt)<-10^(-10))), 
	error('kein zulaessiger Punkt, Ungleichung verletzt.');end;
if ~isempty(find(abs(erg(mt+1:m))>10^(-10))),  
	error('kein zulaessiger Punkt, Gleichung verletzt.');end;
% -------------------------------------------------------------------
abbruch=0;

I=find(abs(erg)<=10^(-10));	%	Aktive Indizes

% ======================================================================
% Beginn der Strategie der aktiven Menge
% ======================================================================
zahl=0;
while (abbruch==0)&(zahl<2*m)
  disp('--------------------------------------')
  zahl=zahl+1;
  disp(['Interationsschritt: ',num2str(zahl)]);
  disp(['     x = ',num2str(x')]);
  disp(['     I = ',num2str(I')]);

  % Loesung des reduzierten Problems
  if isempty(I), mth=0; Ah=[]; bh=[];	% keine Ungleichung aktiv
  else          
    mth=find(I==mt+1)-1;	  % I(:mth) Indizes der aktiven Ungleichungen
    if mt==m, mth=length(I);end;  % I(mth+1:) Indizes der Gleichungen
    Ah=A(I,:); 
    bh=zeros(size(I));		
  end;
  dh=C*x+d; 		
  [dx,l]=reduktion(C,dh,Ah,bh);
  disp(['    dx = ',num2str(dx')]);
  disp(['lambda = ',num2str(l')]);

  if max(abs(dx))<10^(-10)	%	dx=0
		   % Punkt ist Minimum bei den gewaehlten aktiven Ungleichungen
     if mth==0,    % nur Gleichungen aktiv, => Loesung gefunden. 
        abbruch=1; disp('Loesung gefunden');
        if n==2, xlabel('Loesung gefunden');  plot(x(1),x(2),'r*'); end;
     else
        [lq,iq]=min(l(1:mth));	% kleinster Index
        disp(['     q = ',num2str(I(iq))]);

        if (lq>=0); % keine Verbesserung durch weglassen eine Ungleichung mgl. 
	   abbruch=1; disp('Loesung gefunden');
           if n==2, xlabel('Loesung gefunden'); plot(x(1),x(2),'r*'); end;
        else	% Weglessen des Index iq kann Verbesserung bringen 
          J=find(I~=I(iq));
          I=I(J);
          disp(['     I = ',num2str(I)]);
	  
	  % Aufstellen und Loesen des reduzierten Systems
          if isempty(I), mth=0; Ah=[]; bh=[];	% keine Ungleichung aktiv
          else          
            mth=find(I==mt+1)-1;
            if mt==m, mth=length(I); end;
            Ah=A(I,:);
            bh=zeros(size(I));
          end;
          dh=C*x+d;  
          [dx,l]=reduktion(C,dh,Ah,bh);
        end;
     end;
  end;	%endif dx==0

% Bestimmung einer neuen Iterierten mit
% Schrittlaenge a=min(1,aquer),
% Berechnung von aquer, der Schrittlaenge bei der neue Ungleichung aktiv wird. 
% und evt. Bestimmung des Index p der aktiv werdenden Ungleichung
  if abbruch==0   
     s=dx;	%	Richtung
     disp(['     s = ',num2str(s')])

     % Bestimmung von aquer
     Iquer=setdiff([1:1:m],I); % nicht aktive Indizes
     Aquer=A(Iquer,:);
     erg=Aquer*s;	     
     Iquerh=find(erg<0);       
     Iquer=Iquer(Iquerh);      % Indizes die aktiv werden koennen	
     Aquer=A(Iquer,:);		
     bquer=b(Iquer);
     if ~isempty(Iquer);
         [aquer,p]=min((bquer-Aquer*x)./(Aquer*s));
         a=min(1,aquer);
         if a==aquer
              I=union(I,Iquer(p));
	      I=reshape(I,length(I),1);
         end;           
     else
         a=1;
     end;
     
     x=x+a*s;	% neue Iterierte
     disp([' aquer = ',num2str(aquer)])
     disp(['     I = ',num2str(I')]);
     disp(['     x = ', num2str(x')])
     if n==2, plot(x(1),x(2),'*');
        title(['Iterationsschritt:',num2str(zahl)])
        xlabel('weiter mit Mausklick ins Bild');
        ginput(1);	% Warten bis Maus geklickt. 
     end;
  end;	%	endif abbruch==0
end;	%	endwhile

return;












