File : lin-ops.adb

package body Lin.Ops is

  function CompareR(U,V: Trip) return Integer is
--- 1 means U<V, -1 means U>V
  begin
    if U.R<V.R then return  1; end if;
    if U.R>V.R then return -1; end if;
    return 0;
  end CompareR;

  function CompareIJ(U,V: Trip) return Integer is
--- 1 means U<V, -1 means U>V
  begin
    if U.I<V.I then return  1; end if;
    if U.I>V.I then return -1; end if;
    if U.J<V.J then return  1; end if;
    if U.J>V.J then return -1; end if;
    return 0;
  end CompareIJ;

  function CompareJI(U,V: Trip) return Integer is
--- 1 means U<V, -1 means U>V
  begin
    if U.J<V.J then return  1; end if;
    if U.J>V.J then return -1; end if;
    if U.I<V.I then return  1; end if;
    if U.I>V.I then return -1; end if;
    return 0;
  end CompareJI;

  procedure Sort(W: in out Sparse; Compare: in TripCompare) is
--- translated from numerical recipes
--- elements in sorted W are increasing
    L,J,Ir,I: Integer;
    Rw: Trip;
  begin
    Ir := W(0).I;
    if Ir<2 then return; end if;
    L := Ir/2+1;
    loop
      if L>1 then
        L := L-1;
        Rw := W(L);
      else
        Rw := W(Ir);
        W(Ir) := W(1);
        Ir := Ir-1;
        if Ir=1 then
          W(1) := Rw;
          return;
        end if;
      end if;
      I := L;
      J := L+L;
      while J<=Ir loop
        if J<Ir and then Compare(W(J),W(J+1))>0 then
          J := J+1;
        end if;
        if Compare(Rw,W(J))>0 then
          W(I) := W(J);
          I := J;
          J := J+J;
        else
          J := Ir+1;
        end if;
      end loop;
      W(I) := Rw;
    end loop;
  end Sort;

  procedure Cleanup(W: in out Sparse; Force: in Boolean := False) is
  begin
    if not Force and then W(0).J=0 then return; end if;
    if W(0).I=0 then return; end if;
    declare
      Changed: Boolean := False;
      L: Integer := 0;
      Keep: array(W'Range) of Boolean;
    begin
      Sort(W,CompareIJ'Access);
      for J in reverse 1 .. Length(W) loop
        if W(J).R=Zero then
          Keep(J) := False;
          Changed := True;
        elsif J>1 and then CompareIJ(W(J),W(J-1))=0 then
          W(J-1).R := W(J-1).R+W(J).R;
          Keep(J) := False;
          Changed := True;
        else
          Keep(J) := True;
        end if;
      end loop;
      if Changed then
        for J in 1 .. Length(W) loop
          if Keep(J) then
            L := L+1;
            W(L) := W(J);
          end if;
        end loop;
        W(0).I := L;
      end if;
      W(0).J := 0;
    end;
  end Cleanup;

end Lin.Ops;