Annotation of OpenXM_contrib/PHC/Ada/Continuation/multprec_root_refiners.adb, Revision 1.1.1.1
1.1 maekawa 1: with integer_io; use integer_io;
2: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
3: with Multprec_Floating_Numbers_io; use Multprec_Floating_Numbers_io;
4: with Multprec_Complex_Numbers; use Multprec_Complex_Numbers;
5: with Multprec_Complex_Numbers_io; use Multprec_Complex_Numbers_io;
6: with Multprec_Complex_Vectors; use Multprec_Complex_Vectors;
7: with Standard_Natural_Vectors;
8: with Multprec_Complex_Matrices; use Multprec_Complex_Matrices;
9: with Multprec_Complex_Linear_Solvers; use Multprec_Complex_Linear_Solvers;
10: with Multprec_Complex_Norms_Equals; use Multprec_Complex_Norms_Equals;
11: with Multprec_Complex_Solutions_io; use Multprec_Complex_Solutions_io;
12:
13: package body Multprec_Root_Refiners is
14:
15: -- AUXILIARIES :
16:
17: function Is_Real ( sol : Solution; tol : Floating_Number ) return boolean is
18:
19: res : boolean := true;
20:
21: begin
22: for i in sol.v'range loop
23: declare
24: abstmp : Floating_Number := AbsVal(IMAG_PART(sol.v(i)));
25: begin
26: if abstmp > tol
27: then res := false;
28: end if;
29: Clear(abstmp);
30: end;
31: exit when not res;
32: end loop;
33: return res;
34: end Is_Real;
35:
36: function Equal ( s1,s2 : Solution; tol : Floating_Number ) return boolean is
37: begin
38: for i in s1.v'range loop
39: if not Equal(s1.v(i),s2.v(i),tol)
40: then return false;
41: end if;
42: end loop;
43: return true;
44: end Equal;
45:
46: function Complex_Conjugate ( s : Solution ) return Solution is
47:
48: res : Solution(s.n) := s;
49:
50: begin
51: for i in res.v'range loop
52: res.v(i) := Create(REAL_PART(s.v(i)),-IMAG_PART(s.v(i)));
53: end loop;
54: return res;
55: end Complex_Conjugate;
56:
57: function Is_Clustered ( sol : Solution; nb : natural; sols : Solution_List;
58: tol : Floating_Number ) return natural is
59:
60: tmp : Solution_List := sols;
61: cnt : natural := 0;
62:
63: begin
64: while not Is_Null(tmp) loop
65: cnt := cnt + 1;
66: if cnt /= nb
67: then if Equal(sol,Head_Of(tmp).all,tol)
68: then return cnt;
69: end if;
70: end if;
71: tmp := Tail_Of(tmp);
72: end loop;
73: return nb;
74: end Is_Clustered;
75:
76: function Is_Clustered ( sol : Solution; nb : natural; sols : Solution_Array;
77: tol : Floating_Number ) return natural is
78: begin
79: for i in sols'range loop
80: if i /= nb
81: then if Equal(sol,sols(i).all,tol)
82: then return i;
83: end if;
84: end if;
85: end loop;
86: return nb;
87: end Is_Clustered;
88:
89: function Multiplicity ( sol : Solution; sols : Solution_List;
90: tol : Floating_Number ) return natural is
91:
92: tmp : Solution_List := sols;
93: cnt : natural := 0;
94:
95: begin
96: while not Is_Null(tmp) loop
97: if Equal(sol,Head_Of(tmp).all,tol)
98: then cnt := cnt + 1;
99: end if;
100: tmp := Tail_Of(tmp);
101: end loop;
102: return cnt;
103: end Multiplicity;
104:
105: function Multiplicity ( sol : Solution; sols : Solution_Array;
106: tol : Floating_Number ) return natural is
107: cnt : natural := 0;
108:
109: begin
110: for i in sols'range loop
111: if Equal(sol,sols(i).all,tol)
112: then cnt := cnt + 1;
113: end if;
114: end loop;
115: return cnt;
116: end Multiplicity;
117:
118: procedure Write_Bar ( file : in file_type ) is
119: begin
120: for i in 1..75 loop
121: put(file,'=');
122: end loop;
123: new_line(file);
124: end Write_Bar;
125:
126: procedure Write_Info ( file : in file_type; zero : in Solution;
127: i,n,numb : in natural; fail : in boolean ) is
128:
129: -- DESCRIPTION :
130: -- The information concerning the zero is written
131:
132: begin
133: -- Write_Bar(file);
134: put(file,"solution "); put(file,i,1); put(file," : ");
135: put(file," start residual : "); put(file,zero.res,3); new_line(file);
136: put(file,zero);
137: end Write_Info;
138:
139: procedure Root_Accounting
140: ( file : in file_type; ls : in out Link_to_Solution;
141: nb : in natural; sa : in out Solution_Array;
142: fail : in boolean; tolsing,tolclus : in Floating_Number;
143: nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus : in out natural ) is
144:
145: -- DESCRIPTION :
146: -- This procedure does root accounting of the solution sol, w.r.t. the
147: -- solution list sols. Information will be provided concerning the type
148: -- of solution.
149:
150: tolreal : Floating_Number := Create(10.0**(-13));
151:
152: begin
153: if fail
154: then put_line(file," no solution ==");
155: nbfail := nbfail + 1;
156: ls.m := 0;
157: else if Is_Real(ls.all,tolreal)
158: then put(file," real ");
159: nbreal := nbreal + 1;
160: else put(file," complex ");
161: nbcomp := nbcomp + 1;
162: end if;
163: if sa(nb).rco < tolsing
164: then declare
165: m : natural := Multiplicity(ls.all,sa,tolclus);
166: begin
167: if m = 1
168: then m := 0;
169: end if;
170: ls.m := m;
171: end;
172: put_line(file,"singular ==");
173: nbsing := nbsing + 1;
174: else declare
175: nb2 : natural := Is_Clustered(ls.all,nb,sa,tolclus);
176: begin
177: if nb2 = nb
178: then put_line(file,"regular ==");
179: nbreg := nbreg + 1;
180: else put(file,"clustered : ");
181: put(file,nb2,1);
182: put_line(file," ==");
183: nbclus := nbclus + 1;
184: end if;
185: ls.m := 1;
186: end;
187: end if;
188: end if;
189: Clear(tolreal);
190: end Root_Accounting;
191:
192: procedure Root_Accounting
193: ( ls : in out Link_to_Solution; nb : in natural;
194: sa : in out Solution_Array; fail : in boolean;
195: tolsing,tolclus : in Floating_Number ) is
196:
197: -- DESCRIPTION :
198: -- This procedure does root accounting of the solution sol, w.r.t. the
199: -- solution list sols. Information will be provided concerning the type
200: -- of solution.
201:
202: begin
203: if fail
204: then ls.m := 0;
205: elsif sa(nb).rco < tolsing
206: then declare
207: m : natural := Multiplicity(ls.all,sa,tolclus);
208: begin
209: if m = 1
210: then ls.m := 0;
211: else ls.m := m;
212: end if;
213: end;
214: else ls.m := 1;
215: end if;
216: end Root_Accounting;
217:
218: procedure Write_Global_Info
219: ( file : in file_type;
220: tot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus : in natural ) is
221:
222: begin
223: Write_Bar(file);
224: put(file,"A list of "); put(file,tot,1);
225: put_line(file," solutions has been refined :");
226: put(file,"Number of regular solutions : "); put(file,nbreg,1);
227: put_line(file,".");
228: put(file,"Number of singular solutions : "); put(file,nbsing,1);
229: put_line(file,".");
230: put(file,"Number of real solutions : "); put(file,nbreal,1);
231: put_line(file,".");
232: put(file,"Number of complex solutions : "); put(file,nbcomp,1);
233: put_line(file,".");
234: put(file,"Number of clustered solutions : "); put(file,nbclus,1);
235: put_line(file,".");
236: put(file,"Number of failures : "); put(file,nbfail,1);
237: put_line(file,".");
238: Write_Bar(file);
239: end Write_Global_Info;
240:
241: -- TARGET ROUTINES :
242:
243: procedure Silent_Newton
244: ( p_eval : in Eval_Poly_Sys; j_eval : in Eval_Jaco_Mat;
245: zero : in out Solution; epsxa,epsfa : in Floating_Number;
246: numit : in out natural; max : in natural;
247: fail : out boolean ) is
248:
249: n : natural := p_eval'length;
250: jacobian : matrix(1..n,1..n);
251: ipvt : Standard_Natural_Vectors.Vector(1..n);
252: y,deltax : Vector(1..n);
253:
254: begin
255: y := eval(p_eval,zero.v); -- y = f(zero)
256: for i in 1..max loop
257: jacobian := Eval(j_eval,zero.v); -- solve jacobian*deltax = -f(zero)
258: lufco(jacobian,n,ipvt,zero.rco);
259: if Equal(zero.rco,0.0)
260: then fail := (Sum_Norm(y) > epsfa);
261: return; -- singular Jacobian matrix
262: end if;
263: deltax := -y; Clear(y);
264: lusolve(jacobian,n,ipvt,deltax);
265: Add(zero.v,deltax); -- make the updates
266: y := eval(p_eval,zero.v);
267: Clear(zero.err); Clear(zero.res);
268: zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
269: numit := numit + 1;
270: if ( zero.err < epsxa ) or else ( zero.res < epsfa )
271: -- stopping criteria
272: then fail := false; exit;
273: elsif numit >= max
274: then fail := true; exit;
275: end if;
276: Clear(deltax); Clear(jacobian);
277: end loop;
278: jacobian := eval(j_eval,zero.v); -- compute condition number
279: lufco(jacobian,n,ipvt,zero.rco);
280: Clear(jacobian);
281: exception
282: when numeric_error | constraint_error => fail := true; return;
283: end Silent_Newton;
284:
285: procedure Reporting_Newton
286: ( file : in file_type;
287: p_eval : in Eval_Poly_Sys; j_eval : in Eval_Jaco_Mat;
288: zero : in out Solution; epsxa,epsfa : in Floating_Number;
289: numit : in out natural; max : in natural;
290: fail : out boolean ) is
291:
292: n : natural := p_eval'length;
293: jacobian : matrix(1..n,1..n);
294: ipvt : Standard_Natural_Vectors.Vector(1..n);
295: y,deltax : Vector(1..n);
296:
297: begin
298: y := Eval(p_eval,zero.v); -- y = f(zero)
299: for i in 1..max loop
300: jacobian := eval(j_eval,zero.v); -- solve jacobian*deltax = -f(zero)
301: lufco(jacobian,n,ipvt,zero.rco);
302: if Equal(zero.rco,0.0) -- singular jacobian matrix
303: then fail := (Sum_Norm(y) > epsfa); -- accuracy not reached yet
304: return;
305: end if;
306: deltax := -y; Clear(y);
307: lusolve(jacobian,n,ipvt,deltax);
308: Add(zero.v,deltax); -- make the updates
309: y := eval(p_eval,zero.v);
310: Clear(zero.err); Clear(zero.res);
311: zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
312: numit := numit + 1;
313: put(file,"Step "); put(file,numit,4); new_line(file); -- output
314: put(file," |errxa| : "); put(file,zero.err); new_line(file);
315: put(file," |errfa| : "); put(file,zero.res); new_line(file);
316: if ( zero.err < epsxa ) or else ( zero.res < epsfa )
317: -- stopping criteria
318: then fail := false; exit;
319: elsif numit >= max
320: then fail := true; exit;
321: end if;
322: Clear(deltax); Clear(jacobian);
323: end loop;
324: jacobian := eval(j_eval,zero.v); -- compute condition number
325: lufco(jacobian,n,ipvt,zero.rco);
326: Clear(jacobian);
327: exception
328: when numeric_error | constraint_error => fail := true; return;
329: end Reporting_Newton;
330:
331: procedure Silent_Root_Refiner
332: ( p : in Poly_Sys; sols : in out Solution_List;
333: epsxa,epsfa,tolsing : in Floating_Number;
334: numit : in out natural; max : in natural ) is
335:
336: n : natural := p'length;
337: p_eval : Eval_Poly_Sys(1..n) := Create(p);
338: jac : Jaco_Mat(1..n,1..n) := Create(p);
339: jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
340: numb : natural;
341: fail : boolean;
342: sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
343: eval_acc : Vector(1..n);
344:
345: begin
346: for i in sa'range loop
347: numb := 0;
348: Clear(sa(i).res);
349: eval_acc := Eval(p_eval,sa(i).v);
350: sa(i).res := Sum_Norm(eval_acc);
351: Clear(eval_acc);
352: if sa(i).res < 1.0
353: then Silent_Newton(p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
354: else fail := true;
355: end if;
356: Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
357: numit := numit + numb;
358: end loop;
359: clear(p_eval); clear(jac); clear(jac_eval);
360: Shallow_Clear(sols); sols := Create(sa); -- Shallow_Clear(sa);
361: end Silent_Root_Refiner;
362:
363: procedure Silent_Root_Refiner
364: ( p : in Poly_Sys; sols,refsols : in out Solution_List;
365: epsxa,epsfa,tolsing : in Floating_Number;
366: numit : in out natural; max : in natural ) is
367:
368: n : natural := p'length;
369: p_eval : Eval_Poly_Sys(1..n) := Create(p);
370: jac : Jaco_Mat(1..n,1..n) := Create(p);
371: jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
372: numb : natural;
373: fail : boolean;
374: sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
375: refsols_last : Solution_List;
376: eval_acc : Vector(1..n);
377:
378: begin
379: for i in sa'range loop
380: numb := 0;
381: Clear(sa(i).res);
382: eval_acc := Eval(p_eval,sa(i).v);
383: sa(i).res := Sum_Norm(eval_acc);
384: Clear(eval_acc);
385: if sa(i).res < 1.0
386: then Silent_Newton(p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
387: if not fail
388: then Append(refsols,refsols_last,sa(i).all);
389: end if;
390: else fail := true;
391: end if;
392: Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
393: numit := numit + numb;
394: end loop;
395: clear(p_eval); clear(jac); clear(jac_eval);
396: Shallow_Clear(sols); sols := Create(sa); -- Shallow_Clear(sa);
397: end Silent_Root_Refiner;
398:
399: procedure Reporting_Root_Refiner
400: ( file : in file_type;
401: p : in Poly_Sys; sols : in out Solution_List;
402: epsxa,epsfa,tolsing : in Floating_Number;
403: numit : in out natural; max : in natural;
404: wout : in boolean ) is
405:
406: n : natural := p'length;
407: p_eval : Eval_Poly_Sys(1..n) := Create(p);
408: jac : Jaco_Mat(1..n,1..n) := Create(p);
409: jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
410: numb : natural;
411: nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
412: nbtot : constant natural := Length_Of(sols);
413: fail : boolean;
414: sa : Solution_Array(1..nbtot) := Create(sols);
415: eval_acc : Vector(1..n);
416:
417: begin
418: new_line(file);
419: put_line(file,"THE SOLUTIONS :"); new_line(file);
420: put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
421: Write_Bar(file);
422: for i in sa'range loop
423: numb := 0;
424: Clear(sa(i).res);
425: eval_acc := Eval(p_eval,sa(i).v);
426: sa(i).res := Sum_Norm(eval_acc);
427: Clear(eval_acc);
428: if sa(i).res < 1.0
429: then
430: if wout
431: then Reporting_Newton
432: (file,p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
433: else Silent_Newton
434: (p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
435: end if;
436: else
437: fail := true;
438: end if;
439: Write_Info(file,sa(i).all,i,n,numb,fail);
440: Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
441: nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
442: numit := numit + numb;
443: end loop;
444: Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
445: clear(p_eval); clear(jac); clear(jac_eval);
446: Shallow_Clear(sols); sols := Create(sa); -- Shallow_Clear(sa);
447: end Reporting_Root_Refiner;
448:
449: procedure Reporting_Root_Refiner
450: ( file : in file_type;
451: p : in Poly_Sys; sols,refsols : in out Solution_List;
452: epsxa,epsfa,tolsing : in Floating_Number;
453: numit : in out natural; max : in natural;
454: wout : in boolean ) is
455:
456: n : natural := p'length;
457: p_eval : Eval_Poly_Sys(1..n) := Create(p);
458: jac : Jaco_Mat(1..n,1..n) := Create(p);
459: jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
460: numb : natural;
461: nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
462: nbtot : constant natural := Length_Of(sols);
463: fail : boolean;
464: sa : Solution_Array(1..nbtot) := Create(sols);
465: refsols_last : Solution_List;
466: eval_acc : Vector(1..n);
467:
468: begin
469: new_line(file);
470: put_line(file,"THE SOLUTIONS :"); new_line(file);
471: put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
472: Write_Bar(file);
473: for i in sa'range loop
474: numb := 0;
475: Clear(sa(i).res);
476: eval_acc := Eval(p_eval,sa(i).v);
477: sa(i).res := Sum_Norm(eval_acc);
478: Clear(eval_acc);
479: if sa(i).res < 1.0
480: then
481: if wout
482: then Reporting_Newton
483: (file,p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
484: else Silent_Newton
485: (p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
486: end if;
487: if not fail
488: then Append(refsols,refsols_last,sa(i).all);
489: end if;
490: else
491: fail := true;
492: end if;
493: Write_Info(file,sa(i).all,i,n,numb,fail);
494: Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
495: nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
496: numit := numit + numb;
497: end loop;
498: Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
499: clear(p_eval); clear(jac); clear(jac_eval);
500: Shallow_Clear(sols); sols := Create(sa); -- Shallow_Clear(sa);
501: end Reporting_Root_Refiner;
502:
503: end Multprec_Root_Refiners;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to multprec_root_refiners.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Continuation