%Es wird ein Schritt des Simplexverfahrens programmiert, ausgehend von

%  c'x=min, Ax=b, x>=0, A eine m x n Matrix mit m < n und Rang A = m

%bei bekannter Basisindexmenge I, d. h. B := A(I) besteht aus m linear 
%unabh\"angigen Spalten und Bx=b hat eine zul\"assige L\"osung. Das Ergebnis
%ist entweder eine neue Basisindexmenge J, oder das Problem hat wegen
%unbeschr\"ankter Zielfunktion keine L\"osung. Das Programm hat die Form
%
%      function [wert,x,J,info]=simplex(A,b,c,I)
%
%In info gibt es die folgenden Informationen:
 
%info=-1: Problem hat keine L\"osung, Rechnung beenden,
%info=0 : I definiert noch keine L\"osung, Rechnung fortsetzen,
%info=1 : I definiert eine L\"osung, Rechnung beenden.

%Ist die L\"osung erreicht (bei info=1), so enth\"alt x die L\"osung in 
%der urspr\"unglichen Numerierung.
%Programm: Gerhard Opfer, 12. Januar 1999
%===========================================================================
function [wert,x,J,info]=simplex(A,b,c,I);
%wert=0;x=[];J=I; %rein formal, keine Bedeutung
[m,n]=size(A);
if m>=n error('Matrix im Simplexprogramm hat falsche Dimension');end;
%===========================================================================
%                                            Erstes Gleichungssystem l\"osen
B=A(:,I);  y=B'\c(I);  svektor=A'*y;  notI=1:n;  notI(I)=[];
t=c-svektor; t=t(notI);
[mi,i]=min(t);                             %i kann mehrere Indizes enthalten
%===========================================================================
                                    %Pr\"ufung, ob L\"osung bereits erreicht
if mi >= 0, info=1; x=B\b; J=I; wert=c(I)'*x;         %L\"osung ist erreicht 
  xx=zeros(n,1); xx(I')=x; xx(notI')=[0;0];x=xx;
%===========================================================================
  else                                          %L\"osung ist nicht erreicht
  r=notI(i(1));  %Pivotspalte  %Zweites und drittes Gleichungssystem l\"osen
  a=A(:,r);      %r-te Spalte von A
  dx=B\[a,b];d=dx(:,1);x=dx(:,2);
%===========================================================================
                                     %Pr\"ufung, ob L\"osung nicht existiert
  if max(d)<=0,info=-1; x=[]; return, end;         %L\"osung existiert nicht
%===========================================================================
                                                               %weitermachen
  P=find(d>0);
  [delta,s]=min(x(P)./d(P));               %s kann mehrere Indizes enthalten
  zeile=I(P(s(1)));
  p=find(I==zeile); J=I; J(p)=r;               %Vertauschung ohne Sortierung
  info=0;  wert=c(I)'*x;
  xx=zeros(n,1); xx(I')=x; xx(notI')=zeros(n-m,1); x=xx;
end; %if
%===========================================================================
