{$M 65000,0,0}
{$N+ $E}
program Simplex; { Este programa recebe um arquivo texto contendo problemas
                   de programacao linear, no formato: Max cx sujeito a Ax=b
                   e x>=0, e devolve um arquivo com a solucao dos problemas.
                   Autor: Glauber Cintra - 28Abr1993 }
uses crt,dos;

const maxlin = 50;
      maxcol = 50;

type mat = array [1..maxlin,1..maxcol] of comp;
     vetreal = array [1..maxcol] of real;
     vet = array [1..maxcol] of comp;

var A,G : mat;
    C,B,P,Q : vet;
    X,R,Y,S : vetreal;
    arq,saida : text;
    arqsaida : string;
    m,n,res,it,cont : integer;

procedure Le_Problema(var A:mat;var B,C:vet;var m,n:integer);
  { Le, a partir da posicao corrente da cabeca de leitura, o Plin contido
    em 'arq', que precisa estar aberto, devolvendo Plin(A,b,c), m e n }
  var i,j:byte;
  begin
  read(arq,m);
  readln(arq,n);
  for i:=1 to m do
    begin
    for j:=1 to n do
      read(arq,A[i,j]);
    readln(arq,b[i])
    end;
  for i:=1 to n do
    read(arq,c[i])
  end; { Fim de Le_Problema }


{ Predicados }

function Inviavel(var A,G:mat;var B:vet;var S:vetreal;var m,n:integer;
                      pivoant:comp):boolean;
  { Verifica se existe inviabilidade e em caso positivo retorna o
    vetor de inviabilidade G[m,] em S }
  var maior_que_0,menor_que_0:boolean;
      i,j,t:byte;
  begin
  inviavel:=false;
  for i:=1 to m do
    begin
    maior_que_0:=true;
    menor_que_0:=true;
    for j:=1 to n do
      if A[i,j]>0 then
        menor_que_0:=false
        else
        if A[i,j]<0 then
          maior_que_0:=false;
    if (maior_que_0 and (b[i]<0)) then
      begin
      Inviavel:=true;
      for t:=1 to m do
        S[t]:=G[i,t]/pivoant;
      i:=m
      end
      else
      if (menor_que_0 and (b[i]>0)) then
        begin
        Inviavel:=true;
        for t:=1 to m do
          S[t]:=-G[i,t]/pivoant;
        i:=m
        end
    end;
  end; { Fim de Inviavel }


function Ilimitado(var A:mat;var C,q:vet;var R:vetreal;var m,n:integer;
                       pivoant:comp):boolean;
  { Verifica se existe ilimitacao e em caso positivo retorna o vetor
    de ilimitacao em R }
  var menor_que_0:boolean;
      i,j,t:byte;
  begin
  ilimitado:=false;
  for j:=1 to n do
    if c[j]>0 then
      begin
      menor_que_0:=true;
      for i:=1 to m do
        if A[i,j]>0 then
          begin
          menor_que_0:=false;
          i:=m
          end;
      if menor_que_0 then
        begin
        ilimitado:=true;
        for t:=1 to n do
          if t=j then
            R[t]:=1
            else
            if Q[t]<>0 then
              R[t]:=-A[trunc(q[t]),j]/pivoant
              else
              R[t]:=0;
        j:=n
        end
      end
  end; { Fim de Ilimitado }


function M_menos_P_Nulo(var A:mat;var B,P,Q:vet;var m,n:integer):boolean;
  { Verifica se todas as linhas fora de bijecao sao nulas }
  var i,j : byte;
  begin
  M_menos_P_Nulo:=true;
  for i:=1 to m do
    if P[i]=0 then
      for j:=1 to n do
        if abs(A[i,j])>0 then
          begin
          M_menos_P_nulo:=false;
          i:=m;
          j:=n
          end;
  end; { Fim de M_menos_P_Nulo }


function C_Menor_Que_Zero(var c:vet;n:integer):boolean;
  { Verifica se c<=0 }
  var i : byte;
  begin
  c_menor_que_zero:=true;
  for i:=1 to n do
    if c[i]>0 then
      begin
      c_menor_que_zero:=false;
      i:=n
      end
  end;


{ Outras procedures }

