Annotation of OpenXM_contrib/PHC/Ada/Continuation/path_trackers.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_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