Return to multprec_root_refiners.adb CVS log

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

Import the second public release of PHCpack.

OKed by Jan Verschelde.

with integer_io;                         use integer_io;
with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Multprec_Floating_Numbers_io;       use Multprec_Floating_Numbers_io;
with Multprec_Complex_Numbers;           use Multprec_Complex_Numbers;
with Multprec_Complex_Numbers_io;        use Multprec_Complex_Numbers_io;
with Multprec_Complex_Vectors;           use Multprec_Complex_Vectors;
with Standard_Natural_Vectors;
with Multprec_Complex_Matrices;          use Multprec_Complex_Matrices;
with Multprec_Complex_Linear_Solvers;    use Multprec_Complex_Linear_Solvers;
with Multprec_Complex_Norms_Equals;      use Multprec_Complex_Norms_Equals;
with Multprec_Complex_Solutions_io;      use Multprec_Complex_Solutions_io;

package body Multprec_Root_Refiners is

-- AUXILIARIES :

  function Is_Real ( sol : Solution; tol : Floating_Number ) return boolean is

    res : boolean := true;

  begin
    for i in sol.v'range loop
      declare
        abstmp : Floating_Number := AbsVal(IMAG_PART(sol.v(i)));
      begin
        if abstmp > tol
         then res := false;
        end if;
        Clear(abstmp);
      end;
      exit when not res;
    end loop;
    return res;
  end Is_Real;

  function Equal ( s1,s2 : Solution; tol : Floating_Number ) return boolean is
  begin
    for i in s1.v'range loop
      if not Equal(s1.v(i),s2.v(i),tol)
       then return false;
      end if;
    end loop;
    return true;
  end Equal;

  function Complex_Conjugate ( s : Solution ) return Solution is

    res : Solution(s.n) := s;

  begin
    for i in res.v'range loop
      res.v(i) := Create(REAL_PART(s.v(i)),-IMAG_PART(s.v(i)));
    end loop;
    return res;
  end Complex_Conjugate;

  function Is_Clustered ( sol : Solution; nb : natural; sols : Solution_List;
                          tol : Floating_Number ) return natural is

    tmp : Solution_List := sols;
    cnt : natural := 0;

  begin
    while not Is_Null(tmp) loop
      cnt := cnt + 1;
      if cnt /= nb
       then if Equal(sol,Head_Of(tmp).all,tol)
	     then return cnt;
            end if;
      end if;
      tmp := Tail_Of(tmp);
    end loop;
    return nb;
  end Is_Clustered;

  function Is_Clustered ( sol : Solution; nb : natural; sols : Solution_Array;
                          tol : Floating_Number ) return natural is
  begin
    for i in sols'range loop
      if i /= nb
       then if Equal(sol,sols(i).all,tol)
             then return i;
            end if;
      end if;
    end loop;
    return nb;
  end Is_Clustered;

  function Multiplicity ( sol : Solution; sols : Solution_List; 
                          tol : Floating_Number ) return natural is

    tmp : Solution_List := sols;
    cnt : natural := 0;

  begin
    while not Is_Null(tmp) loop
      if Equal(sol,Head_Of(tmp).all,tol)
       then cnt := cnt + 1;
      end if;
      tmp := Tail_Of(tmp);
    end loop;
    return cnt;
  end Multiplicity;

  function Multiplicity ( sol : Solution; sols : Solution_Array;
                          tol : Floating_Number ) return natural is
    cnt : natural := 0;

  begin
    for i in sols'range loop
      if Equal(sol,sols(i).all,tol)
       then cnt := cnt + 1;
      end if;
    end loop;
    return cnt;
  end Multiplicity;

  procedure Write_Bar ( file : in file_type ) is
  begin
    for i in 1..75 loop
      put(file,'=');
    end loop;
    new_line(file);
  end Write_Bar;

  procedure Write_Info ( file : in file_type; zero : in Solution;
                         i,n,numb : in natural; fail : in boolean ) is

  -- DESCRIPTION :
  --   The information concerning the zero is written

  begin
   -- Write_Bar(file);
    put(file,"solution "); put(file,i,1); put(file," :        ");
    put(file," start residual : "); put(file,zero.res,3); new_line(file);
    put(file,zero);
  end Write_Info;
 
  procedure Root_Accounting
               ( file : in file_type; ls : in out Link_to_Solution;
                 nb : in natural; sa : in out Solution_Array;
                 fail : in boolean; tolsing,tolclus : in Floating_Number;
                 nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus : in out natural ) is

  -- DESCRIPTION :
  --   This procedure does root accounting of the solution sol, w.r.t. the
  --   solution list sols.  Information will be provided concerning the type
  --   of solution.

    tolreal : Floating_Number := Create(10.0**(-13));

  begin
    if fail
     then put_line(file," no solution ==");
          nbfail := nbfail + 1;
          ls.m := 0;
     else if Is_Real(ls.all,tolreal)
	   then put(file," real ");
	        nbreal := nbreal + 1;
	   else put(file," complex ");
	        nbcomp := nbcomp + 1;
          end if;
          if sa(nb).rco < tolsing
	   then declare
                  m : natural := Multiplicity(ls.all,sa,tolclus);
                begin
                  if m = 1
                   then m := 0;
		  end if;
		  ls.m := m;
	        end;
	        put_line(file,"singular ==");
	        nbsing := nbsing + 1;
           else declare
                  nb2 : natural := Is_Clustered(ls.all,nb,sa,tolclus);
                begin
                  if nb2 = nb
		   then put_line(file,"regular ==");
		        nbreg := nbreg + 1;
                   else put(file,"clustered : ");
		        put(file,nb2,1);
		        put_line(file," ==");
		        nbclus := nbclus + 1;
                  end if;
		  ls.m := 1;
	        end;	   
          end if;  
    end if;
    Clear(tolreal);
  end Root_Accounting;

  procedure Root_Accounting 
                ( ls : in out Link_to_Solution; nb : in natural;
                  sa : in out Solution_Array; fail : in boolean;
                  tolsing,tolclus : in Floating_Number ) is

  -- DESCRIPTION :
  --   This procedure does root accounting of the solution sol, w.r.t. the
  --   solution list sols.  Information will be provided concerning the type
  --   of solution.

  begin
    if fail
     then ls.m := 0;
     elsif sa(nb).rco < tolsing
         then declare
                m : natural := Multiplicity(ls.all,sa,tolclus);
              begin
                if m = 1
                 then ls.m := 0;
                 else ls.m := m;
                end if;
              end;
         else ls.m := 1;
    end if;
  end Root_Accounting;

  procedure Write_Global_Info
             ( file : in file_type;
               tot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus : in natural ) is

  begin
    Write_Bar(file);
    put(file,"A list of "); put(file,tot,1);
    put_line(file," solutions has been refined :");
    put(file,"Number of regular solutions   : "); put(file,nbreg,1);
    put_line(file,".");
    put(file,"Number of singular solutions  : "); put(file,nbsing,1);
    put_line(file,".");
    put(file,"Number of real solutions      : "); put(file,nbreal,1);
    put_line(file,".");
    put(file,"Number of complex solutions   : "); put(file,nbcomp,1);
    put_line(file,".");
    put(file,"Number of clustered solutions : "); put(file,nbclus,1);
    put_line(file,".");
    put(file,"Number of failures            : "); put(file,nbfail,1);
    put_line(file,".");
    Write_Bar(file);
  end Write_Global_Info;