procedure pivotacao(var A,G:mat;var b,c,p,q:vet;var y:vetreal;
                        var m,n:integer;fase1:boolean;var pivoant:comp);
  { Executa uma pivotacao do simplex }
  var achoupivo,menor : boolean;
      i,j,l,h,k,t,mm : byte;
      sinal : shortint;
      aik,ck : comp;
  begin
  { Escolhe o pivo usando regra de Bland }
  achoupivo:=false;
  if pivoant>0 then sinal:=1 else sinal:=-1;

  if fase1 then
    { Procurando pivo fora da bijecao }
    for i:=1 to m do
      if P[i]=0 then
        for j:=1 to n do
          if Q[j]=0 then
            if ((A[i,j]>0) and (b[i]>=0)) or ((A[i,j]<0) and (b[i]<=0)) then
              begin
              menor:=true;
              for l:=1 to m do
                if (P[l]<>0) and (sinal*A[l,j]>0) and (b[i]/A[i,j]>b[l]/A[l,j])
                  then
                  begin
                  menor:=false;
                  l:=m
                  end;
              if menor then
                begin
                achoupivo:=true;
                h:=i;
                k:=j;
                P[h]:=k;
                Q[k]:=h;
                i:=m;
                j:=n
                end;
              end;

  if not achoupivo then
    { Procurando o pivo nas linhas de bijecao }
    begin
    if fase1 then
      begin
      for i:=1 to m do
        if P[i]=0 then
           begin
           mm:=i;
           i:=m
           end;
      for i:=1 to m do
        if P[i]<>0 then
          for j:=1 to n do
            if (((A[mm,j]>0) and (b[mm]>=0)) or ((A[mm,j]<0) and
              (b[mm]<=0))) and (Q[j]=0) then
              begin
              k:=j;
              j:=n;
              i:=m
              end;
      end
      else
      for i:=1 to m do
        if P[i]<>0 then
          for j:=1 to n do
            if (((A[i,j]>0) and (b[i]>=0)) or ((A[i,j]<0) and
              (b[i]<=0))) and (Q[j]=0) and (c[j]>0) then
              begin
              k:=j;
              j:=n;
              i:=m
              end;
    for i:=1 to m do
      if P[i]<>0 then
        if A[i,k]>0 then
          begin
          menor:=true;
          h:=i;
          for l:=1 to m do
            if (P[l]<>0) and (sinal*A[l,k]>0) and (b[i]/A[i,k]>b[l]/A[l,k])
              then
              begin
              menor:=false;
              l:=m
              end
              else
              if (P[l]<>0) and (l<>h) and (P[l]<P[h]) and (sinal*A[l,k]>0)
                and (b[i]/A[i,k]=b[l]/A[l,k]) then
                h:=l;
          if menor then
            begin
            achoupivo:=true;
            Q[trunc(P[h])]:=0;
            P[h]:=k;
            Q[k]:=h;
            i:=m;
            j:=n
            end;
          end;
    end;
  if A[h,k]<0 then
    begin
    for j:=1 to n do
      begin
      A[h,j]:=-A[h,j];
      if j<=m then
        G[h,j]:=-G[h,j];
      end;
    b[h]:=-b[h]
    end;
  write(saida,h,',',k,' ');

  { Executa pivotacao }
  for i:=1 to m do
    if i<>h then
      begin
      aik:=A[i,k];
      for j:=1 to n do
        begin
        A[i,j]:=(A[h,k]*A[i,j]-aik*A[h,j])/pivoant;
        if j<=m then
          G[i,j]:=(A[h,k]*G[i,j]-aik*G[h,j])/pivoant
        end;
      b[i]:=(A[h,k]*b[i]-aik*b[h])/pivoant
      end;
  ck:=c[k];
  for j:=1 to n do
    begin
    c[j]:=(A[h,k]*c[j]-A[h,j]*ck)/pivoant;
    if j<=m then
      y[j]:=(A[h,k]*y[j]-G[h,j]*ck)/pivoant
    end;
  pivoant:=A[h,k];
  end;

