Annotation of OpenXM_contrib/PHC/Ada/Continuation/predictors.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
2: with Standard_Natural_Vectors;
3: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
4:
5: package body Predictors is
6:
7: -- PREDICTORS for t :
8:
9: procedure Real_Predictor
10: ( t : in out Complex_Number; target : in Complex_Number;
11: h,tol : in double_float; pow : in positive := 1;
12: hh : out double_float ) is
13:
14: nt : Complex_Number;
15:
16: begin
17: if pow = 1
18: then nt := t + Create(h);
19: else nt := Create((h+REAL_PART(t)**pow)**(1.0/double_float(pow)));
20: end if;
21: if REAL_PART(nt) >= REAL_PART(target)
22: or else abs( REAL_PART(nt) - REAL_PART(target) ) <= tol
23: then hh := REAL_PART(target) - REAL_PART(t);
24: t := target;
25: else hh := h;
26: t := nt;
27: end if;
28: end Real_Predictor;
29:
30: procedure Complex_Predictor
31: ( t : in out Complex_Number; target : in Complex_Number;
32: h,tol : in double_float; hh : out double_float;
33: distance : in double_float; trial : in natural ) is
34:
35: nt,target_direction,step : Complex_Number;
36: dre,dim,d,x,y,alfa,absnt : double_float;
37: tr3 : natural range 0..2;
38:
39: begin
40: target_direction := target - t;
41: d := AbsVal(target_direction);
42: if d < tol
43: then nt := target - Create(distance);
44: hh := AbsVal(nt - t);
45: else tr3 := trial mod 3;
46: case tr3 is
47: when 0 => target_direction := target_direction / Create(d);
48: if (d - h) > distance
49: then step := Create(h) * target_direction;
50: else step := (d - distance) * target_direction;
51: end if;
52: nt := t + step;
53: when 1 => dre := REAL_PART(target_direction);
54: dim := IMAG_PART(target_direction);
55: alfa := ARCTAN(dim/dre);
56: x := h * COS(alfa + PI/4.0);
57: y := h * SIN(alfa + PI/4.0);
58: step := Create(x,y);
59: nt := t + step;
60: absnt := AbsVal(nt);
61: if (absnt > AbsVal(target))
62: or else (AbsVal(nt - target) < distance)
63: then step := Create(0.0,y) * Create(h);
64: nt := t + step;
65: end if;
66: when 2 => dre := REAL_PART(target_direction);
67: dim := IMAG_PART(target_direction);
68: alfa := ARCTAN(dim/dre);
69: x := h * COS(alfa - PI/4.0);
70: y := h * SIN(alfa - PI/4.0);
71: step := Create(x,y);
72: nt := t + step;
73: absnt := AbsVal(nt);
74: if (absnt > AbsVal(target))
75: or else (AbsVal(nt - target) < distance)
76: then step := Create(0.0,-y) * Create(h);
77: nt := t + step;
78: end if;
79: end case;
80: hh := AbsVal(step);
81: end if;
82: t := nt;
83: end Complex_Predictor;
84:
85: procedure Circular_Predictor
86: ( t : in out Complex_Number; theta : in out double_float;
87: t0_min_target,target : in Complex_Number;
88: h : in double_float ) is
89:
90: e_i_theta : Complex_Number;
91:
92: begin
93: theta := theta + h;
94: e_i_theta := Create(COS(theta),SIN(theta));
95: t := target + t0_min_target * e_i_theta;
96: end Circular_Predictor;
97:
98: procedure Geometric_Predictor
99: ( t : in out Complex_Number; target : in Complex_Number;
100: h,tol : in double_float ) is
101:
102: nt : Complex_Number;
103:
104: begin
105: nt := target - Create(h)*(target - t);
106: if abs( REAL_PART(nt) - REAL_PART(target) ) <= tol
107: then t := target;
108: else t := nt;
109: end if;
110: end Geometric_Predictor;
111:
112: -- PREDICTORS FOR x :
113:
114: procedure Secant_Predictor
115: ( x : in out Solution_Array; prev_x : in Solution_Array;
116: fac : in Complex_Number; dist_x : in double_float ) is
117:
118: j : natural;
119: xx : Solution_Array(x'range);
120:
121: begin
122: Copy(x,xx);
123: for i in x'range loop
124: x(i).v := x(i).v + fac * ( x(i).v - prev_x(i).v );
125: j := xx'first;
126: Equals(xx,x(i).v,i,dist_x,j);
127: if j /= i
128: then Copy(xx,x);
129: exit;
130: end if;
131: end loop;
132: Clear(xx);
133: end Secant_Predictor;
134:
135: -- PREDICTORS FOR t AND x :
136:
137: procedure Secant_Single_Real_Predictor
138: ( x : in out Vector; prev_x : in Vector;
139: t : in out Complex_Number; prev_t,target : in Complex_Number;
140: h,tol : in double_float; pow : in positive := 1 ) is
141:
142: hh,tmp : double_float;
143: factor : Complex_Number;
144:
145: begin
146: tmp := AbsVal(t - prev_t);
147: Real_Predictor(t,target,h,tol,pow,hh);
148: if tmp > tol
149: then factor := Create(hh/tmp);
150: x := x + factor * ( x - prev_x );
151: end if;
152: end Secant_Single_Real_Predictor;
153:
154: procedure Secant_Multiple_Real_Predictor
155: ( x : in out Solution_Array; prev_x : in Solution_Array;
156: t : in out Complex_Number; prev_t,target : in Complex_Number;
157: h,tol,dist_x : in double_float; pow : in positive := 1 ) is
158:
159: hh,tmp : double_float;
160: factor : Complex_Number;
161:
162: begin
163: tmp := AbsVal(t - prev_t);
164: Real_Predictor(t,target,h,tol,pow,hh);
165: if tmp > tol
166: then factor := Create(hh/tmp);
167: Secant_Predictor(x,prev_x,factor,dist_x);
168: end if;
169: for k in x'range loop
170: x(k).t := t;
171: end loop;
172: end Secant_Multiple_Real_Predictor;
173:
174: procedure Secant_Single_Complex_Predictor
175: ( x : in out Vector; prev_x : in Vector;
176: t : in out Complex_Number; prev_t,target : in Complex_Number;
177: h,tol,dist_t : in double_float; trial : in natural ) is
178:
179: hh,tmp : double_float;
180: factor : Complex_Number;
181:
182: begin
183: tmp := AbsVal(t - prev_t);
184: Complex_Predictor(t,target,h,tol,hh,dist_t,trial);
185: if tmp > tol
186: then factor := Create(hh/tmp);
187: x := x + factor * ( x - prev_x );
188: end if;
189: end Secant_Single_Complex_Predictor;
190:
191: procedure Secant_Multiple_Complex_Predictor
192: ( x : in out Solution_Array; prev_x : in Solution_Array;
193: t : in out Complex_Number; prev_t,target : in Complex_Number;
194: h,tol,dist_x,dist_t : in double_float; trial : in natural ) is
195:
196: hh,tmp : double_float;
197: factor : Complex_Number;
198:
199: begin
200: tmp := AbsVal(t - prev_t);
201: Complex_Predictor(t,target,h,tol,hh,dist_t,trial);
202: if tmp > tol
203: then factor := Create(hh/tmp);
204: Secant_Predictor(x,prev_x,factor,dist_x);
205: end if;
206: for k in x'range loop
207: x(k).t := t;
208: end loop;
209: end Secant_Multiple_Complex_Predictor;
210:
211: procedure Secant_Circular_Predictor
212: ( x : in out Vector; prev_x : in Vector;
213: t : in out Complex_Number; theta : in out double_float;
214: prev_t,t0_min_target,target : in Complex_Number;
215: h,tol : in double_float ) is
216:
217: factor : Complex_Number;
218: tmp : double_float := AbsVal(t-prev_t);
219:
220: begin
221: if tmp <= tol
222: then Circular_Predictor(t,theta,t0_min_target,target,h);
223: else factor := Create(h/tmp);
224: Circular_Predictor(t,theta,t0_min_target,target,h);
225: x := x + factor * ( x - prev_x );
226: end if;
227: end Secant_Circular_Predictor;
228:
229: procedure Secant_Geometric_Predictor
230: ( x : in out Vector; prev_x : in Vector;
231: t : in out Complex_Number; prev_t,target : in Complex_Number;
232: h,tol : in double_float ) is
233:
234: dist_prev,dist : double_float;
235: factor,tmp : Complex_Number;
236:
237: begin
238: dist_prev := AbsVal(t - prev_t); -- old stepsize
239: tmp := t;
240: Geometric_Predictor(t,target,h,tol);
241: dist := AbsVal(t - tmp); -- new stepsize
242: if dist_prev > tol
243: then factor := Create(dist/dist_prev);
244: x := x + factor * ( x - prev_x );
245: end if;
246: end Secant_Geometric_Predictor;
247:
248: procedure Tangent_Single_Real_Predictor
249: ( x : in out Vector; t : in out Complex_Number;
250: target : in Complex_Number; h,tol : in double_float;
251: pow : in positive := 1 ) is
252:
253: n : natural := x'last;
254: info : integer;
255: norm_tan,hh : double_float;
256: rhs : Vector(x'range);
257: j : Matrix(1..n,1..n);
258: ipvt : Standard_Natural_Vectors.Vector(1..n);
259: prev_t : Complex_Number := t;
260:
261: begin
262: Real_Predictor(t,target,h,tol,pow,hh);
263: j := dH(x,prev_t);
264: lufac(j,n,ipvt,info);
265: if info = 0
266: then rhs := dH(x,prev_t);
267: Min(rhs);
268: lusolve(j,n,ipvt,rhs);
269: norm_tan := Norm(rhs);
270: if norm_tan > tol
271: then hh := hh / norm_tan;
272: Mul(rhs,Create(hh));
273: Add(x,rhs);
274: end if;
275: end if;
276: end Tangent_Single_Real_Predictor;
277:
278: procedure Tangent_Multiple_Real_Predictor
279: ( x : in out Solution_Array; t : in out Complex_Number;
280: target : in Complex_Number; h,tol,dist_x : in double_float;
281: nsys : in out natural; pow : in positive := 1 ) is
282:
283: n : natural := x(x'first).n;
284: norm_tan,hh : double_float;
285: rhs : Vector(1..n);
286: info : integer;
287: j : Matrix(1..n,1..n);
288: ipvt : Standard_Natural_Vectors.Vector(1..n);
289: xx : Solution_Array(x'range);
290: jj : natural;
291: prev_t : Complex_Number := t;
292:
293: begin
294: Real_Predictor(t,target,h,tol,pow,hh);
295: Copy(x,xx);
296: for i in x'range loop
297: j := dH(x(i).v,prev_t);
298: lufac(j,n,ipvt,info);
299: if info = 0
300: then rhs := dH(x(i).v,prev_t);
301: Min(rhs);
302: lusolve(j,n,ipvt,rhs);
303: nsys := nsys + 1;
304: norm_tan := Norm(rhs);
305: if norm_tan > tol
306: then hh := hh / norm_tan;
307: Mul(rhs,Create(hh));
308: Add(x(i).v,rhs);
309: end if;
310: jj := xx'first;
311: Equals(xx,x(i).v,i,dist_x,jj);
312: if jj /= i
313: then Copy(xx,x); exit;
314: end if;
315: else Copy(xx,x); exit;
316: end if;
317: end loop;
318: Clear(xx);
319: for k in x'range loop
320: x(k).t := t;
321: end loop;
322: end Tangent_Multiple_Real_Predictor;
323:
324: procedure Tangent_Single_Complex_Predictor
325: ( x : in out Vector; t : in out Complex_Number;
326: target : in Complex_Number;
327: h,tol,dist_t : in double_float; trial : in natural) is
328:
329: n : natural := x'last;
330: info : integer;
331: norm_tan,hh : double_float;
332: rhs : Vector(x'range);
333: j : Matrix(1..n,1..n);
334: ipvt : Standard_Natural_Vectors.Vector(1..n);
335: prev_t : Complex_Number := t;
336:
337: begin
338: Complex_Predictor(t,target,h,tol,hh,dist_t,trial);
339: j := dH(x,prev_t);
340: lufac(j,n,ipvt,info);
341: if info = 0
342: then rhs := dH(x,prev_t);
343: Min(rhs);
344: lusolve(j,n,ipvt,rhs);
345: norm_tan := Norm(rhs);
346: if norm_tan > tol
347: then hh := hh / norm_tan;
348: Mul(rhs,Create(hh));
349: Add(x,rhs);
350: end if;
351: end if;
352: end Tangent_Single_Complex_Predictor;
353:
354: procedure Tangent_Multiple_Complex_Predictor
355: ( x : in out Solution_Array; t : in out Complex_Number;
356: target : in Complex_Number;
357: h,tol,dist_x,dist_t : in double_float;
358: trial : in natural; nsys : in out natural) is
359:
360: n : natural := x(x'first).n;
361: norm_tan,hh : double_float;
362: rhs : Vector(1..n);
363: info : integer;
364: j : Matrix(1..n,1..n);
365: ipvt : Standard_Natural_Vectors.Vector(1..n);
366: xx : Solution_Array(x'range);
367: jj : natural;
368: prev_t : Complex_Number := t;
369:
370: begin
371: Complex_Predictor(t,target,h,tol,hh,dist_t,trial);
372: Copy(x,xx);
373: for i in x'range loop
374: j := dH(x(i).v,prev_t);
375: lufac(j,n,ipvt,info);
376: if info = 0
377: then rhs := dH(x(i).v,prev_t);
378: Min(rhs);
379: lusolve(j,n,ipvt,rhs);
380: nsys := nsys + 1;
381: norm_tan := Norm(rhs);
382: if norm_tan > tol
383: then hh := hh / norm_tan;
384: Mul(rhs,Create(hh));
385: Add(x(i).v,rhs);
386: end if;
387: jj := xx'first;
388: Equals(xx,x(i).v,i,dist_x,jj);
389: if jj /= i
390: then Copy(xx,x); exit;
391: end if;
392: else Copy(xx,x); exit;
393: end if;
394: end loop;
395: Clear(xx);
396: for k in x'range loop
397: x(k).t := t;
398: end loop;
399: end Tangent_Multiple_Complex_Predictor;
400:
401: procedure Tangent_Circular_Predictor
402: ( x : in out Vector; t : in out Complex_Number;
403: theta : in out double_float;
404: t0_min_target,target : in Complex_Number;
405: h,tol : in double_float ) is
406:
407: n : natural := x'last;
408: info : integer;
409: norm_tan : double_float;
410: rhs : Vector(x'range);
411: j : Matrix(1..n,1..n);
412: ipvt : Standard_Natural_Vectors.Vector(1..n);
413:
414: begin
415: lufac(j,n,ipvt,info);
416: if info = 0
417: then rhs := dH(x,t);
418: Min(rhs);
419: lusolve(j,n,ipvt,rhs);
420: end if;
421: Circular_Predictor(t,theta,t0_min_target,target,h);
422: if info = 0
423: then norm_tan := Norm(rhs);
424: if norm_tan > tol
425: then Mul(rhs,Create(h/norm_tan));
426: Add(x,rhs);
427: end if;
428: end if;
429: end Tangent_Circular_Predictor;
430:
431: procedure Tangent_Geometric_Predictor
432: ( x : in out Vector; t : in out Complex_Number;
433: target : in Complex_Number; h,tol : in double_float ) is
434:
435: n : natural := x'last;
436: info : integer;
437: norm_tan,step : double_float;
438: rhs : Vector(x'range);
439: j : Matrix(1..n,1..n);
440: ipvt : Standard_Natural_Vectors.Vector(1..n);
441: prev_t : Complex_Number := t;
442:
443: begin
444: Geometric_Predictor(t,target,h,tol);
445: j := dH(x,prev_t);
446: lufac(j,n,ipvt,info);
447: if info = 0
448: then rhs := dH(x,prev_t);
449: Min(rhs);
450: lusolve(j,n,ipvt,rhs);
451: norm_tan := Norm(rhs);
452: if norm_tan > tol
453: then step := AbsVal(t-prev_t) / norm_tan;
454: Mul(rhs,Create(step));
455: Add(x,rhs);
456: end if;
457: end if;
458: end Tangent_Geometric_Predictor;
459:
460: function Hermite ( t0,t1,t,x0,x1,v0,v1 : Complex_Number )
461: return Complex_Number is
462:
463: -- DESCRIPTION :
464: -- Returns the value of the third degree interpolating polynomial x(t),
465: -- such that x(t0) = x0, x(t1) = x1, x'(t0) = v0 and x'(t1) = v1.
466:
467: -- REQUIRED : t0 /= t1.
468:
469: -- IMPLEMENTATION :
470: -- x(t) = a3*t^3 + a2*t^2 + a1*t + a0,
471: -- The four interpolation conditions lead to a linear system in
472: -- the coefficients of x(t). This system is first solved explicitly
473: -- and then the polynomial x(t) is evaluated.
474:
475: a0,a1,a2,a3,t10,v10 : Complex_Number;
476:
477: begin
478: t10 := t1 - t0;
479: v10 := (x1 - x0)/t10;
480: a3 := (v1 + v0 - Create(2.0)*v10)/(t10**2);
481: a2 := (v10 - v0 - (t1**2 + t1*t0 - Create(2.0)*t0**2)*a3)/t10;
482: a1 := v0 - (Create(3.0)*a3*t0 + Create(2.0)*a2)*t0;
483: a0 := x0 - ((a3*t0 + a2)*t0 + a1)*t0;
484: return (((a3*t + a2)*t + a1)*t + a0);
485: end Hermite;
486:
487: function Hermite ( t0,t1,t : Complex_Number; x0,x1,v0,v1 : Vector )
488: return Vector is
489:
490: -- DESCRIPTION :
491: -- Returns the value of the third degree interpolating polynomial x(t),
492: -- such that x(t0) = x0, x(t1) = x1, x'(t0) = v0 and x'(t1) = v1,
493: -- for every component.
494:
495: -- REQUIRED : t0 /= t1.
496:
497: res : Vector(x0'range);
498:
499: begin
500: for i in res'range loop
501: res(i) := Hermite(t0,t1,t,x0(i),x1(i),v0(i),v1(i));
502: end loop;
503: return res;
504: end Hermite;
505:
506: procedure Hermite_Single_Real_Predictor
507: ( x : in out Vector; prev_x : in Vector;
508: t : in out Complex_Number; prev_t,target : in Complex_Number;
509: v : in out Vector; prev_v : in Vector;
510: h,tol : in double_float; pow : in positive := 1 ) is
511:
512: n : natural := x'last;
513: info : integer;
514: hh : double_float;
515: j : Matrix(1..n,1..n);
516: ipvt : Standard_Natural_Vectors.Vector(1..n);
517: t1 : Complex_Number := t;
518:
519: begin
520: Real_Predictor(t,target,h,tol,pow,hh);
521: j := dH(x,t1);
522: lufac(j,n,ipvt,info);
523: if info = 0
524: then v := dH(x,t1);
525: Min(v);
526: lusolve(j,n,ipvt,v);
527: if AbsVal(prev_t - t1) > tol
528: then x := Hermite(prev_t,t1,t,prev_x,x,prev_v,v);
529: end if;
530: end if;
531: end Hermite_Single_Real_Predictor;
532:
533: end Predictors;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to predictors.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