{$M 65000,0,0}
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 zero = 1e-10;
      maxlin = 50;
      maxcol = 50;

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

var A,G : mat;
    C,X,R,B,Y,S,P,Q : vet;
    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.
    O lay-out do arquivo de entrada, que é:

m n
a_{1,1} a_{1, 2} ... a_{1,n} b_1
a_{2,1} a_{2, 2} ... a_{2,n} b_2
...
a_{m,1} a_{m, 2} ... a_{m,n} b_m
c_1 c_2 ... c_n

  onde m é o número de restrições, n o numero de variáveis, a_{i, j} o
  coeficiente da variável j na restrição i, b_i o lado direito da restrição i
  e c_j o coeficiente da variável j na função objetivo. Está implícito,
  conforme o comentário no início do programa, que a função é de maximização,
  as restrições são todas do tipo = e as variáveis são todas não negativas.

  Por exemplo, o problema:

maximizar 2x - 3y
sujeito a:
  x + 4y = 17
 2x - y = 3.5
  x + y = 12
  x >= 0
  y >= 0

  seria representado por um arquivo de texto com o seguinte conteúdo:

3 2
1 4 17
2 -1 3.5
1 1 12
2 -3

}

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,S:vet;var m,n:integer):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]>zero then
        menor_que_0:=false
        else
        if A[i,j]<-zero then
          maior_que_0:=false;
    if (maior_que_0 and (b[i]<-zero)) then
      begin
      Inviavel:=true;
      for t:=1 to m do
        S[t]:=G[i,t];
      i:=m
      end
      else
      if (menor_que_0 and (b[i]>zero)) then
        begin
        Inviavel:=true;
        for t:=1 to m do
          S[t]:=-G[i,t];
        i:=m
        end
    end;
  end; { Fim de Inviavel }


function Ilimitado(var A:mat;var C,R,q:vet;var m,n:integer):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]>zero 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]
              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])>zero 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 forma_bijecao(var A:mat;var B,P,Q:vet;var m,n:integer);
  { Coloca em P e Q linhas e colunas tal que A[P,Q] e' de bijecao e
    b[P]>=0 }
  var i,j,lin_1 : byte;
      tem_1_1,tem_zeros:boolean;
  begin
  for j:=1 to n do
    begin
    tem_zeros:=true;
    tem_1_1:=false;
    for i:=1 to m do
      if (A[i,j]=1) and (c[j]=0) then
        if tem_1_1 then
          begin
          tem_1_1:=false;
          i:=m
          end
          else
          begin
          tem_1_1:=true;
          lin_1:=i
          end
        else
        if A[i,j]<>0 then
          begin
          tem_zeros:=false;
          i:=m
          end;
    if tem_1_1 and tem_zeros and (P[lin_1]=0) and (b[lin_1]>=0) then
      begin
      P[lin_1]:=j;
      Q[j]:=lin_1
      end
    end;
  end; { Fim de Forma_bijecao }


procedure pivotacao(var A,G:mat;var b,c,y,p,q:vet;var m,n:integer;fase1:boolean);
  { Executa uma pivotacao do simplex }
  var achoupivo,menor : boolean;
      i,j,l,h,k,t,mm : byte;
      pivo,teta : real;
  begin
  { Escolhe o pivo usando regra de Bland }
  achoupivo:=false;

  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]>zero) and (b[i]>=0)) or ((A[i,j]<-zero) and (b[i]<=0)) then
              begin
              menor:=true;
              for l:=1 to m do
                if (P[l]<>0) and (A[l,j]>zero) 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]>zero) and (b[mm]>=0)) or ((A[mm,j]<-zero) 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]>zero) and (b[i]>=0)) or ((A[i,j]<-zero) 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]>zero then
          begin
          menor:=true;
          h:=i;
          for l:=1 to m do
            if (P[l]<>0) and (A[l,k]>zero) 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 (A[l,k]>zero)
                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;
  pivo:=A[h,k];

  { Executa pivotacao }
  for i:=1 to m do
    if i<>h then
      begin
      teta:=A[i,k]/pivo;
      for j:=1 to n do
        begin
        if j=k then
          A[i,j]:=0
          else
          A[i,j]:=A[i,j]-A[h,j]*teta;
        if j<=m then
          G[i,j]:=G[i,j]-G[h,j]*teta
        end;
      b[i]:=b[i]-b[h]*teta
      end;
  teta:=c[k]/pivo;
  for j:=1 to n do
    begin
    if j=k then
      c[j]:=0
      else
      c[j]:=c[j]-A[h,j]*teta;
    if j<=m then
      y[j]:=y[j]-G[h,j]*teta
    end;
  for j:=1 to n do
    begin
    if j=k then
      A[h,j]:=1
      else
      A[h,j]:=A[h,j]/pivo;
    if j<=m then
      G[h,j]:=G[h,j]/pivo
    end;
  b[h]:=b[h]/pivo
  end;


procedure Resolve_Plin(var A:mat;var B,Y,S,C,X,R:vet;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;
  begin
  res:=3;
  it:=1;

  { Verificando se ha' inviabilidade ou montando a bijecao }
  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) then
    res:=1
    else
    begin
    for i:=1 to n do
      Q[i]:=0;
    forma_bijecao(A,b,P,Q,m,n)
    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,y,P,Q,m,n,true);
    inc(it);
    if inviavel(A,G,b,s,m,n) then
      res:=1
    end;

  if (res=3) and ilimitado(A,c,r,q,m,n) 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,y,P,Q,m,n,false);
    inc(it);
    if ilimitado(A,c,r,q,m,n) 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])]
        else
        x[j]:=0;
    for j:=1 to m do
      y[j]:=-y[j]
    end

  end; { Fim de Resolve_Plin }


procedure Grava_Solucao(var res,it,cont:integer;var x,r,y,s:vet);
  var i : byte;
  begin
  write(saida,'Problema ',cont,': ',it,' iteracoes, ');
  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:='saida'
    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  -  28Abr1993');writeln;
writeln(saida,'Glauber Cintra,  28Abr1993');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...');
    resolve_plin(A,B,y,s,C,x,r,m,n,res,it);
    grava_solucao(res,it,cont,x,r,y,s);
    inc(cont)
    end
  end;
close(arq);
close(saida)
end. { Fim do programa principal }