Return to standard_integer_linear_equalities.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices

Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/standard_integer_linear_equalities.adb, Revision 1.1.1.1

1.1       maekawa     1: with Standard_Common_Divisors;           use Standard_Common_Divisors;
                      2: with Standard_Integer_Linear_Solvers;    use Standard_Integer_Linear_Solvers;
                      3:
                      4: package body Standard_Integer_Linear_Equalities is
                      5:
                      6:   procedure Triangulate ( l : in natural; m : in matrix;
                      7:                           index : in natural; ineq : in out matrix ) is
                      8:
                      9:     a,b,lcmab,faca,facb : integer;
                     10:
                     11:   begin
                     12:     if ineq(index,l) /= 0
                     13:      then a := m(l,l); b := ineq(index,l);
                     14:           lcmab := lcm(a,b);
                     15:           if lcmab < 0 then lcmab := -lcmab; end if;
                     16:           faca := lcmab/a;  facb := lcmab/b;
                     17:           if facb < 0
                     18:            then facb := -facb; faca := -faca;
                     19:           end if;
                     20:           for j in l..ineq'last(2) loop
                     21:             ineq(index,j) := facb*ineq(index,j) - faca*m(l,j);
                     22:           end loop;
                     23:     end if;
                     24:   end Triangulate;
                     25:
                     26:   procedure Triangulate ( l : in natural; m : in matrix;
                     27:                           first,last : in natural; ineq : in out matrix ) is
                     28:   begin
                     29:     for i in first..last loop
                     30:       if ineq(i,l) /= 0
                     31:        then Triangulate(l,m,i,ineq);
                     32:       end if;
                     33:     end loop;
                     34:   end Triangulate;
                     35:
                     36:   procedure Triangulate ( index,start : in natural; ineq : in out matrix ) is
                     37:
                     38:     column : natural := start;  -- current column counter
                     39:     firstineq : natural := ineq'first(1);
                     40:     found : boolean;
                     41:     ind2 : natural;
                     42:     a,b,lcmab,faca,facb : integer;
                     43:
                     44:   begin
                     45:     --put_line("INEQUALTITY");
                     46:     --put(" BEFORE UPDATE : ");
                     47:     --for l in ineq'range(2) loop
                     48:     --  put(ineq(index,l),1); put(' ');
                     49:     --end loop;
                     50:     --new_line;
                     51:     loop
                     52:      -- SEARCH FOR FIRST NONZERO ENTRY IN CURRENT INEQUALITY :
                     53:       while column < ineq'last(2) and then ineq(index,column) = 0 loop
                     54:         column := column + 1;
                     55:       end loop;
                     56:       exit when ((ineq(index,column) = 0) or else (column = ineq'last(2)));
                     57:                                         -- nothing to eliminate
                     58:      -- SEARCH FOR INEQUALITY,
                     59:      --    WITH SAME FIRST NONZERO COLUMN, BUT WITH OPPOSITE SIGN :
                     60:       found := false;
                     61:       for k in firstineq..(index-1) loop
                     62:         if ineq(index,column)*ineq(k,column) < 0 -- check for sign
                     63:          then found := true;
                     64:               for l in start..column-1 loop      -- check if same zero pattern
                     65:                -- if ineq(index,l) = 0
                     66:                --  then
                     67:                 found := (ineq(k,l) = 0);
                     68:                -- end if;
                     69:                 exit when not found;
                     70:               end loop;
                     71:               if found then ind2 := k; end if;
                     72:         end if;
                     73:         exit when found;
                     74:       end loop;
                     75:       exit when not found;  -- no possibility for elimination
                     76:      -- if found
                     77:      --  then
                     78:         -- ELIMINATE BY MAKING A POSITIVE LINEAR COMBINATION :
                     79:          a := ineq(index,column);  b := ineq(ind2,column);
                     80:          if a < 0
                     81:           then lcmab := lcm(-a,b);
                     82:                faca := lcmab/(-a); facb := lcmab/b;
                     83:           else lcmab := lcm(a,-b);
                     84:                facb := lcmab/(-b); faca := lcmab/a;
                     85:          end if;
                     86:          if ineq(index,ineq'last(2)) >= 0
                     87:            or else  -- PRESERVE SIGN OF ineq(index,ineq'last(2)) !!!!
                     88:                 (faca*ineq(index,ineq'last(2))
                     89:                             + facb*ineq(ind2,ineq'last(2)) < 0)
                     90:           then
                     91:             for l in start..ineq'last(2) loop
                     92:               ineq(index,l) := faca*ineq(index,l) + facb*ineq(ind2,l);
                     93:             end loop;
                     94:            -- PROCEED :
                     95:             column := column + 1;
                     96:             firstineq := ineq'first(1);
                     97:           else
                     98:            -- TRY TO ELIMINATE WITH OTHER INEQUALITIES :
                     99:             firstineq := ind2 + 1;
                    100:          end if;
                    101:          if (firstineq >= index)
                    102:            -- impossible to eliminate with sign preservation
                    103:           then firstineq := ineq'first(2);
                    104:                column := column + 1;
                    105:          end if;
                    106:      --  else
                    107:      --    column := column + 1;
                    108:      --    firstineq := ineq'first(2);
                    109:      -- end if;
                    110:       exit when (column >= ineq'last(2)-1);
                    111:     end loop;
                    112:     --put(" AFTER UPDATE : ");
                    113:     --for l in ineq'range(2) loop
                    114:     --  put(ineq(index,l),1); put(' ');
                    115:     --end loop;
                    116:     --new_line;
                    117:   end Triangulate;
                    118:
                    119: end Standard_Integer_Linear_Equalities;