Annotation of OpenXM_contrib/PHC/Ada/Continuation/path_trackers.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_io; use Standard_Complex_Numbers_io;
! 4: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
! 5: with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
! 6: with Standard_Floating_VecVecs; use Standard_Floating_VecVecs;
! 7: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
! 8: with Predictors,Correctors; use Predictors,Correctors;
! 9: with Dispatch_Predictors; use Dispatch_Predictors;
! 10: with Continuation_Parameters; use Continuation_Parameters;
! 11: with Process_io; use Process_io;
! 12: with Directions_of_Solution_Paths; use Directions_of_Solution_Paths;
! 13:
! 14: package body Path_Trackers is
! 15:
! 16: -- AUXILIAIRIES :
! 17:
! 18: function At_End ( t,target : Complex_Number; distance,tol : double_float )
! 19: return boolean is
! 20:
! 21: -- DESCRIPTION :
! 22: -- Decides whether at end of continuation stage.
! 23:
! 24: begin
! 25: if distance < tol
! 26: then return Equal(t,target,tol);
! 27: else if AbsVal(t-target) <= distance
! 28: then return true;
! 29: else return false;
! 30: end if;
! 31: end if;
! 32: end At_End;
! 33:
! 34: function Stop ( p : Pred_Pars; t,target : Complex_Number;
! 35: step : double_float ) return boolean is
! 36:
! 37: -- DESCRIPTION :
! 38: -- Returns true if either the step size is smaller than p.minstep, or
! 39: -- or alternatively, in case of geometric predictor, if the distance to
! 40: -- the target has become smaller than p.minstep.
! 41:
! 42: begin
! 43: if step <= p.minstep
! 44: then return true;
! 45: else if (p.predictor_type = 2) or (p.predictor_type = 5)
! 46: then return (AbsVal(t-target) <= p.minstep);
! 47: else return false;
! 48: end if;
! 49: end if;
! 50: end Stop;
! 51:
! 52: function Is_Multiple ( a,b,tol : double_float ) return integer is
! 53:
! 54: -- DESCRIPTION :
! 55: -- Returns a/b if a is a multiple of b, returns 0 in the other case.
! 56:
! 57: fq : double_float;
! 58: iq : integer;
! 59: begin
! 60: if abs(b) < tol
! 61: then return 0;
! 62: else fq := a/b;
! 63: iq := integer(fq);
! 64: if abs(fq - double_float(iq)) <= tol
! 65: then return iq;
! 66: else return 0;
! 67: end if;
! 68: end if;
! 69: end Is_Multiple;
! 70:
! 71: procedure Linear_Step_Control
! 72: ( success : in boolean; p : in Pred_Pars;
! 73: step : in out double_float;
! 74: nsuccess : in out natural; trial : in natural ) is
! 75:
! 76: -- DESCRIPTION :
! 77: -- Control of step size for linear path following.
! 78: -- With geometric prediction, the ratio (=step) will be enlarged
! 79: -- when not success. In order not to break the sequence, the ratio
! 80: -- is not reduced when success.
! 81:
! 82: begin
! 83: if (p.predictor_type = 2) or (p.predictor_type = 5)
! 84: then if success
! 85: then nsuccess := nsuccess + 1;
! 86: else nsuccess := 0;
! 87: if p.expfac < 1.0
! 88: then step := step + p.expfac*(1.0 - step);
! 89: else step := step + (1.0 - step)/p.expfac;
! 90: end if;
! 91: if step >= 1.0
! 92: then step := p.minstep/2.0; -- that ends the game
! 93: end if;
! 94: end if;
! 95: else if success
! 96: then nsuccess := nsuccess + 1;
! 97: if nsuccess > p.success_steps
! 98: then step := step*p.expfac;
! 99: if step > p.maxstep
! 100: then step := p.maxstep;
! 101: end if;
! 102: end if;
! 103: else nsuccess := 0;
! 104: if trial mod 3 = 0
! 105: then step := step*p.redfac;
! 106: end if;
! 107: end if;
! 108: end if;
! 109: end Linear_Step_Control;
! 110:
! 111: procedure Circular_Step_Control
! 112: ( success : in boolean; p : in Pred_Pars; twopi : in double_float;
! 113: step : in out double_float; nsuccess : in out natural ) is
! 114:
! 115: -- DESCRIPTION :
! 116: -- Control of step size for circular path following, note that the
! 117: -- step size should be some multiple of pi.
! 118:
! 119: maxstep : constant double_float := p.maxstep*twopi;
! 120:
! 121: begin
! 122: if success
! 123: then nsuccess := nsuccess + 1;
! 124: if nsuccess > p.success_steps
! 125: then step := step*2.0; -- expansion factor = 2
! 126: if step > maxstep
! 127: then step := maxstep;
! 128: end if;
! 129: end if;
! 130: else nsuccess := 0;
! 131: step := step*0.5; -- reduction factor = 1/2
! 132: end if;
! 133: end Circular_Step_Control;
! 134:
! 135: procedure Set_Corrector_Parameters
! 136: ( c : in out Corr_Pars; eps : double_float ) is
! 137:
! 138: -- DESCRIPTION :
! 139: -- All eps* parameters in c are set to eps.
! 140:
! 141: begin
! 142: c.epsrx := eps; c.epsax := eps; c.epsrf := eps; c.epsaf := eps;
! 143: end Set_Corrector_Parameters;
! 144:
! 145: function End_Game_Corrector_Parameters
! 146: ( current : Corr_Pars; distance,tol : double_float )
! 147: return Corr_Pars is
! 148:
! 149: -- DESCRIPTION :
! 150: -- Returns corrector parameters for the end game of the first or the
! 151: -- second continuation stage, depending on the distance from the target.
! 152:
! 153: res : Corr_Pars := current;
! 154:
! 155: begin
! 156: if distance < tol -- correct further to detect clustering
! 157: then res := current;
! 158: Set_Corrector_Parameters(res,tol);
! 159: else res := Create_End_Game; -- or to move to end game more smoothly
! 160: end if;
! 161: return res;
! 162: end End_Game_Corrector_Parameters;
! 163:
! 164: -- MANAGEMENT OF DATA DURING PATH FOLLOWING :
! 165:
! 166: procedure Linear_Single_Initialize
! 167: ( s : in Solu_Info; p : in Pred_Pars;
! 168: old_t,prev_t : out Complex_Number;
! 169: prev_v,old_solution,prev_solution : out Vector ) is
! 170:
! 171: -- DESCRIPTION :
! 172: -- Initialization for linear path following of one path.
! 173:
! 174: -- ON ENTRY :
! 175: -- s solution at beginning of path;
! 176: -- p predictor parameters.
! 177:
! 178: -- ON RETURN :
! 179: -- old_t back up value for continuation parameter;
! 180: -- prev_t previous value of continuation parameter;
! 181: -- old_solution back up value for solution vector;
! 182: -- prev_solution previous value of solution vector;
! 183:
! 184: begin
! 185: old_t := s.sol.t; old_solution := s.sol.v;
! 186: if p.predictor_type <= 2 -- for all secant predictors
! 187: then prev_t := s.sol.t;
! 188: prev_solution := s.sol.v;
! 189: end if;
! 190: prev_v := (prev_v'range => Create(0.0));
! 191: end Linear_Single_Initialize;
! 192:
! 193: procedure Linear_Single_Management
! 194: ( s : in out Solu_Info; p : in Pred_Pars; c : in Corr_Pars;
! 195: old_t,prev_t : in out Complex_Number;
! 196: old_solution,prev_solution,old_v,prev_v,vv : in out Vector;
! 197: step : in out double_float; nsuccess,trial : in out natural;
! 198: success : in out boolean ) is
! 199:
! 200: -- DESCRIPTION :
! 201: -- Management of data after correction during linear path following.
! 202:
! 203: -- PARAMETERS :
! 204: -- s current solution;
! 205: -- p predictor parameters;
! 206: -- c corrector parameters;
! 207: -- old_t back up value for continuation parameter;
! 208: -- prev_t previous value of continuation parameter;
! 209: -- old_solution back up value for solution vector;
! 210: -- prev_solution previous value of solution vector;
! 211: -- old_v back up value for tangent direction;
! 212: -- prev_v previous value for tangent direction;
! 213: -- vv current tangent direction;
! 214: -- step current step size;
! 215: -- nsuccess number of consecutive successes;
! 216: -- trial number of trials after failure;
! 217: -- success successful correction step.
! 218:
! 219: begin
! 220: success := (((s.resa <= c.epsaf) or else (s.cora <= c.epsax))
! 221: or else ((s.resr <= c.epsrf) or else (s.corr <= c.epsrx)));
! 222: if ((p.predictor_type <= 2) and then success) -- secant predictors
! 223: then prev_t := old_t;
! 224: prev_solution := old_solution;
! 225: end if;
! 226: if ((p.predictor_type = 6) and then success) -- Hermite predictor
! 227: then prev_t := old_t;
! 228: prev_solution := old_solution;
! 229: prev_v := old_v;
! 230: end if;
! 231: if ((p.predictor_type = 1) or (p.predictor_type = 3)) -- complex predictors
! 232: then if success
! 233: then trial := 0;
! 234: else trial := trial + 1;
! 235: end if;
! 236: end if;
! 237: s.nstep := s.nstep + 1;
! 238: if not success
! 239: then s.nfail := s.nfail + 1;
! 240: end if;
! 241: Linear_Step_Control(success,p,step,nsuccess,trial);
! 242: if step < p.minstep
! 243: then return;
! 244: end if;
! 245: if not success
! 246: then s.sol.t := old_t;
! 247: s.sol.v := old_solution;
! 248: else old_t := s.sol.t;
! 249: old_solution := s.sol.v;
! 250: if p.predictor_type = 6
! 251: then old_v := vv;
! 252: end if;
! 253: end if;
! 254: end Linear_Single_Management;
! 255:
! 256: procedure Linear_Multiple_Initialize
! 257: ( s : in Solu_Info_Array; p : in Pred_Pars;
! 258: t,old_t,prev_t : out Complex_Number;
! 259: sa,old_sa,prev_sa : in out Solution_Array ) is
! 260:
! 261: -- DECRIPTION :
! 262: -- Initialization for linear path following for more than one path.
! 263:
! 264: begin
! 265: t := s(s'first).sol.t;
! 266: old_t := s(s'first).sol.t;
! 267: Copy(s,sa); Copy(sa,old_sa);
! 268: case p.predictor_type is
! 269: when 0 | 1 | 2 => prev_t := s(s'first).sol.t; Copy(sa,prev_sa);
! 270: when others => null;
! 271: end case;
! 272: end Linear_Multiple_Initialize;
! 273:
! 274: procedure Linear_Multiple_Management
! 275: ( s : in out Solu_Info_array;
! 276: sa,old_sa,prev_sa : in out Solution_Array;
! 277: t,old_t,prev_t : in out Complex_Number; p : in Pred_Pars;
! 278: step : in out double_float; pivot : in natural;
! 279: nsuccess,trial : in out natural; success : in boolean ) is
! 280:
! 281: -- DESCRIPTION :
! 282: -- Management of data after correction during linear path following.
! 283:
! 284: -- PARAMETERS :
! 285: -- s current solutions with information statistics;
! 286: -- sa current solutions;
! 287: -- old_sa back up value for solutions;
! 288: -- prev_sa previous solutions;
! 289: -- t current value of continuation parameter;
! 290: -- old_t back up value for continuation parameter;
! 291: -- prev_t previous value of continuation parameter;
! 292: -- p predictor parameters;
! 293: -- step current step size;
! 294: -- pivot solution where correction is started;
! 295: -- nsuccess number of consecutive successes;
! 296: -- trial number of trials after failure;
! 297: -- success successful correction step.
! 298:
! 299: begin
! 300: if ((p.predictor_type <= 2) and then success) -- for secant predictors
! 301: then prev_t := old_t; Copy(old_sa,prev_sa);
! 302: end if;
! 303: if ((p.predictor_type = 1) or (p.predictor_type = 3)) -- complex predictors
! 304: then if success
! 305: then trial := 0;
! 306: else trial := trial + 1;
! 307: end if;
! 308: end if;
! 309: for k in s'range loop
! 310: s(k).nstep := s(k).nstep + 1;
! 311: end loop;
! 312: if not success
! 313: then s(pivot).nfail := s(pivot).nfail + 1;
! 314: end if;
! 315: Linear_Step_Control(success,p,step,nsuccess,trial);
! 316: if step < p.minstep then return; end if;
! 317: if success
! 318: then Copy(sa,old_sa); old_t := t;
! 319: else Copy(old_sa,sa); t := old_t;
! 320: end if;
! 321: end Linear_Multiple_Management;
! 322:
! 323: procedure Circular_Management
! 324: ( s : in out Solu_Info; p : in Pred_Pars; c : in Corr_Pars;
! 325: old_t,prev_t : in out Complex_Number;
! 326: start_solution : in Vector;
! 327: old_solution,prev_solution,w_sum,w_all_sum : in out Vector;
! 328: twopi,epslop,tol : in double_float;
! 329: theta,old_theta,step : in out double_float;
! 330: nsuccess,n_sum,n_all_sum,w_c : in out natural;
! 331: max_wc : in natural; stop,success : in out boolean ) is
! 332:
! 333: -- DESCRIPTION :
! 334: -- Management of circular path following.
! 335:
! 336: -- PARAMETERS :
! 337: -- s current solution;
! 338: -- p predictor parameters;
! 339: -- c corrector parameters;
! 340: -- old_t back up value for continuation parameter;
! 341: -- prev_t previous value of continuation parameter;
! 342: -- start_solution solution vector at start of continuation;
! 343: -- old_solution back up value for solution vector;
! 344: -- prev_solution previous value of solution vector;
! 345: -- w_sum sum of cycle;
! 346: -- w_all_sum sum of all cycles;
! 347: -- twopi two times PI;
! 348: -- epslop tolerance to decide whether two vectors are equal;
! 349: -- theta current value of theta;
! 350: -- old_theta back up value for theta;
! 351: -- step current step size;
! 352: -- nsuccess number of consecutive successes;
! 353: -- n_sum number of cycles;
! 354: -- n_all_sum total number of cycles;
! 355: -- w_c winding number;
! 356: -- max_wc upper bound on winding number;
! 357: -- stop true when winding number has been computed;
! 358: -- success successful correction step.
! 359:
! 360: tmp : integer;
! 361:
! 362: begin
! 363: success := (((s.resa <= c.epsaf) or else (s.cora <= c.epsax))
! 364: or else ((s.resr <= c.epsrf) or else (s.corr <= c.epsrx)));
! 365: if p.predictor_type = 0 and then success
! 366: then prev_t := old_t;
! 367: prev_solution := old_solution;
! 368: end if;
! 369: s.nstep := s.nstep + 1;
! 370: if not success then s.nfail := s.nfail + 1; end if;
! 371: if success
! 372: then old_theta := theta; old_t := s.sol.t;
! 373: old_solution := s.sol.v;
! 374: w_all_sum := w_all_sum + s.sol.v;
! 375: n_all_sum := n_all_sum + 1;
! 376: tmp := Is_Multiple(theta,twopi*p.maxstep,tol);
! 377: if tmp /= 0
! 378: then w_sum := w_sum + s.sol.v; n_sum := n_sum + 1;
! 379: end if;
! 380: w_c := Is_Multiple(theta,twopi,tol);
! 381: if w_c /= 0
! 382: then if Equal(s.sol.v,start_solution,epslop)
! 383: then w_sum := w_sum - s.sol.v * Create(0.5);
! 384: w_all_sum := w_all_sum - s.sol.v * Create(0.5);
! 385: stop := true;
! 386: else stop := (w_c >= max_wc);
! 387: end if;
! 388: end if;
! 389: else theta := old_theta;
! 390: s.sol.t := old_t;
! 391: s.sol.v := old_solution;
! 392: end if;
! 393: if not stop
! 394: then Circular_Step_Control(success,p,twopi,step,nsuccess);
! 395: end if;
! 396: end Circular_Management;
! 397:
! 398: -- UPDATE OF PATH DIRECTION :
! 399:
! 400: procedure Update_Direction
! 401: ( proj : in boolean;
! 402: freqcnt,defer,r,m,estm,cntm : in out natural;
! 403: thresm : in natural; er : in out integer;
! 404: t,target : in Complex_Number; x : in Vector;
! 405: dt,s,logs : in out Standard_Floating_Vectors.Vector;
! 406: logx,wvl0,wvl1,wvl2 : in out VecVec;
! 407: v,errv : in out Standard_Floating_Vectors.Vector;
! 408: err : in out double_float ) is
! 409:
! 410: -- DESCRIPTION :
! 411: -- Computes an approximation of the direction of the path.
! 412:
! 413: -- ON ENTRY :
! 414: -- file to write intermediate data on, may be omitted;
! 415: -- proj whether solution vector lies in projective space;
! 416: -- freqcnt counts the frequency of calls;
! 417: -- defer only if freqcnt = defer, calculations are done;
! 418: -- r current order in extrapolation formula;
! 419: -- m current value for multiplicity;
! 420: -- estm current estimated for multiplicity;
! 421: -- cntm number of consecutive times estm has been guessed;
! 422: -- thresm threshold for augmenting m to estm;
! 423: -- er order of extrapolator on the errors;
! 424: -- t current value of continuation parameter;
! 425: -- target target value of continuation parameter;
! 426: -- x current solution vector;
! 427: -- dt vector with distances to target;
! 428: -- s s-values w.r.t. the current value m;
! 429: -- logs logarithms of the s-values;
! 430: -- logx logarithms of the solution vectors;
! 431: -- wvl0 consecutive estimates for previous direction;
! 432: -- wvl1 consecutive estimates for current direction;
! 433: -- wvl2 used as work space for wvl0 and wvl1;
! 434: -- v current approximate direction of the path;
! 435: -- errv vector of errors used to estimate m;
! 436: -- err norm of errv.
! 437:
! 438: -- ON RETURN :
! 439: -- All in-out variables are updated, provided freqcnt = defer.
! 440:
! 441: begin
! 442: if freqcnt >= defer
! 443: then
! 444: freqcnt := 0;
! 445: if proj
! 446: then Projective_Update_Direction
! 447: (r,m,estm,cntm,thresm,er,t,target,x,dt,s,logs,logx,v,errv,err);
! 448: else Affine_Update_Direction
! 449: (r,m,estm,cntm,thresm,er,t,target,x,
! 450: dt,s,logs,logx,wvl0,wvl1,wvl2,v,errv,err);
! 451: end if;
! 452: -- defer := defer + 1; -- that is asking for troubles !
! 453: else
! 454: freqcnt := freqcnt + 1;
! 455: end if;
! 456: end Update_Direction;
! 457:
! 458: procedure Update_Direction
! 459: ( file : in file_type; proj : in boolean;
! 460: freqcnt,defer,r,m,estm,cntm : in out natural;
! 461: thresm : in natural; er : in out integer;
! 462: t,target : in Complex_Number; x : in Vector;
! 463: dt,s,logs : in out Standard_Floating_Vectors.Vector;
! 464: logx,wvl0,wvl1,wvl2 : in out VecVec;
! 465: v,errv : in out Standard_Floating_Vectors.Vector;
! 466: err : in out double_float ) is
! 467:
! 468: -- DESCRIPTION :
! 469: -- Computes an approximation of the direction of the path and produces
! 470: -- intermediate output to the file.
! 471:
! 472: begin
! 473: if freqcnt >= defer
! 474: then
! 475: freqcnt := 0;
! 476: if proj
! 477: then Projective_Update_Direction
! 478: (file,r,m,estm,cntm,thresm,er,t,target,x,
! 479: dt,s,logs,logx,v,errv,err);
! 480: else Affine_Update_Direction
! 481: (file,r,m,estm,cntm,thresm,er,t,target,x,
! 482: dt,s,logs,logx,wvl0,wvl1,wvl2,v,errv,err);
! 483: end if;
! 484: -- defer := defer + 1; -- asking for troubles !
! 485: put(file,"direction : "); put(file,v); new_line(file);
! 486: put(file,"difference to old direction : "); put(file,err);
! 487: new_line(file);
! 488: put(file,"++ current m : "); put(file,m,1); put(file," ++ ");
! 489: put(file,cntm,1); put(file," times estimated m : "); put(file,estm,1);
! 490: put(file," ++ threshold : "); put(file,thresm,1); put_line(file," ++");
! 491: new_line(file);
! 492: else
! 493: freqcnt := freqcnt + 1;
! 494: end if;
! 495: end Update_Direction;
! 496:
! 497: -- LINEAR PATH FOLLOWING FOR ONE PATH :
! 498:
! 499: procedure Linear_Single_Normal_Silent_Continue
! 500: ( s : in out Solu_Info; target : in Complex_Number;
! 501: tol : in double_float; proj : in boolean;
! 502: p : in Pred_Pars; c : in Corr_Pars ) is
! 503:
! 504: old_t,prev_t : Complex_Number;
! 505: old_solution,prev_solution,old_v,prev_v,vv : Vector(s.sol.v'range);
! 506: step : double_float := p.maxstep;
! 507: nsuccess,trial : natural := 0;
! 508: success : boolean := true;
! 509:
! 510: procedure Predictor is new Single_Predictor(Norm,dH,dH);
! 511:
! 512: procedure Corrector ( s : in out Solu_Info; c : in Corr_Pars ) is
! 513:
! 514: procedure Affine_Corrector is
! 515: new Affine_Single_Severe_Normal_Silent_Corrector(Norm,H,dH);
! 516: procedure Projective_Corrector is
! 517: new Projective_Single_Severe_Normal_Silent_Corrector(Norm,H,dH);
! 518:
! 519: begin
! 520: if proj
! 521: then Projective_Corrector(s,c);
! 522: else Affine_Corrector(s,c);
! 523: end if;
! 524: end Corrector;
! 525:
! 526: begin
! 527: Linear_Single_Initialize
! 528: (s,p,old_t,prev_t,prev_v,old_solution,prev_solution);
! 529: while not (At_End(s.sol.t,target,p.dist_target,tol) and success)
! 530: and (s.niter <= c.maxtot) loop
! 531:
! 532: Predictor(s,p,prev_solution,prev_v,vv,prev_t,target,step,tol,trial);
! 533: Corrector(s,c);
! 534: Linear_Single_Management(s,p,c,old_t,prev_t,old_solution,prev_solution,
! 535: old_v,prev_v,vv,step,nsuccess,trial,success);
! 536: if Stop(p,s.sol.t,target,step) then return; end if;
! 537: end loop;
! 538: declare
! 539: cp : Corr_Pars := End_Game_Corrector_Parameters(c,p.dist_target,tol);
! 540: begin
! 541: Corrector(s,cp);
! 542: end;
! 543: end Linear_Single_Normal_Silent_Continue;
! 544:
! 545: procedure Linear_Single_Normal_Reporting_Continue
! 546: ( file : in file_type; s : in out Solu_Info;
! 547: target : in Complex_Number; tol : in double_float;
! 548: proj : in boolean; p : in Pred_Pars; c : in Corr_Pars ) is
! 549:
! 550: old_t,prev_t : Complex_Number;
! 551: old_solution,prev_solution,old_v,prev_v,vv : Vector(s.sol.v'range);
! 552: step : double_float := p.maxstep;
! 553: nsuccess,trial : natural := 0;
! 554: success : boolean := true;
! 555:
! 556: procedure Predictor is new Single_Predictor(Norm,dH,dH);
! 557:
! 558: procedure Corrector ( s : in out Solu_Info; c : in Corr_Pars ) is
! 559:
! 560: procedure Affine_Corrector is
! 561: new Affine_Single_Severe_Normal_Reporting_Corrector(Norm,H,dH);
! 562: procedure Projective_Corrector is
! 563: new Projective_Single_Severe_Normal_Reporting_Corrector(Norm,H,dH);
! 564:
! 565: begin
! 566: if proj
! 567: then Projective_Corrector(file,s,c);
! 568: else Affine_Corrector(file,s,c);
! 569: end if;
! 570: end Corrector;
! 571:
! 572: begin
! 573: Linear_Single_Initialize
! 574: (s,p,old_t,prev_t,prev_v,old_solution,prev_solution);
! 575: sWrite(file,s.sol.all);
! 576: while not (At_End(s.sol.t,target,p.dist_target,tol) and success)
! 577: and (s.niter <= c.maxtot) loop
! 578:
! 579: Predictor(s,p,prev_solution,prev_v,vv,prev_t,target,step,tol,trial);
! 580: if p.predictor_type = 6
! 581: then put_line(file,"previous and current direction : ");
! 582: for i in prev_v'range loop
! 583: put(file,prev_v(i),3,3,3); put(file," ");
! 584: put(file,vv(i),3,3,3); new_line(file);
! 585: end loop;
! 586: end if;
! 587: pWrite(file,step,s.sol.t,s.sol.all);
! 588: Corrector(s,c);
! 589: sWrite(file,s.sol.all);
! 590: Linear_Single_Management(s,p,c,old_t,prev_t,old_solution,prev_solution,
! 591: old_v,prev_v,vv,step,nsuccess,trial,success);
! 592: if Stop(p,s.sol.t,target,step) then return; end if;
! 593: end loop;
! 594: declare
! 595: cp : Corr_Pars := End_Game_Corrector_Parameters(c,p.dist_target,tol);
! 596: begin
! 597: Corrector(s,cp);
! 598: sWrite(file,s.sol.all);
! 599: end;
! 600: end Linear_Single_Normal_Reporting_Continue;
! 601:
! 602: procedure Linear_Single_Conditioned_Silent_Continue
! 603: ( s : in out Solu_Info; target : in Complex_Number;
! 604: tol : in double_float; proj : in boolean;
! 605: rtoric : in natural;
! 606: v : in out Standard_Floating_Vectors.Link_to_Vector;
! 607: errorv : in out double_float;
! 608: p : in Pred_Pars; c : in Corr_Pars ) is
! 609:
! 610: old_t,prev_t : Complex_Number;
! 611: old_solution,prev_solution,old_v,prev_v,vv : Vector(s.sol.v'range);
! 612: step : double_float := p.maxstep;
! 613: nsuccess,trial : natural := 0;
! 614: success : boolean := true;
! 615:
! 616: dls,stp : double_float := 0.0;
! 617: tolsing : constant double_float
! 618: := Continuation_Parameters.tol_endg_inverse_condition;
! 619: diff : Standard_Floating_Vectors.Vector(s.sol.v'range)
! 620: := (s.sol.v'range => 0.0);
! 621:
! 622: r : natural := 0;
! 623: er : integer := 0;
! 624: dt,ds,logs : Standard_Floating_Vectors.Vector(0..rtoric)
! 625: := (0..rtoric => 0.0);
! 626: logx : VecVec(0..rtoric);
! 627: wvl0,wvl1,wvl2 : VecVec(1..rtoric);
! 628: errv : Standard_Floating_Vectors.Vector(0..rtoric) := (0..rtoric => 0.0);
! 629:
! 630: m : natural := 1;
! 631: thresm : natural := p.success_steps;
! 632: estm : natural := m;
! 633: fcnt,cntm : natural := 0;
! 634: defer : natural := thresm;
! 635:
! 636: procedure Predictor is new Single_Predictor(Norm,dH,dH);
! 637:
! 638: procedure Corrector ( s : in out Solu_Info; c : in Corr_Pars ) is
! 639:
! 640: procedure Affine_Corrector is
! 641: new Affine_Single_Severe_Conditioned_Silent_Corrector(Norm,H,dH);
! 642: procedure Projective_Corrector is
! 643: new Projective_Single_Severe_Conditioned_Silent_Corrector(Norm,H,dH);
! 644:
! 645: begin
! 646: if proj
! 647: then Projective_Corrector(s,c);
! 648: else Affine_Corrector(s,c);
! 649: end if;
! 650: end Corrector;
! 651:
! 652: begin
! 653: Linear_Single_Initialize
! 654: (s,p,old_t,prev_t,prev_v,old_solution,prev_solution);
! 655: if rtoric > 0
! 656: then s.sol.m := m;
! 657: end if;
! 658: while not (At_End(s.sol.t,target,p.dist_target,tol) and success)
! 659: and (s.niter <= c.maxtot) loop
! 660:
! 661: if (rtoric > 0)
! 662: then
! 663: if success and then s.rcond > tolsing
! 664: and then (errorv < 100.0) -- avoid divergence
! 665: then Update_Direction
! 666: (proj,fcnt,defer,r,s.sol.m,estm,cntm,thresm,er,s.sol.t,target,
! 667: s.sol.v,dt,ds,logs,logx,wvl0,wvl1,wvl2,v.all,errv,errorv);
! 668: else er := -2;
! 669: end if;
! 670: end if;
! 671:
! 672: Predictor(s,p,prev_solution,prev_v,vv,prev_t,target,step,tol,trial);
! 673: Corrector(s,c);
! 674: Linear_Single_Management(s,p,c,old_t,prev_t,old_solution,prev_solution,
! 675: old_v,prev_v,vv,step,nsuccess,trial,success);
! 676: if Stop(p,s.sol.t,target,step) then return; end if;
! 677: end loop;
! 678: declare
! 679: cp : Corr_Pars := End_Game_Corrector_Parameters(c,p.dist_target,tol);
! 680: begin
! 681: Corrector(s,cp);
! 682: end;
! 683: end Linear_Single_Conditioned_Silent_Continue;
! 684:
! 685: procedure Linear_Single_Conditioned_Reporting_Continue
! 686: ( file : in file_type; s : in out Solu_Info;
! 687: target : in Complex_Number; tol : in double_float;
! 688: proj : in boolean; rtoric : in natural;
! 689: v : in out Standard_Floating_Vectors.Link_to_Vector;
! 690: errorv : in out double_float;
! 691: p : in Pred_Pars; c : in Corr_Pars ) is
! 692:
! 693: old_t,prev_t : Complex_Number;
! 694: step : double_float := p.maxstep;
! 695: nsuccess,trial : natural := 0;
! 696: old_solution,prev_solution,old_v,prev_v,vv : Vector(s.sol.v'range);
! 697: success : boolean := true;
! 698:
! 699: dls,stp : double_float := 0.0;
! 700: tolsing : constant double_float
! 701: := Continuation_Parameters.tol_endg_inverse_condition;
! 702: diff : Standard_Floating_Vectors.Vector(s.sol.v'range)
! 703: := (s.sol.v'range => 0.0);
! 704:
! 705: r : natural := 0;
! 706: er : integer := 0;
! 707: dt,ds,logs : Standard_Floating_Vectors.Vector(0..rtoric)
! 708: := (0..rtoric => 0.0);
! 709: logx : VecVec(0..rtoric);
! 710: wvl0,wvl1,wvl2 : VecVec(1..rtoric);
! 711: errv : Standard_Floating_Vectors.Vector(0..rtoric) := (0..rtoric => 0.0);
! 712:
! 713: m : natural := 1;
! 714: thresm : natural := p.success_steps;
! 715: estm : natural := m;
! 716: fcnt,cntm : natural := 0;
! 717: defer : natural := thresm;
! 718:
! 719: procedure Predictor is new Single_Predictor(Norm,dH,dH);
! 720:
! 721: procedure Corrector ( s : in out Solu_Info; c : in Corr_Pars ) is
! 722:
! 723: procedure Affine_Corrector is
! 724: new Affine_Single_Severe_Conditioned_Reporting_Corrector(Norm,H,dH);
! 725: procedure Projective_Corrector is
! 726: new Projective_Single_Severe_Conditioned_Reporting_Corrector(Norm,H,dH);
! 727:
! 728: begin
! 729: if proj
! 730: then Projective_Corrector(file,s,c);
! 731: else Affine_Corrector(file,s,c);
! 732: end if;
! 733: end Corrector;
! 734:
! 735: begin
! 736: Linear_Single_Initialize
! 737: (s,p,old_t,prev_t,prev_v,old_solution,prev_solution);
! 738: sWrite(file,s.sol.all); -- writing the start solution
! 739: if rtoric > 0
! 740: then s.sol.m := m;
! 741: end if;
! 742: while not (At_End(s.sol.t,target,p.dist_target,tol) and success)
! 743: and (s.niter <= c.maxtot) loop
! 744:
! 745: if (rtoric > 0)
! 746: then
! 747: if success and then s.rcond > tolsing
! 748: and then (errorv < 100.0) -- avoid divergence
! 749: then Update_Direction(file,
! 750: proj,fcnt,defer,r,s.sol.m,estm,cntm,thresm,er,s.sol.t,target,
! 751: s.sol.v,dt,ds,logs,logx,wvl0,wvl1,wvl2,v.all,errv,errorv);
! 752: else er := -2;
! 753: end if;
! 754: end if;
! 755:
! 756: Predictor(s,p,prev_solution,prev_v,vv,prev_t,target,step,tol,trial);
! 757: pWrite(file,step,s.sol.t,s.sol.all);
! 758: Corrector(s,c);
! 759: sWrite(file,s.sol.all);
! 760: Linear_Single_Management(s,p,c,old_t,prev_t,old_solution,prev_solution,
! 761: old_v,prev_v,vv,step,nsuccess,trial,success);
! 762: if Stop(p,s.sol.t,target,step) then return; end if;
! 763: end loop;
! 764: declare
! 765: cp : Corr_Pars := End_Game_Corrector_Parameters(c,p.dist_target,tol);
! 766: begin
! 767: Corrector(s,cp);
! 768: sWrite(file,s.sol.all);
! 769: end;
! 770: end Linear_Single_Conditioned_Reporting_Continue;
! 771:
! 772: -- LINEAR PATH FOLLOWING FOR A NUMBER OF PATHS :
! 773:
! 774: procedure Linear_Multiple_Normal_Silent_Continue
! 775: ( s : in out Solu_Info_Array;
! 776: target : in Complex_Number; tol,dist_sols : in double_float;
! 777: proj : in boolean; p : in Pred_Pars; c : in Corr_Pars ) is
! 778:
! 779: t,old_t,prev_t : Complex_Number;
! 780: sa,old_sa,prev_sa : Solution_Array(s'range);
! 781: step : double_float := p.maxstep;
! 782: cnt,nsuccess,trial : natural := 0;
! 783: pivot : natural := s'first;
! 784: success,fail : boolean := true;
! 785:
! 786: procedure Predictor is new Multiple_Predictor(Norm,dH,dH);
! 787:
! 788: procedure Affine_Corrector is
! 789: new Affine_Multiple_Loose_Normal_Silent_Corrector(Norm,H,dH);
! 790: procedure Projective_Corrector is
! 791: new Projective_Multiple_Loose_Normal_Silent_Corrector(Norm,H,dH);
! 792:
! 793: procedure Correct ( s : in out Solu_Info_Array;
! 794: pivot : in out natural; dist_sols : in double_float;
! 795: c : in Corr_Pars; fail : out boolean ) is
! 796: begin
! 797: Copy(sa,s);
! 798: if proj
! 799: then Projective_Corrector(s,pivot,dist_sols,c,fail);
! 800: else Affine_Corrector(s,pivot,dist_sols,c,fail);
! 801: end if;
! 802: end Correct;
! 803:
! 804: begin
! 805: Linear_Multiple_Initialize(s,p,t,old_t,prev_t,sa,old_sa,prev_sa);
! 806: while not (At_End(t,target,p.dist_target,tol) and success)
! 807: and (s(s'first).niter <= c.maxtot) loop
! 808:
! 809: Predictor(s,p,sa,prev_sa,t,prev_t,target,step,tol,dist_sols,trial);
! 810: Correct(s,pivot,dist_sols,c,fail); Copy(s,sa);
! 811: success := not fail;
! 812: Linear_Multiple_Management(s,sa,old_sa,prev_sa,t,old_t,prev_t,p,step,
! 813: pivot,nsuccess,trial,success);
! 814: if step < p.minstep then return; end if;
! 815: end loop;
! 816: declare
! 817: cp : Corr_Pars := End_Game_Corrector_Parameters(c,p.dist_target,tol);
! 818: begin
! 819: Correct(s,pivot,dist_sols,cp,fail);
! 820: end;
! 821: end Linear_Multiple_Normal_Silent_Continue;
! 822:
! 823: procedure Linear_Multiple_Normal_Reporting_Continue
! 824: ( file : in file_type; s : in out Solu_Info_Array;
! 825: target : in Complex_Number; tol,dist_sols : in double_float;
! 826: proj : in boolean; p : in Pred_Pars; c : in Corr_Pars ) is
! 827:
! 828: t,old_t,prev_t : Complex_Number;
! 829: sa,old_sa,prev_sa : Solution_Array(s'range);
! 830: step : double_float := p.maxstep;
! 831: pivot : natural := s'first;
! 832: cnt,nsuccess,trial : natural := 0;
! 833: success,fail : boolean := true;
! 834:
! 835: procedure Predictor is new Multiple_Predictor(Norm,dH,dH);
! 836:
! 837: procedure Affine_Corrector is
! 838: new Affine_Multiple_Loose_Normal_Reporting_Corrector(Norm,H,dH);
! 839: procedure Projective_Corrector is
! 840: new Projective_Multiple_Loose_Normal_Reporting_Corrector(Norm,H,dH);
! 841:
! 842: procedure Correct ( file : in file_type; s : in out Solu_Info_Array;
! 843: pivot : in out natural; dist_sols : in double_float;
! 844: c : in Corr_Pars; fail : out boolean ) is
! 845: begin
! 846: Copy(sa,s);
! 847: if proj
! 848: then Projective_Corrector(file,s,pivot,dist_sols,c,fail);
! 849: else Affine_Corrector(file,s,pivot,dist_sols,c,fail);
! 850: end if;
! 851: end Correct;
! 852:
! 853: begin
! 854: Linear_Multiple_Initialize(s,p,t,old_t,prev_t,sa,old_sa,prev_sa);
! 855: for k in s'range loop -- write the start solutions
! 856: sWrite(file,sa(k).all);
! 857: end loop;
! 858: while not (At_End(t,target,p.dist_target,tol) and success)
! 859: and (s(s'first).niter <= c.maxtot) loop
! 860:
! 861: Predictor(s,p,sa,prev_sa,t,prev_t,target,step,tol,dist_sols,trial);
! 862: pWrite(file,step,t);
! 863: Correct(file,s,pivot,dist_sols,c,fail); Copy(s,sa);
! 864: success := not fail;
! 865: Linear_Multiple_Management(s,sa,old_sa,prev_sa,t,old_t,prev_t,p,step,
! 866: pivot,nsuccess,trial,success);
! 867: if step < p.minstep then return; end if;
! 868: end loop;
! 869: declare
! 870: cp : Corr_Pars := End_Game_Corrector_Parameters(c,p.dist_target,tol);
! 871: begin
! 872: Correct(file,s,pivot,dist_sols,cp,fail);
! 873: end;
! 874: end Linear_Multiple_Normal_Reporting_Continue;
! 875:
! 876: procedure Linear_Multiple_Conditioned_Silent_Continue
! 877: ( s : in out Solu_Info_Array;
! 878: target : in Complex_Number; tol,dist_sols : in double_float;
! 879: proj : in boolean; p : in Pred_Pars; c : in Corr_Pars ) is
! 880:
! 881: t,old_t,prev_t : Complex_Number;
! 882: sa,old_sa,prev_sa : Solution_Array(s'range);
! 883: step : double_float := p.maxstep;
! 884: pivot : natural := s'first;
! 885: cnt,nsuccess,trial : natural := 0;
! 886: success,fail : boolean := true;
! 887:
! 888: procedure Predictor is new Multiple_Predictor(Norm,dH,dH);
! 889:
! 890: procedure Affine_Corrector is
! 891: new Affine_Multiple_Loose_Conditioned_Silent_Corrector(Norm,H,dH);
! 892: procedure Projective_Corrector is
! 893: new Projective_Multiple_Loose_Conditioned_Silent_Corrector(Norm,H,dH);
! 894:
! 895: procedure Correct ( s : in out Solu_Info_Array;
! 896: pivot : in out natural; dist_sols : in double_float;
! 897: c : in Corr_Pars; fail : out boolean ) is
! 898: begin
! 899: Copy(sa,s);
! 900: if proj
! 901: then Projective_Corrector(s,pivot,dist_sols,c,fail);
! 902: else Affine_Corrector(s,pivot,dist_sols,c,fail);
! 903: end if;
! 904: end Correct;
! 905:
! 906: begin
! 907: Linear_Multiple_Initialize(s,p,t,old_t,prev_t,sa,old_sa,prev_sa);
! 908: while not (At_End(t,target,p.dist_target,tol) and success)
! 909: and (s(s'first).niter <= c.maxtot) loop
! 910:
! 911: Predictor(s,p,sa,prev_sa,t,prev_t,target,step,tol,dist_sols,trial);
! 912: Correct(s,pivot,dist_sols,c,fail); Copy(s,sa);
! 913: success := not fail;
! 914: Linear_Multiple_Management(s,sa,old_sa,prev_sa,t,old_t,prev_t,p,step,
! 915: pivot,nsuccess,trial,success);
! 916: if step < p.minstep then return; end if;
! 917: end loop;
! 918: declare
! 919: cp : Corr_Pars := End_Game_Corrector_Parameters(c,p.dist_target,tol);
! 920: begin
! 921: Correct(s,pivot,dist_sols,cp,fail);
! 922: end;
! 923: end Linear_Multiple_Conditioned_Silent_Continue;
! 924:
! 925: procedure Linear_Multiple_Conditioned_Reporting_Continue
! 926: ( file : in file_type; s : in out Solu_Info_Array;
! 927: target : in Complex_Number; tol,dist_sols : in double_float;
! 928: proj : in boolean; p : in Pred_Pars; c : in Corr_Pars ) is
! 929:
! 930: t,old_t,prev_t : Complex_Number;
! 931: sa,old_sa,prev_sa : Solution_Array(s'range);
! 932: step : double_float := p.maxstep;
! 933: pivot : natural := s'first;
! 934: cnt,nsuccess,trial : natural := 0;
! 935: success,fail : boolean := true;
! 936:
! 937: procedure Predictor is new Multiple_Predictor(Norm,dH,dH);
! 938:
! 939: procedure Affine_Corrector is
! 940: new Affine_Multiple_Loose_Conditioned_Reporting_Corrector(Norm,H,dH);
! 941: procedure Projective_Corrector is
! 942: new Projective_Multiple_Loose_Conditioned_Reporting_Corrector(Norm,H,dH);
! 943:
! 944: procedure Correct ( file : in file_type; s : in out Solu_Info_Array;
! 945: pivot : in out natural; dist_sols : in double_float;
! 946: c : in Corr_Pars; fail : out boolean ) is
! 947: begin
! 948: Copy(sa,s);
! 949: if proj
! 950: then Projective_Corrector(file,s,pivot,dist_sols,c,fail);
! 951: else Affine_Corrector(file,s,pivot,dist_sols,c,fail);
! 952: end if;
! 953: end Correct;
! 954:
! 955: begin
! 956: Linear_Multiple_Initialize(s,p,t,old_t,prev_t,sa,old_sa,prev_sa);
! 957: for k in s'range loop -- write start solutions
! 958: sWrite(file,sa(k).all);
! 959: end loop;
! 960: while not (At_End(t,target,p.dist_target,tol) and success)
! 961: and (s(s'first).niter <= c.maxtot) loop
! 962:
! 963: Predictor(s,p,sa,prev_sa,t,prev_t,target,step,tol,dist_sols,trial);
! 964: pWrite(file,step,t);
! 965: Correct(file,s,pivot,dist_sols,c,fail); Copy(s,sa);
! 966: success := not fail;
! 967: Linear_Multiple_Management(s,sa,old_sa,prev_sa,t,old_t,prev_t,p,step,
! 968: pivot,nsuccess,trial,success);
! 969: if step < p.minstep then return; end if;
! 970: end loop;
! 971: declare
! 972: cp : Corr_Pars := End_Game_Corrector_Parameters(c,p.dist_target,tol);
! 973: begin
! 974: Correct(file,s,pivot,dist_sols,cp,fail);
! 975: end;
! 976: end Linear_Multiple_Conditioned_Reporting_Continue;
! 977:
! 978: -- CIRCULAR PATH FOLLOWING FOR ONE PATH :
! 979:
! 980: procedure Circular_Single_Normal_Reporting_Continue
! 981: ( file : in file_type; s : in out Solu_Info;
! 982: target : in Complex_Number; tol,epslop : in double_float;
! 983: wc : out natural; max_wc : in natural; sum,all_sum : out Vector;
! 984: proj : in boolean; p : in Pred_Pars; c : in Corr_Pars ) is
! 985:
! 986: old_t,prev_t : Complex_Number;
! 987: t0_min_target : Complex_Number := s.sol.t - target;
! 988: theta,old_theta : double_float := 0.0;
! 989: twopi : constant double_float := 2.0*PI;
! 990: step : double_float := twopi*p.maxstep;
! 991: old_solution,prev_solution,start_solution : Vector(s.sol.v'range);
! 992: w_c,nsuccess : natural := 0;
! 993: success : boolean := true;
! 994: stop : boolean := false;
! 995: w_sum,w_all_sum : Vector(s.sol.v'range) := Create(0.5)*s.sol.v;
! 996: n_sum,n_all_sum : natural := 0;
! 997:
! 998: procedure T_C_Predictor is new Tangent_Circular_Predictor(Norm,dH,dH);
! 999:
! 1000: procedure Affine_Corrector is
! 1001: new Affine_Single_Severe_Normal_Reporting_Corrector(Norm,H,dH);
! 1002: procedure Projective_Corrector is
! 1003: new Projective_Single_Severe_Normal_Reporting_Corrector(Norm,H,dH);
! 1004:
! 1005: begin
! 1006: old_t := s.sol.t; old_solution := s.sol.v; -- INITIALIZATION
! 1007: start_solution := s.sol.v;
! 1008: if p.predictor_type = 0
! 1009: then prev_t := s.sol.t; prev_solution := s.sol.v;
! 1010: end if;
! 1011: sWrite(file,s.sol.all); -- write the start solution
! 1012: while (s.niter <= c.maxtot) loop
! 1013: case p.predictor_type is
! 1014: when 0 => Secant_Circular_Predictor(s.sol.v,prev_solution,s.sol.t,theta,
! 1015: prev_t,t0_min_target,target,step,tol);
! 1016: when 2 => T_C_Predictor(s.sol.v,s.sol.t,theta, t0_min_target,target,
! 1017: step,tol);
! 1018: s.nsyst := s.nsyst + 1;
! 1019: when others => null; -- these choices make no sense !!!
! 1020: end case;
! 1021: pWrite(file,step,s.sol.t,s.sol.all);
! 1022: if proj
! 1023: then Projective_Corrector(file,s,c);
! 1024: else Affine_Corrector(file,s,c);
! 1025: end if;
! 1026: sWrite(file,s.sol.all);
! 1027: Circular_Management
! 1028: (s,p,c,old_t,prev_t,start_solution,old_solution,prev_solution,
! 1029: w_sum,w_all_sum,twopi,epslop,tol,theta,old_theta,
! 1030: step,nsuccess,n_sum,n_all_sum,w_c,max_wc,stop,success);
! 1031: exit when stop;
! 1032: if step < p.minstep then return; end if;
! 1033: end loop;
! 1034: wc := w_c;
! 1035: if n_all_sum /= 0
! 1036: then all_sum := w_all_sum*Create(1.0/double_float(n_all_sum));
! 1037: end if;
! 1038: if n_sum /= 0
! 1039: then sum := w_sum*Create(1.0/double_float(n_sum));
! 1040: elsif n_all_sum /= 0
! 1041: then all_sum := w_all_sum*Create(1.0/double_float(n_all_sum));
! 1042: end if;
! 1043: end Circular_Single_Normal_Reporting_Continue;
! 1044:
! 1045: procedure Circular_Single_Conditioned_Reporting_Continue
! 1046: ( file : in file_type; s : in out Solu_Info;
! 1047: target : in Complex_Number; tol,epslop : in double_float;
! 1048: wc : out natural; max_wc : in natural; sum,all_sum : out Vector;
! 1049: proj : in boolean; p : in Pred_Pars; c : in Corr_Pars ) is
! 1050:
! 1051: old_t,prev_t : Complex_Number;
! 1052: theta,old_theta : double_float := 0.0;
! 1053: twopi : constant double_float := 2.0*PI;
! 1054: step : double_float := twopi*p.maxstep;
! 1055: t0_min_target : Complex_Number := s.sol.t - target;
! 1056: old_solution,prev_solution,start_solution : Vector(s.sol.v'range);
! 1057: w_c,nsuccess : natural := 0;
! 1058: success : boolean := true;
! 1059: stop : boolean := false;
! 1060: w_sum,w_all_sum : Vector(s.sol.v'range) := Create(0.5)*s.sol.v;
! 1061: n_sum,n_all_sum : natural := 0;
! 1062:
! 1063: procedure T_C_Predictor is new Tangent_Circular_Predictor(Norm,dH,dH);
! 1064:
! 1065: procedure Affine_Corrector is
! 1066: new Affine_Single_Severe_Conditioned_Reporting_Corrector(Norm,H,dH);
! 1067: procedure Projective_Corrector is
! 1068: new Projective_Single_Severe_Conditioned_Reporting_Corrector(Norm,H,dH);
! 1069:
! 1070: begin
! 1071: old_t := s.sol.t; old_solution := s.sol.v; -- INITIALIZATION
! 1072: start_solution := s.sol.v;
! 1073: if p.predictor_type = 0
! 1074: then prev_t := s.sol.t; prev_solution := old_solution;
! 1075: end if;
! 1076: sWrite(file,s.sol.all); -- writing the start solution
! 1077: while s.niter <= c.maxtot loop
! 1078: case p.predictor_type is
! 1079: when 0 => Secant_Circular_Predictor(s.sol.v,prev_solution,s.sol.t,theta,
! 1080: prev_t,t0_min_target,target,step,tol);
! 1081: when 2 => T_C_Predictor(s.sol.v,s.sol.t,theta,t0_min_target,
! 1082: target,step,tol);
! 1083: s.nsyst := s.nsyst + 1;
! 1084: when others => null; -- these choices make no sense !!!
! 1085: end case;
! 1086: pWrite(file,step,s.sol.t,s.sol.all);
! 1087: if proj
! 1088: then Projective_Corrector(file,s,c);
! 1089: else Affine_Corrector(file,s,c);
! 1090: end if;
! 1091: sWrite(file,s.sol.all);
! 1092: Circular_Management
! 1093: (s,p,c,old_t,prev_t,start_solution,old_solution,prev_solution,
! 1094: w_sum,w_all_sum,twopi,epslop,tol,theta,old_theta,
! 1095: step,nsuccess,n_sum,n_all_sum,w_c,max_wc,stop,success);
! 1096: exit when stop;
! 1097: if step < p.minstep then return; end if;
! 1098: end loop;
! 1099: wc := w_c;
! 1100: if n_all_sum /= 0
! 1101: then all_sum := w_all_sum*Create(1.0/double_float(n_all_sum));
! 1102: end if;
! 1103: if n_sum /= 0
! 1104: then sum := w_sum*Create(1.0/double_float(n_sum));
! 1105: elsif n_all_sum /= 0
! 1106: then sum := w_all_sum*Create(1.0/double_float(n_all_sum));
! 1107: end if;
! 1108: end Circular_Single_Conditioned_Reporting_Continue;
! 1109:
! 1110: end Path_Trackers;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to path_trackers.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