File: [local] / OpenXM_contrib / PHC / Ada / Continuation / drivers_for_poly_continuation.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 integer_io; use integer_io;
with Communications_with_User; use Communications_with_User;
with Timing_Package; use Timing_Package;
with File_Scanning; use File_Scanning;
with Standard_Floating_Numbers; use Standard_Floating_Numbers;
with Standard_Complex_Numbers; use Standard_Complex_Numbers;
with Standard_Complex_Numbers_io; use Standard_Complex_Numbers_io;
with Numbers_io; use Numbers_io;
with Standard_Floating_Vectors; use Standard_Floating_Vectors;
with Standard_Floating_VecVecs; use Standard_Floating_VecVecs;
with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
with Symbol_Table,Symbol_Table_io; use Symbol_Table;
with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
with Standard_Complex_Solutions_io; use Standard_Complex_Solutions_io;
with Homotopy;
with Projective_Transformations; use Projective_Transformations;
with Drivers_for_Homotopy_Creation; use Drivers_for_Homotopy_Creation;
with Continuation_Parameters;
with Continuation_Parameters_io;
with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
with Process_io; use Process_io;
with Drivers_for_Path_Directions; use Drivers_for_Path_Directions;
package body Drivers_for_Poly_Continuation is
-- AUXILIARIES :
procedure Continue ( file : in file_type; sols : in out Solution_List;
proj,report : in boolean;
target : in Complex_Number ) is
-- DESCRIPTION :
-- Instantiates the path-trackers.
timer : timing_widget;
procedure Sil_Cont is
new Silent_Continue(Max_Norm,
Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
procedure Rep_Cont is
new Reporting_Continue(Max_Norm,
Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
begin
tstart(timer);
if report
then Rep_Cont(file,sols,proj,target);
else Sil_Cont(sols,proj,target);
end if;
tstop(timer);
new_line(file); print_times(file,timer,"continuation");
end Continue;
procedure Ask_Symbol is
-- DESCRIPTION :
-- This procedure asks for the symbol to display the additional unknown.
sb : Symbol;
begin
put("Give symbol to display additional unknown : ");
sb := (sb'range => ' ');
Symbol_Table.Enlarge(1);
Symbol_Table_io.Get(sb);
Symbol_Table.Add(sb);
end Ask_Symbol;
-- TARGET ROUTINES :
procedure Driver_for_Process_io ( file : in file_type; oc : out natural ) is
ans : character;
m : array(0..8) of string(1..65);
begin
put_line("MENU for Output Information during Continuation : ");
m(0):=" 0 : no intermediate output information during continuation ";
m(1):=" 1 : only the final solutions at the end of the paths ";
m(2):=" 2 : intermediate solutions at each step along the paths ";
m(3):=" 3 : information of the predictor: t and step length ";
m(4):=" 4 : information of the corrector: corrections and residuals ";
m(5):=" 5 : intermediate solutions and information of the predictor ";
m(6):=" 6 : intermediate solutions and information of the corrector ";
m(7):=" 7 : information of predictor and corrector ";
m(8):=" 8 : intermediate solutions, info of predictor and corrector ";
for i in m'range loop
put_line(m(i));
end loop;
put("Type a number between 0 and 8 to select output information : ");
Ask_Alternative(ans,"012345678");
new_line(file);
put_line(file,"OUTPUT INFORMATION DURING CONTINUATION :");
case ans is
when '0' => Set_output_code(nil); oc := 0; put_line(file,m(0));
when '1' => Set_output_code(nil); oc := 1; put_line(file,m(1));
when '2' => Set_output_code(s); oc := 2; put_line(file,m(2));
when '3' => Set_output_code(p); oc := 3; put_line(file,m(3));
when '4' => Set_output_code(c); oc := 4; put_line(file,m(4));
when '5' => Set_output_code(sp); oc := 5; put_line(file,m(5));
when '6' => Set_output_code(sc); oc := 6; put_line(file,m(6));
when '7' => Set_output_code(pc); oc := 7; put_line(file,m(7));
when '8' => Set_output_code(spc); oc := 8; put_line(file,m(8));
when others => null;
end case;
end Driver_for_Process_io;
procedure Driver_for_Continuation_Parameters ( file : in file_type ) is
nb : natural := 0;
procedure Begin_Banner ( file : in file_type ) is
begin
put_line(file,"****************** CURRENT CONTINUATION PARAMETERS "
& "*****************");
end Begin_Banner;
procedure End_Banner ( file : in file_type ) is
begin
put_line(file,"***************************************************"
& "*****************");
end End_Banner;
begin
loop
Begin_Banner(Standard_Output);
Continuation_Parameters_io.put;
End_Banner(Standard_Output);
Continuation_Parameters_io.get(nb);
exit when (nb = 0);
end loop;
new_line(file);
Begin_Banner(file);
Continuation_Parameters_io.put(file);
End_Banner(file);
end Driver_for_Continuation_Parameters;
procedure Check_Continuation_Parameter
( sols : in out Solution_List ) is
ans : character;
tre,tim : double_float;
begin
if not Is_Null(sols)
then
if Head_Of(sols).t = Create(1.0)
then put_line("The first solution has continuation parameter t = 1.0.");
put("Do you want to change t ? (y/n) "); Ask_Yes_or_No(ans);
if ans = 'y'
then put("Give real part of t : "); Read_Double_Float(tre);
put("Give imaginary part of t : "); Read_Double_Float(tim);
Set_Continuation_Parameter(sols,Create(tre,tim));
end if;
end if;
end if;
end Check_Continuation_Parameter;
procedure Driver_for_Polynomial_Continuation
( file : in file_type; p : in Poly_Sys;
sols : out Solution_List; target : out Complex_Number ) is
infile : file_type;
pp,q : Poly_Sys(p'range);
t : Complex_Number;
qsols : Solution_List;
found,proj : boolean;
procedure Read_Start_System is
firstfail : boolean := true;
n : natural;
begin
put_line("Reading the name of the file for start system.");
Read_Name_and_Open_File(infile);
get(infile,n,q);
exception
when others => put("The system on the file is not correct.");
if firstfail
then put_line(" Try again..."); Close(infile);
firstfail := false;
Read_Start_System;
end if;
end Read_Start_System;
begin
new_line;
Read_Start_System;
put_line(file,"THE START SYSTEM : "); put(file,q); new_line(file);
Scan_and_Skip(infile,"SOLUTIONS",found);
if found
then get(infile,qsols);
else new_line; Read(qsols);
end if;
Close(infile);
Check_Continuation_Parameter(qsols);
put_line(file,"THE START SOLUTIONS : ");
put(file,Length_Of(qsols),Head_Of(qsols).n,qsols); new_line(file);
Copy(p,pp);
Driver_for_Homotopy_Construction(file,pp,q,qsols,t);
proj := (Number_of_Unknowns(q(q'first)) > q'last);
if proj
then Ask_Symbol;
end if;
new_line;
Driver_for_Polynomial_Continuation(file,qsols,proj,t);
-- Homotopy.Clear; --> clearing here creates difficulties for root refiner
sols := qsols;
target := t;
end Driver_for_Polynomial_Continuation;
procedure Driver_for_Polynomial_Continuation
( file : in file_type; p : in Poly_Sys; k : in natural;
target : in Complex_Number; sols : out Solution_list ) is
qsols : Solution_List;
begin
new_line; Read(qsols);
put_line(file,"THE START SOLUTIONS :");
put(file,Length_Of(qsols),Head_Of(qsols).n,qsols); new_line(file);
Homotopy.Create(p,k);
put_line(file,"HOMOTOPY PARAMETERS :");
put(file," k : "); put(file,k,2); new_line(file);
put(file," a : "); put(file,target); new_line(file);
Driver_for_Polynomial_Continuation(file,qsols,false,target);
-- Homotopy.Clear; --> clearing here creates difficulties for root refiner
sols := qsols;
end Driver_for_Polynomial_Continuation;
procedure Driver_for_Polynomial_Continuation
( file : in file_type; sols : in out Solution_List;
proj : in boolean;
target : Complex_Number := Create(1.0) ) is
oc : natural;
timer : timing_widget;
report : boolean;
n : constant natural := Head_Of(sols).n;
nv : constant natural := Length_Of(sols);
v : Link_to_VecVec;
errv : Link_to_Vector;
begin
new_line;
Driver_for_Continuation_Parameters(file);
if Continuation_Parameters.endext_order > 0
then Init_Path_Directions(n,nv,v,errv);
end if;
new_line;
Driver_for_Process_io(file,oc);
report := (oc /= 0);
new_line;
put_line("No more input expected. See output file for results.");
new_line;
if Continuation_Parameters.endext_order > 0
then Toric_Continue(file,sols,proj,report,v.all,errv.all,target);
Write_Directions(file,v.all,errv.all);
else Continue(file,sols,proj,report,target);
end if;
end Driver_for_Polynomial_Continuation;
end Drivers_for_Poly_Continuation;