-- TARGET ROUTINES :

  procedure Silent_Newton
               ( p_eval : in Eval_Poly_Sys; j_eval : in Eval_Jaco_Mat;
                 zero : in out Solution; epsxa,epsfa : in Floating_Number; 
                 numit : in out natural; max : in natural;
                 fail : out boolean ) is

    n : natural := p_eval'length;
    jacobian : matrix(1..n,1..n);
    ipvt : Standard_Natural_Vectors.Vector(1..n);
    y,deltax : Vector(1..n);

  begin
    y := eval(p_eval,zero.v);               -- y = f(zero)
    for i in 1..max loop
      jacobian := Eval(j_eval,zero.v);      -- solve jacobian*deltax = -f(zero)
      lufco(jacobian,n,ipvt,zero.rco);
      if Equal(zero.rco,0.0)
       then fail := (Sum_Norm(y) > epsfa);
            return;                         -- singular Jacobian matrix
      end if;
      deltax := -y; Clear(y);
      lusolve(jacobian,n,ipvt,deltax);  
      Add(zero.v,deltax);                   -- make the updates
      y := eval(p_eval,zero.v);
      Clear(zero.err);              Clear(zero.res);
      zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
      numit := numit + 1;
      if ( zero.err < epsxa ) or else ( zero.res < epsfa ) 
                                            -- stopping criteria
       then fail := false; exit;
       elsif numit >= max
           then fail := true; exit;
      end if;
      Clear(deltax); Clear(jacobian);
    end loop;
    jacobian := eval(j_eval,zero.v);        -- compute condition number
    lufco(jacobian,n,ipvt,zero.rco);
    Clear(jacobian);
  exception
    when numeric_error | constraint_error => fail := true; return;
  end Silent_Newton;

  procedure Reporting_Newton
               ( file : in file_type;
                 p_eval : in Eval_Poly_Sys; j_eval : in Eval_Jaco_Mat;
                 zero : in out Solution; epsxa,epsfa : in Floating_Number; 
                 numit : in out natural; max : in natural;
                 fail : out boolean ) is

    n : natural := p_eval'length;
    jacobian : matrix(1..n,1..n);
    ipvt : Standard_Natural_Vectors.Vector(1..n);
    y,deltax : Vector(1..n);

  begin
    y := Eval(p_eval,zero.v);              -- y = f(zero)
    for i in 1..max loop
      jacobian := eval(j_eval,zero.v);     -- solve jacobian*deltax = -f(zero)
      lufco(jacobian,n,ipvt,zero.rco);
      if Equal(zero.rco,0.0)               -- singular jacobian matrix
       then fail := (Sum_Norm(y) > epsfa); -- accuracy not reached yet
            return;
      end if;
      deltax := -y; Clear(y);
      lusolve(jacobian,n,ipvt,deltax);  
      Add(zero.v,deltax);                  -- make the updates
      y := eval(p_eval,zero.v);
      Clear(zero.err);              Clear(zero.res);
      zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
      numit := numit + 1;
      put(file,"Step "); put(file,numit,4); new_line(file);      -- output
      put(file," |errxa| : "); put(file,zero.err); new_line(file);
      put(file," |errfa| : "); put(file,zero.res); new_line(file);
      if ( zero.err < epsxa ) or else ( zero.res < epsfa ) 
                                                  -- stopping criteria
       then fail := false; exit;
       elsif numit >= max
           then fail := true; exit;
      end if;
      Clear(deltax); Clear(jacobian);
    end loop;
    jacobian := eval(j_eval,zero.v);              -- compute condition number
    lufco(jacobian,n,ipvt,zero.rco);
    Clear(jacobian);
  exception
    when numeric_error | constraint_error => fail := true; return;
  end Reporting_Newton;

  procedure Silent_Root_Refiner
               ( p : in Poly_Sys; sols : in out Solution_List;
                 epsxa,epsfa,tolsing : in Floating_Number;
                 numit : in out natural; max : in natural ) is

    n : natural := p'length;
    p_eval : Eval_Poly_Sys(1..n) := Create(p);
    jac : Jaco_Mat(1..n,1..n) := Create(p);
    jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
    numb : natural;
    fail : boolean;
    sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
    eval_acc : Vector(1..n);

  begin
    for i in sa'range loop
      numb := 0;
      Clear(sa(i).res);
      eval_acc := Eval(p_eval,sa(i).v);
      sa(i).res := Sum_Norm(eval_acc);
      Clear(eval_acc);
      if sa(i).res < 1.0
       then Silent_Newton(p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
       else fail := true;
      end if;
      Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
      numit := numit + numb;
    end loop;
    clear(p_eval); clear(jac); clear(jac_eval);
    Shallow_Clear(sols); sols := Create(sa); -- Shallow_Clear(sa);
  end Silent_Root_Refiner;

  procedure Silent_Root_Refiner
               ( p : in Poly_Sys; sols,refsols : in out Solution_List;
                 epsxa,epsfa,tolsing : in Floating_Number;
                 numit : in out natural; max : in natural ) is

    n : natural := p'length;
    p_eval : Eval_Poly_Sys(1..n) := Create(p);
    jac : Jaco_Mat(1..n,1..n) := Create(p);
    jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
    numb : natural;
    fail : boolean;
    sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
    refsols_last : Solution_List;
    eval_acc : Vector(1..n);

  begin
    for i in sa'range loop
      numb := 0;
      Clear(sa(i).res);
      eval_acc := Eval(p_eval,sa(i).v);
      sa(i).res := Sum_Norm(eval_acc);
      Clear(eval_acc);
      if sa(i).res < 1.0
       then Silent_Newton(p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
            if not fail
             then Append(refsols,refsols_last,sa(i).all);
            end if;
       else fail := true;
      end if;
      Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
      numit := numit + numb;
    end loop;
    clear(p_eval); clear(jac); clear(jac_eval);
    Shallow_Clear(sols); sols := Create(sa); -- Shallow_Clear(sa);
  end Silent_Root_Refiner;

  procedure Reporting_Root_Refiner
               ( file : in file_type;
                 p : in Poly_Sys; sols : in out Solution_List;
                 epsxa,epsfa,tolsing : in Floating_Number;
                 numit : in out natural; max : in natural;
                 wout : in boolean ) is

    n : natural := p'length;
    p_eval : Eval_Poly_Sys(1..n) := Create(p);
    jac : Jaco_Mat(1..n,1..n) := Create(p);
    jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
    numb : natural;
    nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
    nbtot : constant natural := Length_Of(sols);
    fail : boolean;
    sa : Solution_Array(1..nbtot) := Create(sols);
    eval_acc : Vector(1..n);

  begin
    new_line(file);
    put_line(file,"THE SOLUTIONS :"); new_line(file);
    put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
    Write_Bar(file);
    for i in sa'range loop 
      numb := 0;
      Clear(sa(i).res);
      eval_acc := Eval(p_eval,sa(i).v);
      sa(i).res := Sum_Norm(eval_acc);
      Clear(eval_acc);
      if sa(i).res < 1.0
       then
         if wout
          then Reporting_Newton
                 (file,p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
          else Silent_Newton
                 (p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
         end if;
       else 
         fail := true;
      end if;
      Write_Info(file,sa(i).all,i,n,numb,fail);
      Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
                      nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
      numit := numit + numb;
    end loop;
    Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
    clear(p_eval); clear(jac); clear(jac_eval);
    Shallow_Clear(sols); sols := Create(sa); -- Shallow_Clear(sa);
  end Reporting_Root_Refiner;

  procedure Reporting_Root_Refiner
               ( file : in file_type;
                 p : in Poly_Sys; sols,refsols : in out Solution_List;
                 epsxa,epsfa,tolsing : in Floating_Number;
                 numit : in out natural; max : in natural;
                 wout : in boolean ) is

    n : natural := p'length;
    p_eval : Eval_Poly_Sys(1..n) := Create(p);
    jac : Jaco_Mat(1..n,1..n) := Create(p);
    jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
    numb : natural;
    nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
    nbtot : constant natural := Length_Of(sols);
    fail : boolean;
    sa : Solution_Array(1..nbtot) := Create(sols);
    refsols_last : Solution_List;
    eval_acc : Vector(1..n);

  begin
    new_line(file);
    put_line(file,"THE SOLUTIONS :"); new_line(file);
    put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
    Write_Bar(file);
    for i in sa'range loop 
      numb := 0;
      Clear(sa(i).res);
      eval_acc := Eval(p_eval,sa(i).v);
      sa(i).res := Sum_Norm(eval_acc);
      Clear(eval_acc);
      if sa(i).res < 1.0
       then
         if wout
          then Reporting_Newton
                 (file,p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
          else Silent_Newton
                 (p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
         end if;
         if not fail
          then Append(refsols,refsols_last,sa(i).all);
         end if;
       else 
         fail := true;
      end if;
      Write_Info(file,sa(i).all,i,n,numb,fail);
      Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
                      nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
      numit := numit + numb;
    end loop;
    Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
    clear(p_eval); clear(jac); clear(jac_eval);
    Shallow_Clear(sols); sols := Create(sa); -- Shallow_Clear(sa);
  end Reporting_Root_Refiner;

end Multprec_Root_Refiners;