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