![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to binomial_system_solvers.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift|
Return to binomial_system_solvers.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift |
1.1 maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
2: with Standard_Complex_Numbers_Polar; use Standard_Complex_Numbers_Polar;
3: with Standard_Integer_Vectors;
4: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
5: with Transforming_Solutions; use Transforming_Solutions;
6:
7: package body Binomial_System_Solvers is
8:
9: -- AN INTERMEDIATE ROUTINE : --
10:
11: procedure Factorize ( v : in out VecVec; n : in natural;
12: t : in out Transfo ) is
13:
14: -- DESCRIPTION :
15: -- This routines factorizes the binomial system defined by the
16: -- degrees in the Vector v;
17:
18: -- ON ENTRY :
19: -- v defines the degrees of binomial system
20: -- n the number of unknowns to be eliminated
21:
22: -- ON RETURN :
23: -- v the factorized binomial system
24: -- t the transformations used to factorize
25:
26: tt : Transfo;
27: pivot : integer;
28: tmp : Standard_Integer_Vectors.Link_to_Vector;
29:
30: begin
31: t := Id(v'last);
32: for i in 1..n loop
33: pivot := 0; -- search for pivot
34: for j in i..v'last loop
35: if v(j)(i) /= 0
36: then pivot := j;
37: exit;
38: end if;
39: end loop;
40: if pivot /= 0 -- interchange if necessary
41: then if pivot /= i
42: then tmp := v(i);
43: v(i) := v(pivot);
44: v(pivot) := tmp;
45: end if;
46: tt := Rotate(v(i),i);
47: Apply(tt,v(i));
48: for j in (i+1)..v'last loop
49: Apply(tt,v(j));
50: v(j).all(i) := 0;
51: end loop;
52: Mult1(t,tt);
53: Clear(tt);
54: end if;
55: end loop;
56: end Factorize;
57:
58: -- AUXILIARY ROUTINES :
59:
60: procedure Update ( sols1,sols2 : in out Solution_List ) is
61:
62: tmp : Solution_List := sols2;
63: ls : Link_to_Solution;
64:
65: begin
66: while not Is_Null(tmp) loop
67: ls := Head_Of(tmp);
68: declare
69: nls : Link_to_Solution := new Solution'(ls.all);
70: begin
71: Construct(nls,sols1);
72: end;
73: tmp := Tail_Of(tmp);
74: end loop;
75: Clear(sols2);
76: end Update;
77:
78: procedure Solve ( v : in VecVec; cv : in Vector;
79: i,n : in natural; sols : in out Solution_List;
80: acc : in out Vector ) is
81:
82: t : Transfo;
83: pivot,d : integer;
84: a : Complex_Number;
85: c : Vector(cv'range) := cv;
86: rv : Vector(cv'range);
87: tmp : Standard_Integer_Vectors.Link_to_Vector;
88: nv : VecVec(v'range);
89: temp : array(v'range) of integer;
90: newsols : Solution_List;
91:
92: begin
93: for j in i..v'last loop
94: nv(j) := new Standard_Integer_Vectors.Vector'(v(j).all);
95: end loop;
96: pivot := 0; -- search for pivot
97: for j in i..v'last loop
98: if nv(j)(i) /= 0
99: then pivot := j;
100: exit;
101: end if;
102: end loop;
103: if pivot /= 0 -- interchange if necessary
104: then if pivot /= i
105: then tmp := nv(i);
106: nv(i) := nv(pivot);
107: nv(pivot) := tmp;
108: a := c(i);
109: c(i) := c(pivot);
110: c(pivot) := a;
111: end if;
112: t := Rotate(nv(i),i);
113: Apply(t,nv(i));
114: for j in (i+1)..nv'last loop
115: Apply(t,nv(j));
116: temp(j) := nv(j)(i);
117: nv(j)(i) := 0;
118: end loop;
119: d := nv(i)(i); -- compute the solutions
120: if d < 0
121: then d := -d;
122: rv(i) := Create(1.0)/c(i);
123: else rv(i) := c(i);
124: end if;
125: for j in 1..d loop
126: a := Root(rv(i),d,j);
127: acc(i) := a;
128: for k in (i+1)..nv'last loop
129: rv(k) := c(k)*a**(-temp(k));
130: end loop;
131: if i < n
132: then Solve(nv,rv,i+1,n,newsols,acc);
133: else declare -- i = n
134: ls : Link_to_Solution;
135: s : Solution(n);
136: begin
137: s.t := Create(0.0);
138: s.m := 1;
139: s.v := acc;
140: ls := new Solution'(s);
141: Construct(ls,newsols);
142: end;
143: end if;
144: Transform(t,newsols);
145: Update(sols,newsols);
146: end loop;
147: Clear(t);
148: end if;
149: for j in i..nv'last loop
150: Standard_Integer_Vectors.Clear(nv(j));
151: end loop;
152: end Solve;
153:
154: -- THE SOLVER :
155:
156: procedure Solve ( v : in VecVec; cv : in Vector;
157: n : in natural; sols : in out Solution_List ) is
158:
159: wv : VecVec(v'range); -- workspace
160: acc : Vector(cv'range); -- accumulator
161:
162: begin
163: for i in v'range loop
164: wv(i) := new Standard_Integer_Vectors.Vector'(v(i).all);
165: end loop;
166: Solve(wv,cv,v'first,v'last,sols,acc);
167: Clear(wv);
168: end Solve;
169:
170: -- COMPUTING THE RESIDUALS :
171:
172: procedure Residuals ( v : in VecVec; cv : in Vector;
173: n : in natural; sols : in Solution_List;
174: res : out Vector ) is
175: nb : natural;
176: x : Vector(cv'range);
177: tmp : Solution_List;
178:
179: function Eval ( v : VecVec; cv,x : Vector ) return Vector is
180:
181: -- DESCRIPTION :
182: -- Computes the value of the vector x in the system defined by
183: -- v and cv.
184:
185: y : Standard_Complex_Vectors.Vector(x'range);
186:
187: begin
188: for i in x'range loop
189: y(i) := Create(1.0);
190: for j in v(i)'range loop
191: y(i) := y(i)*x(j)**v(i)(j);
192: end loop;
193: y(i) := y(i) - cv(i);
194: end loop;
195: return y;
196: end Eval;
197:
198: begin
199: tmp := sols;
200: nb := 1;
201: while not Is_Null(tmp) loop
202: x := Head_Of(tmp).v;
203: res(nb) := Create(Max_Norm(Eval(v,cv,x)));
204: tmp := Tail_Of(tmp);
205: nb := nb + 1;
206: end loop;
207: end Residuals;
208:
209: end Binomial_System_Solvers;