function [rozw,zopt]=sympleks(A,b,c,krok)
%wywolanie funkcji:
%[rozw]=sympleks(A,b,c)
% wejscie:
% A - macierz ukladu ograniczen,
% b - wektor kolumnowy prawych stron ograniczen, 
% c - wektor kolumnowy wspolczynnikow funkcji celu,
% A, b i c powiiny tworzyc postac standardowa zadania
% krok- jest opcjonalnym parametrem 
% krok=1 (default) - praca krokowa, w przeciwnym razie podanie wyniku koncowego 
%wyjscie:
% rozw - rozwiazanie optymalne 
% zopt - optymalna wartosc funkcji celu

%inicjalizacja zmiennych i rezerwacja pamieci
epsi=1e-12;%numeryczna dokladnosc rozwiazania
if nargin<4, krok=1;end
[M,N]=size(A);
xb=zeros(M,1);
kanon=1;
%sprawdzenie kanonicznosci ukladu rownan
I=eye(M);
Ar=[A,b].*repmat(2*((b>=0)-0.5),1,N+1);
for i=1:M
    Ii=I(:,i);
    for j=1:N
        if Ar(i,j)>0
           if all(Ar(:,j)./Ar(i,j)==Ii)
               xb(i)=j; Ar(i,:)=Ar(i,:)./Ar(i,j); break
           end
        end
    end
end
if any(xb==0), kanon=0; end
Y = [Ar;-c',0];

%FAZA I
if ~kanon
   if krok 
    disp('Zadanie nie jest w postaci kanonicznej!!!')
    disp('tworzymy zadanie pomocnicze poprzez dodanie zmiennych sztucznych')
   end
   Ar=[Ar(:,1:N),I(:,xb==0),Ar(:,N+1)];
   dl=length(xb(xb==0)); ind=N+1:N+dl; %dlugosc i indeksy wektora zmiennych sztucznych 
   xb(xb==0)=ind; 
   Y = [Ar;-c',zeros(1,dl),0;zeros(1,N),-ones(1,dl),0];
   if krok
    disp('------------------------')
    disp('POCZATEK FAZY I:')
    disp('tabliczka sympleksowa (w kolejnych wierszach: ograniczenia, oryg. funkcja celu Z, pomocnicza funkcja celu W)');
    disp('Y=')
    disp(Y)
    disp('normalizacja pomocniczej funkcji celu W');
   end
   Y(end,:)=Y(end,:)+sum(Y(xb>N,:),1);
   if krok, disp('Y='), disp(Y), disp('baza'), disp(xb'), end
      
   iter=1;
   while any(Y(end,1:end-1)>epsi)
     if krok,  
      disp('------------------------')
      disp(['I faza, iteracja ',num2str(iter)])
      disp('------------------------')
     end 
     [maxk,kw]=max(Y(end,1:end-1)); %wybor kolumny wiodacej
     indd=Y(1:M,kw)>eps;  %tablica wspolczynnikow przyrostowych 
     if ~any(indd), error('rozw. I fazy jest nieograniczone.'), end
     przyrosty=Inf*ones(M,1); przyrosty(indd)=Y(indd,end)./Y(indd,kw);
     [minw,ww]=min(przyrosty);
     Y(ww,:)=Y(ww,:)./Y(ww,kw);
     pom=Y(:,kw); pom(ww)=0; pom=repmat(Y(ww,:),M+2,1).*repmat(pom,1,N+dl+1);
     Y=Y-pom; xb(ww)=kw;
     if krok 
      disp('Y='), disp(Y), disp('baza'), disp(xb')
      pause
     end
     iter=iter+1;
   end  
   if krok, disp('KONIEC FAZY I:');end
   if Y(end,end)>epsi, error('W>0 --> Pusty zbior rozw. dopuszczalnych'), end
   indb=ones(1,N+dl+1);indb(ind)=0;indb=logical(indb);
%   indb=zeros(1,N+dl+1);indb(ind)=1,indb=logical(indb)
   Y=Y(1:M+1,indb);
   if krok 
      disp('Y='), disp(Y), disp('baza'), disp(xb')
      pause
    disp('-----------------------------')
    disp('-----------------------------')
    disp('POCZATEK FAZY II:');
   end
end

%FAZA II
iter=1;
while any(Y(end,1:end-1)>0)
   if krok 
    disp('------------------------')
    disp(['II faza, iteracja ',num2str(iter)])
    disp('------------------------')
   end
   [maxk,kw]=max(Y(end,1:end-1)); %wybor kolumny wiodacej
   indd=Y(1:M,kw)>0;  %tablica wspolczynnikow przyrostowych 
   if ~any(indd), error('Rozwiazanie jest nieograniczone'), end
   przyrosty=Inf*ones(M,1); przyrosty(indd)=Y(indd,end)./Y(indd,kw);
   [minw,ww]=min(przyrosty);
   Y(ww,:)=Y(ww,:)./Y(ww,kw);
   pom=Y(:,kw); pom(ww)=0; pom=repmat(Y(ww,:),M+1,1).*repmat(pom,1,N+1);
   Y=Y-pom; xb(ww)=kw;
   if krok 
      disp('Y='), disp(Y), disp('baza'), disp(xb')
      pause
   end
   iter=iter+1;
end  
rozw=zeros(1,N);rozw(xb)=Y(1:M,end);zopt=Y(end,end);
if krok, disp('rozwiazanie optymalne'), disp('x='), disp(rozw), end