Return to predictors.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Continuation

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;