procedure Resolve_Plin(var A:mat;var B,C:vet;var X,R,Y,S:vetreal;
                           var m,n,res,it:integer);
  { Recebe A,b,c e devolve o numero de iteracoes 'it', x,y,r,s sendo x a
    solucao do Plin, y a solucao do Lpin, r o vetor de ilimitacao do Plin,
    s a vetor de inviabilidade do Plin, e ainda, se Plin e' inviavel res=1,
    se e' ilimitado res=2 e se e' viavel res=3 }
  var i,j : byte;
      pivoant : comp;
  begin
  res:=3;
  it:=1;
  pivoant:=1;

  { Inicializando G,P,Q,y e Verificando se ha' inviabilidade }
  for i:=1 to m do
    begin
    P[i]:=0;
    y[i]:=0;
    for j:=1 to m do
      if i=j then
        g[i,j]:=1
        else
        g[i,j]:=0
    end;
  if inviavel(A,G,b,s,m,n,pivoant) then
    res:=1
    else
    begin
    for i:=1 to n do
      Q[i]:=0;;
    end;

  { Primeira Fase do Simplex }
  while (res=3) and not M_menos_P_nulo(A,b,P,Q,m,n) do
    begin
    pivotacao(A,G,b,c,P,Q,y,m,n,true,pivoant);
    inc(it);
    if inviavel(A,G,b,s,m,n,pivoant) then
      res:=1
    end;

  if (res=3) and ilimitado(A,c,q,r,m,n,pivoant) then
    res:=2;

  { Segunda Fase do Simplex }
  while (res=3) and not c_menor_que_zero(c,n) do
    begin
    pivotacao(A,G,b,c,P,Q,y,m,n,false,pivoant);
    inc(it);
    if ilimitado(A,c,q,r,m,n,pivoant) then
      res:=2
    end;
  if res<>1 then
    begin
    for j:=1 to n do
      if Q[j]<>0 then
        x[j]:=b[trunc(Q[j])]/pivoant
        else
        x[j]:=0;
    for j:=1 to m do
      y[j]:=-y[j]/pivoant
    end

  end; { Fim de Resolve_Plin }


procedure Grava_Solucao(var res,it,cont:integer;var x,r,y,s:vetreal);
  var i : byte;
  begin
  case res of
    1:
      begin
      writeln(saida,'inviavel');
      write(saida,'s = ');
      for i:=1 to m do
        write(saida,S[i]:3:3,'  ')
      end;
    2:
      begin
      writeln(saida,'ilimitado');
      write(saida,'x = ');
      for i:=1 to n do
        write(saida,x[i]:3:3,'  ');
      writeln(saida);
      write(saida,'r = ');
      for i:=1 to n do
        write(saida,r[i]:3:3,'  ')
      end;
    3:
      begin
      writeln(saida,'soluvel');
      write(saida,'x = ');
      for i:=1 to n do
        write(saida,x[i]:3:3,'  ');
      writeln(saida);
      write(saida,'y = ');
      for i:=1 to m do
        write(saida,y[i]:3:3,'  ')
      end
    end;
  writeln(saida);writeln(saida)
  end; { Fim de Grava_Solucao }


{ Programa Principal }

begin
if paramcount=0 then
  begin
  writeln('Sintaxe: simplex <arquivo-entrada> [<arquivo-saida>]');
  writeln('se <arquivo-saida> nao for especificado "saida" sera assumido');
  halt
  end
  else
  if paramcount=1 then
    arqsaida:='saida1'
    else
    arqsaida:=paramstr(2);
assign(arq,paramstr(1));
assign(saida,arqsaida);
{$I-}
reset(arq);
rewrite(saida);
if ioresult <> 0 then
  begin
  writeln('Nao foi possivel abrir arquivo. Execucao abortada.');
  halt
  end;
{$I+}
writeln;writeln('Solucionador de problemas de programacao linear do tipo:');
writeln('Max c.x   sujeito A.x = b  e x >= 0');
writeln('Autor: Glauber Cintra  -  4Jun1993');writeln;
writeln(saida,'Glauber Cintra,  4Jun1993');writeln(saida);
cont:=1;
while not eof(arq) do
  begin
  le_problema(A,B,C,m,n);
  if (m>0) and (n>0) then
    begin
    writeln('Problema ',cont,' sendo resolvido. Aguarde...');
    writeln(saida,'Problema ',cont);
    write(saida,'Sequencia de Pivos: ');
    resolve_plin(A,B,C,x,r,y,s,m,n,res,it);
    writeln(saida);
    write(saida,it,' iteracoes, ');
    grava_solucao(res,it,cont,x,r,y,s);
    inc(cont)
    end
  end;
close(arq);
close(saida)
end. { Fim do programa principal }