Annotation of OpenXM_contrib/PHC/Ada/Continuation/standard_root_refiners.adb, Revision 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