
(*        
   INPUT	: The associated datafile for Maxflow algorithm is called
		  "MaxflowDatafile".

		  "MaxflowDatafile" consists of 
		  1.  # of nodes in a directed network,N,
		  2.  # od edges in the network,M,
		  3.  WEIGHT MATRIX of NxN

		  FIRST NUMBER in "MaxflowDatafile" is # of nodes, N.

		  SECOND NUMBER is # of edges, M.

		  THIRD set of data is the WEIGHT MATRIX.  Entry (i,j)
	   	  has max flow of k units.  Non existent  edge has
		  max flow value of 0

   Algorithm	: This maxflow algorithm finds a flow pattern that
		  yields the maximum flow value from a specified node
		  S (NODE 1) TO another specified node T(node N) in a
		  given weighted network with N nodes and M edges.

		  Maximum # of nodes in a network is set to be maxnode= 50
		  Maximum # of edges is set to maxedge = 50.
	  	  One can modify them according to their needs.

   Output	: Output is a NxN maximum flow pattern.

   NOTE		: The data structure used here is the WEIGHT MATRIX.
		  If the input network is SPARSE, it would be better
		  to use another data structure such as FORWARD STAR.  *)


program MaximumFlow (input,output, MaxflowDatafile,MaxflowOutfile);

const	maxnode   =   50;
        maxedge   =   50;
        maxelement = maxedge*3;

type	CHARFILE =  file of char;
        ARRNN	 =  array[1..maxnode,1..maxnode] of integer;
        ARRN	 =  array[1..maxnode] of integer;
        ARRMAX   =  array[1..maxelement] of integer;

var	N,  M,  Nextint	: 	integer;
	S,	T	:	integer;
	MAX		:	integer;
	CAPA		:	ARRNN;
	FLOW		:	ARRNN;
        MaxflowDatafile :       CHARFILE;
        MaxflowOutfile  :	CHARFILE;



procedure Infile (var  N,  M,  Nextint  : integer;
		  var  CAPA		: ARRNN);

var row, column : integer;

begin
  reset(MaxflowDatafile);
  readln(MaxflowDatafile, Nextint);
  N := Nextint;
  readln(MaxflowDatafile, Nextint);
  M := Nextint;
  for row := 1 to N do
  begin
    for column := 1 to N do
    begin
      read(MaxflowDatafile, Nextint);
      CAPA[row,column] := Nextint;
    end;
    readln(MaxflowDatafile);
  end;
end;




