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