![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dictionaries.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports|
Return to dictionaries.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports |
1.1 maekawa 1: package body Dictionaries is
2:
3: -- INITIALIZERS :
4:
5: procedure Init_Basis
6: ( in_bas,out_bas : in out Standard_Integer_Vectors.Vector ) is
7:
8: n : constant natural := out_bas'last;
9:
10: begin
11: for i in in_bas'range loop
12: in_bas(i) := n+i; -- slack variable for ith constraint
13: end loop;
14: for i in out_bas'range loop
15: out_bas(i) := i;
16: end loop;
17: end Init_Basis;
18:
19: function Init_Primal_Dictionary
20: ( c : Standard_Floating_Vectors.Vector; a : Matrix )
21: return Matrix is
22:
23: dic : Matrix(0..a'last(1),a'range(2));
24:
25: begin
26: for i in c'range loop
27: dic(0,i) := -c(i);
28: end loop;
29: for i in a'range(1) loop
30: for j in a'range(2) loop
31: dic(i,j) := a(i,j);
32: end loop;
33: end loop;
34: return dic;
35: end Init_Primal_Dictionary;
36:
37: function Init_Dual_Dictionary
38: ( c : Standard_Floating_Vectors.Vector; a : Matrix )
39: return Matrix is
40:
41: dic : Matrix(0..a'last(1),a'range(2));
42:
43: begin
44: for i in c'range loop
45: dic(0,i) := -c(i);
46: end loop;
47: for i in a'range(1) loop
48: for j in a'range(2) loop
49: dic(i,j) := -a(i,j);
50: end loop;
51: end loop;
52: return dic;
53: end Init_Dual_Dictionary;
54:
55: procedure Primal_Init
56: ( c : in Standard_Floating_Vectors.Vector;
57: a : in Matrix; dic : out Matrix;
58: in_bas,out_bas : in out Standard_Integer_Vectors.Vector ) is
59: begin
60: dic := Init_Primal_Dictionary(c,a);
61: Init_Basis(in_bas,out_bas);
62: end Primal_Init;
63:
64: procedure Dual_Init
65: ( c : in Standard_Floating_Vectors.Vector;
66: a : in Matrix; dic : out Matrix;
67: in_bas,out_bas : in out Standard_Integer_Vectors.Vector ) is
68: begin
69: dic := Init_Dual_Dictionary(c,a);
70: Init_Basis(in_bas,out_bas);
71: end Dual_Init;
72:
73: -- MODIFIERS :
74:
75: procedure Primal_Update
76: ( dic : in out Matrix;
77: in_bas,out_bas : in out Standard_Integer_Vectors.Vector;
78: eps : in double_float; unbounded : out boolean ) is
79:
80: column_index : natural := 0;
81: min : double_float := 0.0;
82: row_index : natural := 0;
83: temp : double_float;
84: tmp : integer;
85:
86: begin
87: for i in (dic'first(2)+1)..dic'last(2) loop -- which unknown enters?
88: if dic(0,i) < min
89: then min := dic(0,i); column_index := i;
90: end if;
91: end loop;
92: if column_index /= 0 -- otherwise optimality is reached
93: then
94: for i in (dic'first(1)+1)..dic'last(1) loop -- which unknown leaves?
95: temp := dic(i,column_index);
96: if abs(temp) > eps
97: then temp := dic(i,0)/temp;
98: if temp > 0.0
99: then if row_index = 0
100: then row_index := i; min := temp;
101: elsif temp < min
102: then row_index := i; min := temp;
103: end if;
104: end if;
105: end if;
106: end loop;
107: if row_index = 0 -- solution is unbounded
108: then
109: unbounded := true;
110: else
111: unbounded := false;
112: temp := dic(row_index,column_index); -- pivot
113: for j in dic'range(2) loop -- divide pivot row
114: dic(row_index,j) := dic(row_index,j) / temp;
115: end loop;
116: for i in dic'range(1) loop -- update other rows, except pivot column
117: if i /= row_index
118: then
119: for j in dic'range(2) loop
120: if j /= column_index
121: then
122: dic(i,j) := dic(i,j)-dic(row_index,j)*dic(i,column_index);
123: end if;
124: end loop;
125: end if;
126: end loop;
127: for i in dic'range(1) loop -- update pivot column
128: if i = row_index
129: then dic(i,column_index) := 1.0 / temp;
130: else dic(i,column_index) := - dic(i,column_index) / temp;
131: end if;
132: end loop;
133: tmp := in_bas(row_index); -- update basis information
134: in_bas(row_index) := out_bas(column_index);
135: out_bas(column_index) := tmp;
136: end if;
137: end if;
138: end Primal_Update;
139:
140: procedure Dual_Update
141: ( dic : in out Matrix;
142: in_bas,out_bas : in out Standard_Integer_Vectors.Vector;
143: eps : in double_float; feasible : out boolean ) is
144:
145: column_index : natural := 0;
146: min : double_float := 0.0;
147: row_index : natural := 0;
148: temp : double_float;
149: tmp : integer;
150:
151: begin
152: for i in (dic'first(1)+1)..dic'last(1) loop -- which unknown leaves?
153: if dic(i,0) < min
154: then min := dic(i,0); row_index := i;
155: end if;
156: end loop;
157: if row_index /= 0 -- otherwise optimality is reached
158: then
159: for i in (dic'first(2)+1)..dic'last(2) loop -- which unknown enters?
160: temp := dic(row_index,i);
161: if abs(temp) > eps and then (temp < 0.0)
162: then temp := dic(0,i)/temp;
163: if column_index = 0
164: then column_index := i; min := temp;
165: elsif temp < min
166: then column_index := i; min := temp;
167: end if;
168: end if;
169: end loop;
170: if column_index = 0 -- problem is infeasible
171: then
172: feasible := false;
173: else
174: feasible := true;
175: temp := dic(row_index,column_index); -- pivot
176: for j in dic'range(2) loop -- divide pivot row
177: dic(row_index,j) := dic(row_index,j) / temp;
178: end loop;
179: for i in dic'range(1) loop -- update other rows, except pivot column
180: if i /= row_index
181: then
182: for j in dic'range(2) loop
183: if j /= column_index
184: then
185: dic(i,j) := dic(i,j)-dic(row_index,j)*dic(i,column_index);
186: end if;
187: end loop;
188: end if;
189: end loop;
190: for i in dic'range(1) loop -- update the pivot column
191: if i = row_index
192: then dic(i,column_index) := 1.0 / temp;
193: else dic(i,column_index) := - dic(i,column_index) / temp;
194: end if;
195: end loop;
196: tmp := in_bas(row_index); -- update the basis information
197: in_bas(row_index) := out_bas(column_index);
198: out_bas(column_index) := tmp;
199: end if;
200: end if;
201: end Dual_Update;
202:
203: -- SELECTORS :
204:
205: function Primal_Optimal ( dic : Matrix; eps : double_float ) return boolean is
206: begin
207: for i in (dic'first(2)+1)..dic'last(2) loop
208: if abs(dic(0,i)) > eps
209: then if dic(0,i) < 0.0
210: then return false;
211: end if;
212: end if;
213: end loop;
214: return true;
215: end Primal_Optimal;
216:
217: function Dual_Optimal ( dic : Matrix; eps : double_float ) return boolean is
218: begin
219: for i in (dic'first(1)+1)..dic'last(1) loop
220: if abs(dic(i,0)) > eps
221: then if dic(i,0) < 0.0
222: then return false;
223: end if;
224: end if;
225: end loop;
226: return true;
227: end Dual_Optimal;
228:
229: function Optimum ( dic : Matrix ) return double_float is
230: begin
231: return dic(dic'first(1),dic'first(2));
232: end Optimum;
233:
234: function Primal_Solution
235: ( dic : Matrix;
236: in_bas,out_bas : Standard_Integer_Vectors.Vector)
237: return Standard_Floating_Vectors.Vector is
238:
239: res : Standard_Floating_Vectors.Vector((dic'first(2)+1)..dic'last(2));
240: n : constant natural := dic'last(2);
241:
242: begin
243: for i in in_bas'range loop
244: if in_bas(i) <= n
245: then res(in_bas(i)) := dic(i,0);
246: end if;
247: end loop;
248: for i in out_bas'range loop
249: if out_bas(i) <= n
250: then res(out_bas(i)) := 0.0;
251: end if;
252: end loop;
253: return res;
254: end Primal_Solution;
255:
256: function Dual_Solution
257: ( dic : Matrix;
258: in_bas,out_bas : Standard_Integer_Vectors.Vector)
259: return Standard_Floating_Vectors.Vector is
260:
261: res : Standard_Floating_Vectors.Vector((dic'first(1)+1)..dic'last(1));
262: n : constant natural := dic'last(2);
263:
264: begin
265: for i in in_bas'range loop
266: if in_bas(i) > n
267: then res(in_bas(i)-n) := dic(i,0);
268: end if;
269: end loop;
270: for i in out_bas'range loop
271: if out_bas(i) > n
272: then res(out_bas(i)-n) := dic(0,i);
273: end if;
274: end loop;
275: return res;
276: end Dual_Solution;
277:
278: end Dictionaries;