Return to standard_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_io;       use Standard_Floating_Numbers_io;
with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Complex_Numbers_io;        use Standard_Complex_Numbers_io;
with Standard_Complex_Vectors;           use Standard_Complex_Vectors;
with Standard_Natural_Vectors;
with Standard_Complex_Matrices;          use Standard_Complex_Matrices;
with Standard_Complex_Linear_Solvers;    use Standard_Complex_Linear_Solvers;
with Standard_Complex_Norms_Equals;      use Standard_Complex_Norms_Equals;
with Standard_Complex_Solutions_io;      use Standard_Complex_Solutions_io;

package body Standard_Root_Refiners is

-- AUXILIARIES :

  function Is_Real ( sol : Solution; tol : double_float ) return boolean is
  begin
    for i in sol.v'range loop
      if ABS(IMAG_PART(sol.v(i))) > tol
       then return false;
      end if;
    end loop;
    return true;
  end Is_Real;

  function Equal ( s1,s2 : Solution; tol : double_float ) 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 : double_float ) 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 : double_float ) 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 : double_float ) 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 : double_float ) 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,2,3,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 double_float;
                 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.

  begin
    if fail
     then put_line(file," no solution ==");
          nbfail := nbfail + 1;
          ls.m := 0;
     else if Is_Real(ls.all,10.0**(-13))
	   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;
  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 double_float ) 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 double_float; 
                 numit : in out natural; max : in natural;
                 fail : out boolean ) is

    n : natural := zero.n;
    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 1.0 + zero.rco = 1.0
       then fail := (Sum_Norm(y) > epsfa);
            return;                         -- singular Jacobian matrix
      end if;
      deltax := -y;
      lusolve(jacobian,n,ipvt,deltax);  
      Add(zero.v,deltax);                   -- make the updates
      y := eval(p_eval,zero.v);
      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;
    end loop;
    jacobian := eval(j_eval,zero.v);        -- compute condition number
    lufco(jacobian,n,ipvt,zero.rco);
  exception
    when numeric_error | constraint_error => fail := true; return;
  end Silent_Newton;

  procedure Silent_Newton
               ( p_eval : in Standard_Complex_Poly_SysFun.Evaluator;
                 j_eval : in Standard_Complex_Jaco_Matrices.Evaluator;
                 zero : in out Solution; epsxa,epsfa : in double_float;
                 numit : in out natural; max : in natural;
                 fail : out boolean ) is

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

  begin
    y := p_eval(zero.v);                    -- y = f(zero)
    for i in 1..max loop
      jacobian := j_eval(zero.v);           -- solve jacobian*deltax = -f(zero)
      lufco(jacobian,n,ipvt,zero.rco);
      if 1.0 + zero.rco = 1.0
       then fail := (Sum_Norm(y) > epsfa);
            return;                         -- singular Jacobian matrix
      end if;
      deltax := -y;
      lusolve(jacobian,n,ipvt,deltax);
      Add(zero.v,deltax);                   -- make the updates
      y := p_eval(zero.v);
      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;
    end loop;
    jacobian := j_eval(zero.v);            -- compute condition number
    lufco(jacobian,n,ipvt,zero.rco);
  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 double_float; 
                 numit : in out natural; max : in natural;
                 fail : out boolean ) is

    n : natural := zero.n;
    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 1.0 + zero.rco = 1.0              -- singular jacobian matrix
       then fail := (Sum_Norm(y) > epsfa); -- accuracy not reached yet
            return;
      end if;
      deltax := -y;
      lusolve(jacobian,n,ipvt,deltax);  
      Add(zero.v,deltax);                  -- make the updates
      y := eval(p_eval,zero.v);
      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;
    end loop;
    jacobian := eval(j_eval,zero.v);              -- compute condition number
    lufco(jacobian,n,ipvt,zero.rco);
  exception
    when numeric_error | constraint_error => fail := true; return;
  end Reporting_Newton;

  procedure Reporting_Newton
               ( file : in file_type;
                 p_eval : in Standard_Complex_Poly_SysFun.Evaluator;
                 j_eval : in Standard_Complex_Jaco_Matrices.Evaluator;
                 zero : in out Solution; epsxa,epsfa : in double_float;
                 numit : in out natural; max : in natural;
                 fail : out boolean ) is

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

  begin
    y := p_eval(zero.v);                   -- y = f(zero)
    for i in 1..max loop
      jacobian := j_eval(zero.v);          -- solve jacobian*deltax = -f(zero)
      lufco(jacobian,n,ipvt,zero.rco);
      if 1.0 + zero.rco = 1.0              -- singular jacobian matrix
       then fail := (Sum_Norm(y) > epsfa); -- accuracy not reached yet
            return;
      end if;
      deltax := -y;
      lusolve(jacobian,n,ipvt,deltax);
      Add(zero.v,deltax);                  -- make the updates
      y := p_eval(zero.v);
      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;
    end loop;
    jacobian := j_eval(zero.v);                   -- compute condition number
    lufco(jacobian,n,ipvt,zero.rco);
  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 double_float;
                 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);

  begin
    for i in sa'range loop
      numb := 0;
      sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
      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);
    Deep_Clear(sols); sols := Create(sa); Clear(sa);
  end Silent_Root_Refiner;

  procedure Silent_Root_Refiner
               ( p : in Standard_Complex_Poly_SysFun.Evaluator;
                 j : in Standard_Complex_Jaco_Matrices.Evaluator;
                 sols : in out Solution_List;
                 epsxa,epsfa,tolsing : in double_float;
                 numit : in out natural; max : in natural ) is

    n : constant natural := Head_Of(sols).n;
    numb : natural;
    fail : boolean;
    sa : Solution_Array(1..Length_Of(sols)) := Create(sols);

  begin
    for i in sa'range loop
      numb := 0;
      sa(i).res := Sum_Norm(p(sa(i).v));
      if sa(i).res < 1.0
       then Silent_Newton(p,j,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;
    Deep_Clear(sols); sols := Create(sa); Clear(sa);
  end Silent_Root_Refiner;

  procedure Silent_Root_Refiner
               ( p : in Poly_Sys; sols,refsols : in out Solution_List;
                 epsxa,epsfa,tolsing : in double_float;
                 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;

  begin
    for i in sa'range loop
      numb := 0;
      sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
      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);
    Deep_Clear(sols); sols := Create(sa); Clear(sa);
  end Silent_Root_Refiner;

  procedure Silent_Root_Refiner
               ( p : in Standard_Complex_Poly_SysFun.Evaluator;
                 j : in Standard_Complex_Jaco_Matrices.Evaluator;
                 sols,refsols : in out Solution_List;
                 epsxa,epsfa,tolsing : in double_float;
                 numit : in out natural; max : in natural ) is

    n : constant natural := Head_Of(sols).n;
    numb : natural;
    fail : boolean;
    sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
    refsols_last : Solution_List;

  begin
    for i in sa'range loop
      numb := 0;
      sa(i).res := Sum_Norm(p(sa(i).v));
      if sa(i).res < 1.0
       then Silent_Newton(p,j,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;
    Deep_Clear(sols); sols := Create(sa); 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 double_float;
                 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);

  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;
      sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
      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);
    Deep_Clear(sols); sols := Create(sa); Clear(sa);
  end Reporting_Root_Refiner;

  procedure Reporting_Root_Refiner
               ( file : in file_type;
                 p : in Standard_Complex_Poly_SysFun.Evaluator;
                 j : in Standard_Complex_Jaco_Matrices.Evaluator;
                 sols : in out Solution_List;
                 epsxa,epsfa,tolsing : in double_float;
                 numit : in out natural; max : in natural;
                 wout : in boolean ) is

    n : constant natural := Head_Of(sols).n;
    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);

  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;
      sa(i).res := Sum_Norm(p(sa(i).v));
      if sa(i).res < 1.0
       then
         if wout
          then Reporting_Newton(file,p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
          else Silent_Newton(p,j,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);
    Deep_Clear(sols); sols := Create(sa); 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 double_float;
                 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;

  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;
      sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
      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);
    Deep_Clear(sols); sols := Create(sa); Clear(sa);
  end Reporting_Root_Refiner;

  procedure Reporting_Root_Refiner
               ( file : in file_type;
                 p : in Standard_Complex_Poly_SysFun.Evaluator;
                 j : in Standard_Complex_Jaco_Matrices.Evaluator;
                 sols,refsols : in out Solution_List;
                 epsxa,epsfa,tolsing : in double_float;
                 numit : in out natural; max : in natural;
                 wout : in boolean ) is

    n : constant natural := Head_Of(sols).n;
    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;

  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;
      sa(i).res := Sum_Norm(p(sa(i).v));
      if sa(i).res < 1.0
       then
         if wout
          then Reporting_Newton(file,p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
          else Silent_Newton(p,j,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);
    Deep_Clear(sols); sols := Create(sa); Clear(sa);
  end Reporting_Root_Refiner;

end Standard_Root_Refiners;