File: [local] / OpenXM_contrib / PHC / Ada / Continuation / predictors.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:22 2000 UTC (23 years, 8 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
with Standard_Natural_Vectors;
with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
package body Predictors is
-- PREDICTORS for t :
procedure Real_Predictor
( t : in out Complex_Number; target : in Complex_Number;
h,tol : in double_float; pow : in positive := 1;
hh : out double_float ) is
nt : Complex_Number;
begin
if pow = 1
then nt := t + Create(h);
else nt := Create((h+REAL_PART(t)**pow)**(1.0/double_float(pow)));
end if;
if REAL_PART(nt) >= REAL_PART(target)
or else abs( REAL_PART(nt) - REAL_PART(target) ) <= tol
then hh := REAL_PART(target) - REAL_PART(t);
t := target;
else hh := h;
t := nt;
end if;
end Real_Predictor;
procedure Complex_Predictor
( t : in out Complex_Number; target : in Complex_Number;
h,tol : in double_float; hh : out double_float;
distance : in double_float; trial : in natural ) is
nt,target_direction,step : Complex_Number;
dre,dim,d,x,y,alfa,absnt : double_float;
tr3 : natural range 0..2;
begin
target_direction := target - t;
d := AbsVal(target_direction);
if d < tol
then nt := target - Create(distance);
hh := AbsVal(nt - t);
else tr3 := trial mod 3;
case tr3 is
when 0 => target_direction := target_direction / Create(d);
if (d - h) > distance
then step := Create(h) * target_direction;
else step := (d - distance) * target_direction;
end if;
nt := t + step;
when 1 => dre := REAL_PART(target_direction);
dim := IMAG_PART(target_direction);
alfa := ARCTAN(dim/dre);
x := h * COS(alfa + PI/4.0);
y := h * SIN(alfa + PI/4.0);
step := Create(x,y);
nt := t + step;
absnt := AbsVal(nt);
if (absnt > AbsVal(target))
or else (AbsVal(nt - target) < distance)
then step := Create(0.0,y) * Create(h);
nt := t + step;
end if;
when 2 => dre := REAL_PART(target_direction);
dim := IMAG_PART(target_direction);
alfa := ARCTAN(dim/dre);
x := h * COS(alfa - PI/4.0);
y := h * SIN(alfa - PI/4.0);
step := Create(x,y);
nt := t + step;
absnt := AbsVal(nt);
if (absnt > AbsVal(target))
or else (AbsVal(nt - target) < distance)
then step := Create(0.0,-y) * Create(h);
nt := t + step;
end if;
end case;
hh := AbsVal(step);
end if;
t := nt;
end Complex_Predictor;
procedure Circular_Predictor
( t : in out Complex_Number; theta : in out double_float;
t0_min_target,target : in Complex_Number;
h : in double_float ) is
e_i_theta : Complex_Number;
begin
theta := theta + h;
e_i_theta := Create(COS(theta),SIN(theta));
t := target + t0_min_target * e_i_theta;
end Circular_Predictor;
procedure Geometric_Predictor
( t : in out Complex_Number; target : in Complex_Number;
h,tol : in double_float ) is
nt : Complex_Number;
begin
nt := target - Create(h)*(target - t);
if abs( REAL_PART(nt) - REAL_PART(target) ) <= tol
then t := target;
else t := nt;
end if;
end Geometric_Predictor;
-- PREDICTORS FOR x :
procedure Secant_Predictor
( x : in out Solution_Array; prev_x : in Solution_Array;
fac : in Complex_Number; dist_x : in double_float ) is
j : natural;
xx : Solution_Array(x'range);
begin
Copy(x,xx);
for i in x'range loop
x(i).v := x(i).v + fac * ( x(i).v - prev_x(i).v );
j := xx'first;
Equals(xx,x(i).v,i,dist_x,j);
if j /= i
then Copy(xx,x);
exit;
end if;
end loop;
Clear(xx);
end Secant_Predictor;
-- PREDICTORS FOR t AND x :
procedure Secant_Single_Real_Predictor
( x : in out Vector; prev_x : in Vector;
t : in out Complex_Number; prev_t,target : in Complex_Number;
h,tol : in double_float; pow : in positive := 1 ) is
hh,tmp : double_float;
factor : Complex_Number;
begin
tmp := AbsVal(t - prev_t);
Real_Predictor(t,target,h,tol,pow,hh);
if tmp > tol
then factor := Create(hh/tmp);
x := x + factor * ( x - prev_x );
end if;
end Secant_Single_Real_Predictor;
procedure Secant_Multiple_Real_Predictor
( x : in out Solution_Array; prev_x : in Solution_Array;
t : in out Complex_Number; prev_t,target : in Complex_Number;
h,tol,dist_x : in double_float; pow : in positive := 1 ) is
hh,tmp : double_float;
factor : Complex_Number;
begin
tmp := AbsVal(t - prev_t);
Real_Predictor(t,target,h,tol,pow,hh);
if tmp > tol
then factor := Create(hh/tmp);
Secant_Predictor(x,prev_x,factor,dist_x);
end if;
for k in x'range loop
x(k).t := t;
end loop;
end Secant_Multiple_Real_Predictor;
procedure Secant_Single_Complex_Predictor
( x : in out Vector; prev_x : in Vector;
t : in out Complex_Number; prev_t,target : in Complex_Number;
h,tol,dist_t : in double_float; trial : in natural ) is
hh,tmp : double_float;
factor : Complex_Number;
begin
tmp := AbsVal(t - prev_t);
Complex_Predictor(t,target,h,tol,hh,dist_t,trial);
if tmp > tol
then factor := Create(hh/tmp);
x := x + factor * ( x - prev_x );
end if;
end Secant_Single_Complex_Predictor;
procedure Secant_Multiple_Complex_Predictor
( x : in out Solution_Array; prev_x : in Solution_Array;
t : in out Complex_Number; prev_t,target : in Complex_Number;
h,tol,dist_x,dist_t : in double_float; trial : in natural ) is
hh,tmp : double_float;
factor : Complex_Number;
begin
tmp := AbsVal(t - prev_t);
Complex_Predictor(t,target,h,tol,hh,dist_t,trial);
if tmp > tol
then factor := Create(hh/tmp);
Secant_Predictor(x,prev_x,factor,dist_x);
end if;
for k in x'range loop
x(k).t := t;
end loop;
end Secant_Multiple_Complex_Predictor;
procedure Secant_Circular_Predictor
( x : in out Vector; prev_x : in Vector;
t : in out Complex_Number; theta : in out double_float;
prev_t,t0_min_target,target : in Complex_Number;
h,tol : in double_float ) is
factor : Complex_Number;
tmp : double_float := AbsVal(t-prev_t);
begin
if tmp <= tol
then Circular_Predictor(t,theta,t0_min_target,target,h);
else factor := Create(h/tmp);
Circular_Predictor(t,theta,t0_min_target,target,h);
x := x + factor * ( x - prev_x );
end if;
end Secant_Circular_Predictor;
procedure Secant_Geometric_Predictor
( x : in out Vector; prev_x : in Vector;
t : in out Complex_Number; prev_t,target : in Complex_Number;
h,tol : in double_float ) is
dist_prev,dist : double_float;
factor,tmp : Complex_Number;
begin
dist_prev := AbsVal(t - prev_t); -- old stepsize
tmp := t;
Geometric_Predictor(t,target,h,tol);
dist := AbsVal(t - tmp); -- new stepsize
if dist_prev > tol
then factor := Create(dist/dist_prev);
x := x + factor * ( x - prev_x );
end if;
end Secant_Geometric_Predictor;
procedure Tangent_Single_Real_Predictor
( x : in out Vector; t : in out Complex_Number;
target : in Complex_Number; h,tol : in double_float;
pow : in positive := 1 ) is
n : natural := x'last;
info : integer;
norm_tan,hh : double_float;
rhs : Vector(x'range);
j : Matrix(1..n,1..n);
ipvt : Standard_Natural_Vectors.Vector(1..n);
prev_t : Complex_Number := t;
begin
Real_Predictor(t,target,h,tol,pow,hh);
j := dH(x,prev_t);
lufac(j,n,ipvt,info);
if info = 0
then rhs := dH(x,prev_t);
Min(rhs);
lusolve(j,n,ipvt,rhs);
norm_tan := Norm(rhs);
if norm_tan > tol
then hh := hh / norm_tan;
Mul(rhs,Create(hh));
Add(x,rhs);
end if;
end if;
end Tangent_Single_Real_Predictor;
procedure Tangent_Multiple_Real_Predictor
( x : in out Solution_Array; t : in out Complex_Number;
target : in Complex_Number; h,tol,dist_x : in double_float;
nsys : in out natural; pow : in positive := 1 ) is
n : natural := x(x'first).n;
norm_tan,hh : double_float;
rhs : Vector(1..n);
info : integer;
j : Matrix(1..n,1..n);
ipvt : Standard_Natural_Vectors.Vector(1..n);
xx : Solution_Array(x'range);
jj : natural;
prev_t : Complex_Number := t;
begin
Real_Predictor(t,target,h,tol,pow,hh);
Copy(x,xx);
for i in x'range loop
j := dH(x(i).v,prev_t);
lufac(j,n,ipvt,info);
if info = 0
then rhs := dH(x(i).v,prev_t);
Min(rhs);
lusolve(j,n,ipvt,rhs);
nsys := nsys + 1;
norm_tan := Norm(rhs);
if norm_tan > tol
then hh := hh / norm_tan;
Mul(rhs,Create(hh));
Add(x(i).v,rhs);
end if;
jj := xx'first;
Equals(xx,x(i).v,i,dist_x,jj);
if jj /= i
then Copy(xx,x); exit;
end if;
else Copy(xx,x); exit;
end if;
end loop;
Clear(xx);
for k in x'range loop
x(k).t := t;
end loop;
end Tangent_Multiple_Real_Predictor;
procedure Tangent_Single_Complex_Predictor
( x : in out Vector; t : in out Complex_Number;
target : in Complex_Number;
h,tol,dist_t : in double_float; trial : in natural) is
n : natural := x'last;
info : integer;
norm_tan,hh : double_float;
rhs : Vector(x'range);
j : Matrix(1..n,1..n);
ipvt : Standard_Natural_Vectors.Vector(1..n);
prev_t : Complex_Number := t;
begin
Complex_Predictor(t,target,h,tol,hh,dist_t,trial);
j := dH(x,prev_t);
lufac(j,n,ipvt,info);
if info = 0
then rhs := dH(x,prev_t);
Min(rhs);
lusolve(j,n,ipvt,rhs);
norm_tan := Norm(rhs);
if norm_tan > tol
then hh := hh / norm_tan;
Mul(rhs,Create(hh));
Add(x,rhs);
end if;
end if;
end Tangent_Single_Complex_Predictor;
procedure Tangent_Multiple_Complex_Predictor
( x : in out Solution_Array; t : in out Complex_Number;
target : in Complex_Number;
h,tol,dist_x,dist_t : in double_float;
trial : in natural; nsys : in out natural) is
n : natural := x(x'first).n;
norm_tan,hh : double_float;
rhs : Vector(1..n);
info : integer;
j : Matrix(1..n,1..n);
ipvt : Standard_Natural_Vectors.Vector(1..n);
xx : Solution_Array(x'range);
jj : natural;
prev_t : Complex_Number := t;
begin
Complex_Predictor(t,target,h,tol,hh,dist_t,trial);
Copy(x,xx);
for i in x'range loop
j := dH(x(i).v,prev_t);
lufac(j,n,ipvt,info);
if info = 0
then rhs := dH(x(i).v,prev_t);
Min(rhs);
lusolve(j,n,ipvt,rhs);
nsys := nsys + 1;
norm_tan := Norm(rhs);
if norm_tan > tol
then hh := hh / norm_tan;
Mul(rhs,Create(hh));
Add(x(i).v,rhs);
end if;
jj := xx'first;
Equals(xx,x(i).v,i,dist_x,jj);
if jj /= i
then Copy(xx,x); exit;
end if;
else Copy(xx,x); exit;
end if;
end loop;
Clear(xx);
for k in x'range loop
x(k).t := t;
end loop;
end Tangent_Multiple_Complex_Predictor;
procedure Tangent_Circular_Predictor
( x : in out Vector; t : in out Complex_Number;
theta : in out double_float;
t0_min_target,target : in Complex_Number;
h,tol : in double_float ) is
n : natural := x'last;
info : integer;
norm_tan : double_float;
rhs : Vector(x'range);
j : Matrix(1..n,1..n);
ipvt : Standard_Natural_Vectors.Vector(1..n);
begin
lufac(j,n,ipvt,info);
if info = 0
then rhs := dH(x,t);
Min(rhs);
lusolve(j,n,ipvt,rhs);
end if;
Circular_Predictor(t,theta,t0_min_target,target,h);
if info = 0
then norm_tan := Norm(rhs);
if norm_tan > tol
then Mul(rhs,Create(h/norm_tan));
Add(x,rhs);
end if;
end if;
end Tangent_Circular_Predictor;
procedure Tangent_Geometric_Predictor
( x : in out Vector; t : in out Complex_Number;
target : in Complex_Number; h,tol : in double_float ) is
n : natural := x'last;
info : integer;
norm_tan,step : double_float;
rhs : Vector(x'range);
j : Matrix(1..n,1..n);
ipvt : Standard_Natural_Vectors.Vector(1..n);
prev_t : Complex_Number := t;
begin
Geometric_Predictor(t,target,h,tol);
j := dH(x,prev_t);
lufac(j,n,ipvt,info);
if info = 0
then rhs := dH(x,prev_t);
Min(rhs);
lusolve(j,n,ipvt,rhs);
norm_tan := Norm(rhs);
if norm_tan > tol
then step := AbsVal(t-prev_t) / norm_tan;
Mul(rhs,Create(step));
Add(x,rhs);
end if;
end if;
end Tangent_Geometric_Predictor;
function Hermite ( t0,t1,t,x0,x1,v0,v1 : Complex_Number )
return Complex_Number is
-- DESCRIPTION :
-- Returns the value of the third degree interpolating polynomial x(t),
-- such that x(t0) = x0, x(t1) = x1, x'(t0) = v0 and x'(t1) = v1.
-- REQUIRED : t0 /= t1.
-- IMPLEMENTATION :
-- x(t) = a3*t^3 + a2*t^2 + a1*t + a0,
-- The four interpolation conditions lead to a linear system in
-- the coefficients of x(t). This system is first solved explicitly
-- and then the polynomial x(t) is evaluated.
a0,a1,a2,a3,t10,v10 : Complex_Number;
begin
t10 := t1 - t0;
v10 := (x1 - x0)/t10;
a3 := (v1 + v0 - Create(2.0)*v10)/(t10**2);
a2 := (v10 - v0 - (t1**2 + t1*t0 - Create(2.0)*t0**2)*a3)/t10;
a1 := v0 - (Create(3.0)*a3*t0 + Create(2.0)*a2)*t0;
a0 := x0 - ((a3*t0 + a2)*t0 + a1)*t0;
return (((a3*t + a2)*t + a1)*t + a0);
end Hermite;
function Hermite ( t0,t1,t : Complex_Number; x0,x1,v0,v1 : Vector )
return Vector is
-- DESCRIPTION :
-- Returns the value of the third degree interpolating polynomial x(t),
-- such that x(t0) = x0, x(t1) = x1, x'(t0) = v0 and x'(t1) = v1,
-- for every component.
-- REQUIRED : t0 /= t1.
res : Vector(x0'range);
begin
for i in res'range loop
res(i) := Hermite(t0,t1,t,x0(i),x1(i),v0(i),v1(i));
end loop;
return res;
end Hermite;
procedure Hermite_Single_Real_Predictor
( x : in out Vector; prev_x : in Vector;
t : in out Complex_Number; prev_t,target : in Complex_Number;
v : in out Vector; prev_v : in Vector;
h,tol : in double_float; pow : in positive := 1 ) is
n : natural := x'last;
info : integer;
hh : double_float;
j : Matrix(1..n,1..n);
ipvt : Standard_Natural_Vectors.Vector(1..n);
t1 : Complex_Number := t;
begin
Real_Predictor(t,target,h,tol,pow,hh);
j := dH(x,t1);
lufac(j,n,ipvt,info);
if info = 0
then v := dH(x,t1);
Min(v);
lusolve(j,n,ipvt,v);
if AbsVal(prev_t - t1) > tol
then x := Hermite(prev_t,t1,t,prev_x,x,prev_v,v);
end if;
end if;
end Hermite_Single_Real_Predictor;
end Predictors;