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;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>