Annotation of OpenXM_contrib/PHC/Ada/Continuation/standard_root_refiners.adb, Revision 1.1.1.1
1.1 maekawa 1: with integer_io; use integer_io;
2: with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
3: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
4: with Standard_Complex_Numbers_io; use Standard_Complex_Numbers_io;
5: with Standard_Complex_Vectors; use Standard_Complex_Vectors;
6: with Standard_Natural_Vectors;
7: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
8: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
9: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
10: with Standard_Complex_Solutions_io; use Standard_Complex_Solutions_io;
11:
12: package body Standard_Root_Refiners is
13:
14: -- AUXILIARIES :
15:
16: function Is_Real ( sol : Solution; tol : double_float ) return boolean is
17: begin
18: for i in sol.v'range loop
19: if ABS(IMAG_PART(sol.v(i))) > tol
20: then return false;
21: end if;
22: end loop;
23: return true;
24: end Is_Real;
25:
26: function Equal ( s1,s2 : Solution; tol : double_float ) return boolean is
27: begin
28: for i in s1.v'range loop
29: if not Equal(s1.v(i),s2.v(i),tol)
30: then return false;
31: end if;
32: end loop;
33: return true;
34: end Equal;
35:
36: function Complex_Conjugate ( s : Solution ) return Solution is
37:
38: res : Solution(s.n) := s;
39:
40: begin
41: for i in res.v'range loop
42: res.v(i) := Create(REAL_PART(s.v(i)),-IMAG_PART(s.v(i)));
43: end loop;
44: return res;
45: end Complex_Conjugate;
46:
47: function Is_Clustered ( sol : Solution; nb : natural; sols : Solution_List;
48: tol : double_float ) return natural is
49:
50: tmp : Solution_List := sols;
51: cnt : natural := 0;
52:
53: begin
54: while not Is_Null(tmp) loop
55: cnt := cnt + 1;
56: if cnt /= nb
57: then if Equal(sol,Head_Of(tmp).all,tol)
58: then return cnt;
59: end if;
60: end if;
61: tmp := Tail_Of(tmp);
62: end loop;
63: return nb;
64: end Is_Clustered;
65:
66: function Is_Clustered ( sol : Solution; nb : natural; sols : Solution_Array;
67: tol : double_float ) return natural is
68: begin
69: for i in sols'range loop
70: if i /= nb
71: then if Equal(sol,sols(i).all,tol)
72: then return i;
73: end if;
74: end if;
75: end loop;
76: return nb;
77: end Is_Clustered;
78:
79: function Multiplicity ( sol : Solution; sols : Solution_List;
80: tol : double_float ) return natural is
81:
82: tmp : Solution_List := sols;
83: cnt : natural := 0;
84:
85: begin
86: while not Is_Null(tmp) loop
87: if Equal(sol,Head_Of(tmp).all,tol)
88: then cnt := cnt + 1;
89: end if;
90: tmp := Tail_Of(tmp);
91: end loop;
92: return cnt;
93: end Multiplicity;
94:
95: function Multiplicity ( sol : Solution; sols : Solution_Array;
96: tol : double_float ) return natural is
97: cnt : natural := 0;
98:
99: begin
100: for i in sols'range loop
101: if Equal(sol,sols(i).all,tol)
102: then cnt := cnt + 1;
103: end if;
104: end loop;
105: return cnt;
106: end Multiplicity;
107:
108: procedure Write_Bar ( file : in file_type ) is
109: begin
110: for i in 1..75 loop
111: put(file,'=');
112: end loop;
113: new_line(file);
114: end Write_Bar;
115:
116: procedure Write_Info ( file : in file_type; zero : in Solution;
117: i,n,numb : in natural; fail : in boolean ) is
118:
119: -- DESCRIPTION :
120: -- The information concerning the zero is written
121:
122: begin
123: -- Write_Bar(file);
124: put(file,"solution "); put(file,i,1); put(file," : ");
125: put(file," start residual : "); put(file,zero.res,2,3,3); new_line(file);
126: put(file,zero);
127: end Write_Info;
128:
129: procedure Root_Accounting
130: ( file : in file_type; ls : in out Link_to_Solution;
131: nb : in natural; sa : in out Solution_Array;
132: fail : in boolean; tolsing,tolclus : in double_float;
133: nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus : in out natural ) is
134:
135: -- DESCRIPTION :
136: -- This procedure does root accounting of the solution sol, w.r.t. the
137: -- solution list sols. Information will be provided concerning the type
138: -- of solution.
139:
140: begin
141: if fail
142: then put_line(file," no solution ==");
143: nbfail := nbfail + 1;
144: ls.m := 0;
145: else if Is_Real(ls.all,10.0**(-13))
146: then put(file," real ");
147: nbreal := nbreal + 1;
148: else put(file," complex ");
149: nbcomp := nbcomp + 1;
150: end if;
151: if sa(nb).rco < tolsing
152: then declare
153: m : natural := Multiplicity(ls.all,sa,tolclus);
154: begin
155: if m = 1
156: then m := 0;
157: end if;
158: ls.m := m;
159: end;
160: put_line(file,"singular ==");
161: nbsing := nbsing + 1;
162: else declare
163: nb2 : natural := Is_Clustered(ls.all,nb,sa,tolclus);
164: begin
165: if nb2 = nb
166: then put_line(file,"regular ==");
167: nbreg := nbreg + 1;
168: else put(file,"clustered : ");
169: put(file,nb2,1);
170: put_line(file," ==");
171: nbclus := nbclus + 1;
172: end if;
173: ls.m := 1;
174: end;
175: end if;
176: end if;
177: end Root_Accounting;
178:
179: procedure Root_Accounting
180: ( ls : in out Link_to_Solution; nb : in natural;
181: sa : in out Solution_Array; fail : in boolean;
182: tolsing,tolclus : in double_float ) is
183:
184: -- DESCRIPTION :
185: -- This procedure does root accounting of the solution sol, w.r.t. the
186: -- solution list sols. Information will be provided concerning the type
187: -- of solution.
188:
189: begin
190: if fail
191: then ls.m := 0;
192: elsif sa(nb).rco < tolsing
193: then declare
194: m : natural := Multiplicity(ls.all,sa,tolclus);
195: begin
196: if m = 1
197: then ls.m := 0;
198: else ls.m := m;
199: end if;
200: end;
201: else ls.m := 1;
202: end if;
203: end Root_Accounting;
204:
205: procedure Write_Global_Info
206: ( file : in file_type;
207: tot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus : in natural ) is
208:
209: begin
210: Write_Bar(file);
211: put(file,"A list of "); put(file,tot,1);
212: put_line(file," solutions has been refined :");
213: put(file,"Number of regular solutions : "); put(file,nbreg,1);
214: put_line(file,".");
215: put(file,"Number of singular solutions : "); put(file,nbsing,1);
216: put_line(file,".");
217: put(file,"Number of real solutions : "); put(file,nbreal,1);
218: put_line(file,".");
219: put(file,"Number of complex solutions : "); put(file,nbcomp,1);
220: put_line(file,".");
221: put(file,"Number of clustered solutions : "); put(file,nbclus,1);
222: put_line(file,".");
223: put(file,"Number of failures : "); put(file,nbfail,1);
224: put_line(file,".");
225: Write_Bar(file);
226: end Write_Global_Info;
227:
228: -- TARGET ROUTINES :
229:
230: procedure Silent_Newton
231: ( p_eval : in Eval_Poly_Sys; j_eval : in Eval_Jaco_Mat;
232: zero : in out Solution; epsxa,epsfa : in double_float;
233: numit : in out natural; max : in natural;
234: fail : out boolean ) is
235:
236: n : natural := zero.n;
237: jacobian : matrix(1..n,1..n);
238: ipvt : Standard_Natural_Vectors.Vector(1..n);
239: y,deltax : Vector(1..n);
240:
241: begin
242: y := eval(p_eval,zero.v); -- y = f(zero)
243: for i in 1..max loop
244: jacobian := eval(j_eval,zero.v); -- solve jacobian*deltax = -f(zero)
245: lufco(jacobian,n,ipvt,zero.rco);
246: if 1.0 + zero.rco = 1.0
247: then fail := (Sum_Norm(y) > epsfa);
248: return; -- singular Jacobian matrix
249: end if;
250: deltax := -y;
251: lusolve(jacobian,n,ipvt,deltax);
252: Add(zero.v,deltax); -- make the updates
253: y := eval(p_eval,zero.v);
254: zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
255: numit := numit + 1;
256: if ( zero.err < epsxa ) or else ( zero.res < epsfa )
257: -- stopping criteria
258: then fail := false; exit;
259: elsif numit >= max
260: then fail := true; exit;
261: end if;
262: end loop;
263: jacobian := eval(j_eval,zero.v); -- compute condition number
264: lufco(jacobian,n,ipvt,zero.rco);
265: exception
266: when numeric_error | constraint_error => fail := true; return;
267: end Silent_Newton;
268:
269: procedure Silent_Newton
270: ( p_eval : in Standard_Complex_Poly_SysFun.Evaluator;
271: j_eval : in Standard_Complex_Jaco_Matrices.Evaluator;
272: zero : in out Solution; epsxa,epsfa : in double_float;
273: numit : in out natural; max : in natural;
274: fail : out boolean ) is
275:
276: n : natural := zero.n;
277: jacobian : matrix(1..n,1..n);
278: ipvt : Standard_Natural_Vectors.Vector(1..n);
279: y,deltax : Vector(1..n);
280:
281: begin
282: y := p_eval(zero.v); -- y = f(zero)
283: for i in 1..max loop
284: jacobian := j_eval(zero.v); -- solve jacobian*deltax = -f(zero)
285: lufco(jacobian,n,ipvt,zero.rco);
286: if 1.0 + zero.rco = 1.0
287: then fail := (Sum_Norm(y) > epsfa);
288: return; -- singular Jacobian matrix
289: end if;
290: deltax := -y;
291: lusolve(jacobian,n,ipvt,deltax);
292: Add(zero.v,deltax); -- make the updates
293: y := p_eval(zero.v);
294: zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
295: numit := numit + 1;
296: if ( zero.err < epsxa ) or else ( zero.res < epsfa )
297: -- stopping criteria
298: then fail := false; exit;
299: elsif numit >= max
300: then fail := true; exit;
301: end if;
302: end loop;
303: jacobian := j_eval(zero.v); -- compute condition number
304: lufco(jacobian,n,ipvt,zero.rco);
305: exception
306: when numeric_error | constraint_error => fail := true; return;
307: end Silent_Newton;
308:
309: procedure Reporting_Newton
310: ( file : in file_type;
311: p_eval : in Eval_Poly_Sys; j_eval : in Eval_Jaco_Mat;
312: zero : in out Solution; epsxa,epsfa : in double_float;
313: numit : in out natural; max : in natural;
314: fail : out boolean ) is
315:
316: n : natural := zero.n;
317: jacobian : matrix(1..n,1..n);
318: ipvt : Standard_Natural_Vectors.Vector(1..n);
319: y,deltax : Vector(1..n);
320:
321: begin
322: y := Eval(p_eval,zero.v); -- y = f(zero)
323: for i in 1..max loop
324: jacobian := eval(j_eval,zero.v); -- solve jacobian*deltax = -f(zero)
325: lufco(jacobian,n,ipvt,zero.rco);
326: if 1.0 + zero.rco = 1.0 -- singular jacobian matrix
327: then fail := (Sum_Norm(y) > epsfa); -- accuracy not reached yet
328: return;
329: end if;
330: deltax := -y;
331: lusolve(jacobian,n,ipvt,deltax);
332: Add(zero.v,deltax); -- make the updates
333: y := eval(p_eval,zero.v);
334: zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
335: numit := numit + 1;
336: put(file,"Step "); put(file,numit,4); new_line(file); -- output
337: put(file," |errxa| : "); put(file,zero.err); new_line(file);
338: put(file," |errfa| : "); put(file,zero.res); new_line(file);
339: if ( zero.err < epsxa ) or else ( zero.res < epsfa )
340: -- stopping criteria
341: then fail := false; exit;
342: elsif numit >= max
343: then fail := true; exit;
344: end if;
345: end loop;
346: jacobian := eval(j_eval,zero.v); -- compute condition number
347: lufco(jacobian,n,ipvt,zero.rco);
348: exception
349: when numeric_error | constraint_error => fail := true; return;
350: end Reporting_Newton;
351:
352: procedure Reporting_Newton
353: ( file : in file_type;
354: p_eval : in Standard_Complex_Poly_SysFun.Evaluator;
355: j_eval : in Standard_Complex_Jaco_Matrices.Evaluator;
356: zero : in out Solution; epsxa,epsfa : in double_float;
357: numit : in out natural; max : in natural;
358: fail : out boolean ) is
359:
360: n : natural := zero.n;
361: jacobian : matrix(1..n,1..n);
362: ipvt : Standard_Natural_Vectors.Vector(1..n);
363: y,deltax : Vector(1..n);
364:
365: begin
366: y := p_eval(zero.v); -- y = f(zero)
367: for i in 1..max loop
368: jacobian := j_eval(zero.v); -- solve jacobian*deltax = -f(zero)
369: lufco(jacobian,n,ipvt,zero.rco);
370: if 1.0 + zero.rco = 1.0 -- singular jacobian matrix
371: then fail := (Sum_Norm(y) > epsfa); -- accuracy not reached yet
372: return;
373: end if;
374: deltax := -y;
375: lusolve(jacobian,n,ipvt,deltax);
376: Add(zero.v,deltax); -- make the updates
377: y := p_eval(zero.v);
378: zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
379: numit := numit + 1;
380: put(file,"Step "); put(file,numit,4); new_line(file); -- output
381: put(file," |errxa| : "); put(file,zero.err); new_line(file);
382: put(file," |errfa| : "); put(file,zero.res); new_line(file);
383: if ( zero.err < epsxa ) or else ( zero.res < epsfa )
384: -- stopping criteria
385: then fail := false; exit;
386: elsif numit >= max
387: then fail := true; exit;
388: end if;
389: end loop;
390: jacobian := j_eval(zero.v); -- compute condition number
391: lufco(jacobian,n,ipvt,zero.rco);
392: exception
393: when numeric_error | constraint_error => fail := true; return;
394: end Reporting_Newton;
395:
396: procedure Silent_Root_Refiner
397: ( p : in Poly_Sys; sols : in out Solution_List;
398: epsxa,epsfa,tolsing : in double_float;
399: numit : in out natural; max : in natural ) is
400:
401: n : natural := p'length;
402: p_eval : Eval_Poly_Sys(1..n) := Create(p);
403: jac : Jaco_Mat(1..n,1..n) := Create(p);
404: jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
405: numb : natural;
406: fail : boolean;
407: sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
408:
409: begin
410: for i in sa'range loop
411: numb := 0;
412: sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
413: if sa(i).res < 1.0
414: then Silent_Newton(p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
415: else fail := true;
416: end if;
417: Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
418: numit := numit + numb;
419: end loop;
420: clear(p_eval); clear(jac); clear(jac_eval);
421: Deep_Clear(sols); sols := Create(sa); Clear(sa);
422: end Silent_Root_Refiner;
423:
424: procedure Silent_Root_Refiner
425: ( p : in Standard_Complex_Poly_SysFun.Evaluator;
426: j : in Standard_Complex_Jaco_Matrices.Evaluator;
427: sols : in out Solution_List;
428: epsxa,epsfa,tolsing : in double_float;
429: numit : in out natural; max : in natural ) is
430:
431: n : constant natural := Head_Of(sols).n;
432: numb : natural;
433: fail : boolean;
434: sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
435:
436: begin
437: for i in sa'range loop
438: numb := 0;
439: sa(i).res := Sum_Norm(p(sa(i).v));
440: if sa(i).res < 1.0
441: then Silent_Newton(p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
442: else fail := true;
443: end if;
444: Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
445: numit := numit + numb;
446: end loop;
447: Deep_Clear(sols); sols := Create(sa); Clear(sa);
448: end Silent_Root_Refiner;
449:
450: procedure Silent_Root_Refiner
451: ( p : in Poly_Sys; sols,refsols : in out Solution_List;
452: epsxa,epsfa,tolsing : in double_float;
453: numit : in out natural; max : in natural ) is
454:
455: n : natural := p'length;
456: p_eval : Eval_Poly_Sys(1..n) := Create(p);
457: jac : Jaco_Mat(1..n,1..n) := Create(p);
458: jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
459: numb : natural;
460: fail : boolean;
461: sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
462: refsols_last : Solution_List;
463:
464: begin
465: for i in sa'range loop
466: numb := 0;
467: sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
468: if sa(i).res < 1.0
469: then Silent_Newton(p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
470: if not fail
471: then Append(refsols,refsols_last,sa(i).all);
472: end if;
473: else fail := true;
474: end if;
475: Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
476: numit := numit + numb;
477: end loop;
478: clear(p_eval); clear(jac); clear(jac_eval);
479: Deep_Clear(sols); sols := Create(sa); Clear(sa);
480: end Silent_Root_Refiner;
481:
482: procedure Silent_Root_Refiner
483: ( p : in Standard_Complex_Poly_SysFun.Evaluator;
484: j : in Standard_Complex_Jaco_Matrices.Evaluator;
485: sols,refsols : in out Solution_List;
486: epsxa,epsfa,tolsing : in double_float;
487: numit : in out natural; max : in natural ) is
488:
489: n : constant natural := Head_Of(sols).n;
490: numb : natural;
491: fail : boolean;
492: sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
493: refsols_last : Solution_List;
494:
495: begin
496: for i in sa'range loop
497: numb := 0;
498: sa(i).res := Sum_Norm(p(sa(i).v));
499: if sa(i).res < 1.0
500: then Silent_Newton(p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
501: if not fail
502: then Append(refsols,refsols_last,sa(i).all);
503: end if;
504: else fail := true;
505: end if;
506: Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
507: numit := numit + numb;
508: end loop;
509: Deep_Clear(sols); sols := Create(sa); Clear(sa);
510: end Silent_Root_Refiner;
511:
512: procedure Reporting_Root_Refiner
513: ( file : in file_type;
514: p : in Poly_Sys; sols : in out Solution_List;
515: epsxa,epsfa,tolsing : in double_float;
516: numit : in out natural; max : in natural;
517: wout : in boolean ) is
518:
519: n : natural := p'length;
520: p_eval : Eval_Poly_Sys(1..n) := Create(p);
521: jac : Jaco_Mat(1..n,1..n) := Create(p);
522: jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
523: numb : natural;
524: nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
525: nbtot : constant natural := Length_Of(sols);
526: fail : boolean;
527: sa : Solution_Array(1..nbtot) := Create(sols);
528:
529: begin
530: new_line(file);
531: put_line(file,"THE SOLUTIONS :"); new_line(file);
532: put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
533: Write_Bar(file);
534: for i in sa'range loop
535: numb := 0;
536: sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
537: if sa(i).res < 1.0
538: then
539: if wout
540: then Reporting_Newton
541: (file,p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
542: else Silent_Newton
543: (p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
544: end if;
545: else
546: fail := true;
547: end if;
548: Write_Info(file,sa(i).all,i,n,numb,fail);
549: Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
550: nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
551: numit := numit + numb;
552: end loop;
553: Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
554: clear(p_eval); clear(jac); clear(jac_eval);
555: Deep_Clear(sols); sols := Create(sa); Clear(sa);
556: end Reporting_Root_Refiner;
557:
558: procedure Reporting_Root_Refiner
559: ( file : in file_type;
560: p : in Standard_Complex_Poly_SysFun.Evaluator;
561: j : in Standard_Complex_Jaco_Matrices.Evaluator;
562: sols : in out Solution_List;
563: epsxa,epsfa,tolsing : in double_float;
564: numit : in out natural; max : in natural;
565: wout : in boolean ) is
566:
567: n : constant natural := Head_Of(sols).n;
568: numb : natural;
569: nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
570: nbtot : constant natural := Length_Of(sols);
571: fail : boolean;
572: sa : Solution_Array(1..nbtot) := Create(sols);
573:
574: begin
575: new_line(file);
576: put_line(file,"THE SOLUTIONS :"); new_line(file);
577: put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
578: Write_Bar(file);
579: for i in sa'range loop
580: numb := 0;
581: sa(i).res := Sum_Norm(p(sa(i).v));
582: if sa(i).res < 1.0
583: then
584: if wout
585: then Reporting_Newton(file,p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
586: else Silent_Newton(p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
587: end if;
588: else
589: fail := true;
590: end if;
591: Write_Info(file,sa(i).all,i,n,numb,fail);
592: Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
593: nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
594: numit := numit + numb;
595: end loop;
596: Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
597: Deep_Clear(sols); sols := Create(sa); Clear(sa);
598: end Reporting_Root_Refiner;
599:
600: procedure Reporting_Root_Refiner
601: ( file : in file_type;
602: p : in Poly_Sys; sols,refsols : in out Solution_List;
603: epsxa,epsfa,tolsing : in double_float;
604: numit : in out natural; max : in natural;
605: wout : in boolean ) is
606:
607: n : natural := p'length;
608: p_eval : Eval_Poly_Sys(1..n) := Create(p);
609: jac : Jaco_Mat(1..n,1..n) := Create(p);
610: jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
611: numb : natural;
612: nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
613: nbtot : constant natural := Length_Of(sols);
614: fail : boolean;
615: sa : Solution_Array(1..nbtot) := Create(sols);
616: refsols_last : Solution_List;
617:
618: begin
619: new_line(file);
620: put_line(file,"THE SOLUTIONS :"); new_line(file);
621: put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
622: Write_Bar(file);
623: for i in sa'range loop
624: numb := 0;
625: sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
626: if sa(i).res < 1.0
627: then
628: if wout
629: then Reporting_Newton
630: (file,p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
631: else Silent_Newton
632: (p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
633: end if;
634: if not fail
635: then Append(refsols,refsols_last,sa(i).all);
636: end if;
637: else
638: fail := true;
639: end if;
640: Write_Info(file,sa(i).all,i,n,numb,fail);
641: Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
642: nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
643: numit := numit + numb;
644: end loop;
645: Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
646: clear(p_eval); clear(jac); clear(jac_eval);
647: Deep_Clear(sols); sols := Create(sa); Clear(sa);
648: end Reporting_Root_Refiner;
649:
650: procedure Reporting_Root_Refiner
651: ( file : in file_type;
652: p : in Standard_Complex_Poly_SysFun.Evaluator;
653: j : in Standard_Complex_Jaco_Matrices.Evaluator;
654: sols,refsols : in out Solution_List;
655: epsxa,epsfa,tolsing : in double_float;
656: numit : in out natural; max : in natural;
657: wout : in boolean ) is
658:
659: n : constant natural := Head_Of(sols).n;
660: numb : natural;
661: nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
662: nbtot : constant natural := Length_Of(sols);
663: fail : boolean;
664: sa : Solution_Array(1..nbtot) := Create(sols);
665: refsols_last : Solution_List;
666:
667: begin
668: new_line(file);
669: put_line(file,"THE SOLUTIONS :"); new_line(file);
670: put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
671: Write_Bar(file);
672: for i in sa'range loop
673: numb := 0;
674: sa(i).res := Sum_Norm(p(sa(i).v));
675: if sa(i).res < 1.0
676: then
677: if wout
678: then Reporting_Newton(file,p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
679: else Silent_Newton(p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
680: end if;
681: if not fail
682: then Append(refsols,refsols_last,sa(i).all);
683: end if;
684: else
685: fail := true;
686: end if;
687: Write_Info(file,sa(i).all,i,n,numb,fail);
688: Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
689: nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
690: numit := numit + numb;
691: end loop;
692: Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
693: Deep_Clear(sols); sols := Create(sa); Clear(sa);
694: end Reporting_Root_Refiner;
695:
696: end Standard_Root_Refiners;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to standard_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