procedure MAXFLOW(
       N,S,T,MAX:integer;
   var CAPA,FLOW:ARRNN);

   var LCAPA,LFLOW,EL         :ARRNN;
       LABELS,LAYERS,PRED,SUCC:ARRN;
       DATA,LINK              :ARRMAX;
       STPATH                 :boolean;
       AV,I,J,L,Dif           :integer;

   procedure ADD(var A:ARRN;I,X:integer);
      var P    :integer;
          FOUND:boolean;
   begin
       FOUND:=false;
       P:=A[I];
       while (not FOUND) and (P <> 0) do
          if DATA[P] = X then FOUND:=true
          else P:=LINK[P];
       if not FOUND then begin
          P:=AV;  AV:=LINK[AV];
          DATA[P]:=X;  LINK[P]:=A[I];
          A[I]:=P
       end  (* if not FOUND *)
   end;  (* ADD *)

   procedure DELETE(var A:ARRN;I,X:integer);
      var P,Q:integer;
   begin
      P:= A[I];
      if DATA[P] = X then A[I]:=LINK[P]
      else begin
         while DATA[P] <> X do begin
            Q:=P;  P:=LINK[P]
         end;
         LINK[Q]:=LINK[P]
      end
   end;  (* DELETE *)

   procedure LAYER(var L:integer;var LCAPA,EL:ARRNN;
                   var STPATH:boolean);
      var I,J,P,Q,X,Y,V,W:integer;
   begin
      STPATH:=true;
      for V:=1 to N do begin
         LABELS[V]:=-1;
         SUCC[V]:=0;  PRED[V]:=0
      end;
      P:=AV;  AV:=LINK[AV];
      DATA[P]:=S;  LINK[P]:=0;
      LAYERS[1]:=P;  LABELS[S]:=1;
      for X:=1 to N do
         for Y:=1 to N do LCAPA[X,Y]:=CAPA[X,Y];
      for X:=1 to N do
         for Y:=1 to N do begin
            EL[X,Y]:=0;
            if FLOW[X,Y] > 0 then begin
               LCAPA[X,Y]:=CAPA[X,Y]-FLOW[X,Y];
               LCAPA[Y,X]:=FLOW[X,Y]+LCAPA[Y,X]
            end
         end;  (* for Y, X *)                   (* INITIALIZATION OVER *)
      I:=1;                        (* forWARD TRAVERSAL and LABELING *)
      while (LAYERS[I] <> 0) and (LABELS[T] = -1) do begin
         LAYERS[I+1]:=0;  P:=LAYERS[I];
         while P <> 0 do begin
            X:=DATA[P];
            for Y:=1 to N do
               if ((LABELS[Y] = -1) or (LABELS[Y] = I+1))
                  and (LCAPA[X,Y] > 0) then begin
                  ADD(LAYERS,I+1,Y);
                  LABELS[Y]:=I+1;
                  ADD(SUCC,X,Y);  ADD(PRED,Y,X);
                  EL[X,Y]:=1
               end  (* if ((LABELS[Y] = -1) ... *)
               else LCAPA[X,Y]:=0;
            P:=LINK[P]
         end;  (* while P <> 0 *)
         I:=I+1
      end;  (* while (LAYERS[I] <> 0) ... *)
      L:=I;                        (* BACKWARD TRAVERSAL and PRUNING *)
      if LABELS[T] = -1 then STPATH:=false
      else begin
         J:=I;
         while J <> 1 do begin
            P:=LAYERS[J];
            while P <> 0 do begin
               W:=DATA[P];
               if (SUCC[W] = 0) and (W <> T) then begin
                  Q:=PRED[W];
                  while Q <> 0 do begin
                     X:=DATA[Q];  EL[X,W]:=0;
                     DELETE(SUCC,X,W);
                     LCAPA[X,W]:=0;  Q:=LINK[Q]
                  end;
                  DELETE(LAYERS,J,W);
                  PRED[W]:=0;  LABELS[W]:=-1
               end;  (* if (SUCC[W] = 0) ... *)
               P:=LINK[P]
            end;  (* while P <> 0 *)
            J:=J-1
         end  (* while J <> 1 *)
      end  (* else: LABELS[T] <> -1 *)
   end;  (* LAYER *)

   procedure SATURATE(var L:integer;var LFLOW,EL:ARRNN);
      var INPOT,OUTPOT,POTEN,INFLOW,OUTFLOW:ARRN;
          V,X,Y,I,J,K,P,R,RLAYER           :integer;
          FLAG                             :boolean;

      procedure REFNODE(var L,R,RLAYER:integer);
         var V:integer;
      begin
         POTEN[S]:=OUTPOT[S];  POTEN[T]:=INPOT[T];
         if POTEN[S] < POTEN[T] then  begin
            R:=S;  RLAYER:=1
         end
         else  begin
            R:=T;  RLAYER:=L
         end;
         for V:=1 to N do
            if (LABELS[V] <> -1) and (V <> S) and (V <> T) then
            begin
               if INPOT[V] < OUTPOT[V] then POTEN[V]:=INPOT[V]
               else POTEN[V]:=OUTPOT[V];
               if POTEN[V] < POTEN[R] then  begin
                  R:=V;
                  RLAYER:=LABELS[V]
               end
            end  (* if (LABELS[V] <> -1) ..., for V *)
      end;  (* REFNODE *)

      procedure PUSH(I:integer);
         var P,Q,U,V,AVACAP:integer;
      begin
         P:=LAYERS[I];
         while P <> 0 do begin
            U:=DATA[P];
            if OUTFLOW[U] > 0 then begin
               Q:=SUCC[U];
               while (OUTFLOW[U] > 0) and (Q <> 0) do begin
                  V:=DATA[Q];
                  if LABELS[V] <> -1 then begin
                     AVACAP:=LCAPA[U,V]-LFLOW[U,V];
                     if AVACAP > 0 then begin
                        (* Send MIN(AVACAP,OUTFLOW[U]) THROUGH (U,W) *)
                        if AVACAP > OUTFLOW[U] then
                           AVACAP:=OUTFLOW[U];
                        LFLOW[U,V]:=LFLOW[U,V]+AVACAP;
                        OUTFLOW[U]:=OUTFLOW[U]-AVACAP;
                        OUTFLOW[V]:=OUTFLOW[V]+AVACAP;
                        OUTPOT[U]:=OUTPOT[U]-AVACAP;
                        INPOT[V]:=INPOT[V]-AVACAP
                     end  (* if AVACAP > 0 *)
                  end;  (* if LABELS[V] <> -1 *)
                  Q:=LINK[Q]
               end  (* while (OUTFLOW[U] > 0) ... *)
            end;  (* if OUTFLOW[U] > 0 *)
            P:=LINK[P]
         end  (* while P <> 0 *)
      end;  (* PUSH *)

      procedure PULL(J:integer);
         var P,Q,U,V,AVACAP:integer;
      begin
         P:=LAYERS[J];
         while P <> 0 do begin
            U:=DATA[P];
            if INFLOW[U] > 0 then begin
               Q:=PRED[U];
               while (INFLOW[U] > 0) and (Q <> 0) do begin
                  V:=DATA[Q];
                  if LABELS[V] <> -1 then begin
                     AVACAP:=LCAPA[V,U]-LFLOW[V,U];
                     if AVACAP > 0 then begin
                         (* Send MIN(AVACAP,INFLOW[U]) THROUGH (V,U) *)
                        if AVACAP > INFLOW[U] then
                           AVACAP:=INFLOW[U];
                        LFLOW[V,U]:=LFLOW[V,U]+AVACAP;
                        INFLOW[U]:=INFLOW[U]-AVACAP;
                        INFLOW[V]:=INFLOW[V]+AVACAP;
                        OUTPOT[V]:=OUTPOT[V]-AVACAP;
                        INPOT[U]:=INPOT[U]-AVACAP
                     end  (* if AVACAP > 0 *)
                  end;  (* if LABELS[V] <> -1 *)
                  Q:=LINK[Q]
               end  (* while (INFLOW[U] > 0) ... *)
            end;  (* if INFLOW[U] > 0 *)
            P:=LINK[P]
         end  (* while P <> 0 *)
      end;  (* PULL *)
   begin                                         (* BODY OF SATURATE *)
      for V:=1 to N do begin
         INPOT[V]:=0;  OUTPOT[V]:=0
      end;
      for I:=1 to N do
         for J:=1 to N do
            if EL[I,J] = 1 then begin
               INPOT[J]:=INPOT[J]+LCAPA[I,J];
               OUTPOT[I]:=OUTPOT[I]+LCAPA[I,J]
            end;  (* for J, I *)
      for I:=1 to N do
         for J:=1 to N do LFLOW[I,J]:=0;
      FLAG:=true;                             (* INITIALIZATION OVER *)
      while FLAG do begin                 (* while GL IS UNSATURATED *)
         REFNODE(L,R,RLAYER);
         if POTEN[R] <> 0 then begin
            INFLOW[R]:=POTEN[R];  OUTFLOW[R]:=POTEN[R];
            for V:=1 to N do
               if (LABELS[V] <> -1) and (V <> R) then begin
                  INFLOW[V]:=0;  OUTFLOW[V]:=0
               end;
            for K:=RLAYER to L-1 do PUSH(K);
            for J:=RLAYER downto 2 do PULL(J)
         end;  (* if POTEN[R] <> 0 *)
         if (POTEN[R] <> 0) or ((R <> S) and (R <> T)) then begin
            LABELS[R]:=-1;  P:=SUCC[R];
            while P <> 0 do begin
               Y:=DATA[P];
               INPOT[Y]:=INPOT[Y]-(LCAPA[R,Y]-LFLOW[R,Y]);
               EL[R,Y]:=0;
               DELETE(SUCC,R,Y);  DELETE(PRED,Y,R);
               P:=LINK[P]
            end;  (* while P <> 0 *)
            P:=PRED[R];
            while P <> 0 do begin
               X:=DATA[P];
               OUTPOT[X]:=OUTPOT[X]-(LCAPA[X,R]-LFLOW[X,R]);
               EL[X,R]:=0;
               DELETE(SUCC,X,R);  DELETE(PRED,R,X);
               P:=LINK[P]
            end  (* while P <> 0 *)
         end  (* if (POTEN[R] <> 0) ... *)
         else FLAG:=false
      end  (* while FLAG *)
   end;  (* SATURATE *)

   procedure INITIALIZE(MAX:integer);
      var I:integer;
   begin
      for I:=1 to MAX do LINK[I]:=I+1;
      AV:=1
   end;  (* INITIALIZE *)

begin                                                   (* MAIN BODY *)
   for I:=1 to N do
      for J:=1 to N do FLOW[I,J]:=0;
   INITIALIZE(MAX);                           (* INITIALIZATION OVER *)
   LAYER(L,LCAPA,EL,STPATH);
   while STPATH do begin
      SATURATE(L,LFLOW,EL);
      for I:=1 to N do
         for J:=1 to N do
            if LFLOW[I,J] > 0 then begin
               Dif:=LFLOW[I,J]-FLOW[J,I];
               if Dif > 0 then begin
                  FLOW[I,J]:=FLOW[I,J]+Dif;
                  FLOW[J,I]:=0
               end
               else FLOW[J,I]:=-Dif
            end;  (* if LFLOW[I,J] > 0, for J, I *)
      INITIALIZE(MAX);
      LAYER(L,LCAPA,EL,STPATH)
   end  (* while STPATH *)
end;  (* MAXFLOW *)





procedure Outfile(N	: integer;
		  FLOW	: ARRNN   );

var row, column : integer;

begin
  rewrite(MaxflowOutfile);
  writeln(MaxflowOutfile,' Max-Flow pattern is ');
  for row := 1 to N do
  begin
    for column := 1 to N do
       write(MaxflowOutfile,FLOW[row,column]:5);
    writeln(MaxflowOutfile);
  end;
end;




begin (* main  *)
  Infile(N,M,Nextint,CAPA);
  S := 1;
  T := N;
  MAX := 3*M;
  MAXFLOW(N,S,T,MAX,CAPA,FLOW);
  Outfile(N,FLOW);
end.
