File: [local] / OpenXM_contrib / PHC / Ada / Continuation / driver_for_winding_numbers.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_Floating_Numbers_io; use Standard_Floating_Numbers_io;
with Standard_Complex_Numbers; use Standard_Complex_Numbers;
with Numbers_io; use Numbers_io;
with Standard_Complex_Vectors; use Standard_Complex_Vectors;
with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
with Homotopy;
with Standard_Complex_Solutions_io; use Standard_Complex_Solutions_io;
with Continuation_Parameters; use Continuation_Parameters;
with Continuation_Data; use Continuation_Data;
with Path_Trackers; use Path_Trackers;
with Drivers_for_Homotopy_Creation; use Drivers_for_Homotopy_Creation;
with Drivers_for_Poly_Continuation; use Drivers_for_Poly_Continuation;
procedure Driver_for_Winding_Numbers
( file : in file_type; p : in Poly_Sys;
sols : in out Solution_List ) is
infile : file_type;
timer : timing_widget;
pp,q : Poly_Sys(p'range);
found,proj : boolean := false;
target : Complex_Number;
ans : character;
oc,max_wc : natural;
procedure Write_Statistics_and_Condition
( file : in file_type; i,nstep,nfail,niter,nsyst : in natural;
rcond : in double_float ) is
-- DESCRIPTION :
-- Writes the computing statistics of the ith path on file, followed
-- by the estimate for the inverse of the condition number.
begin
put(file,"========================================");
put_line(file,"===================================");
put(file,"== "); put(file,i,1); put(file," = ");
put(file," #step : "); put(file,nstep,3);
put(file," #fail : "); put(file,nfail,2);
put(file," #iter : "); put(file,niter,3);
if nsyst /= niter
then put(file," #syst : "); put(file,nsyst,3);
end if;
put(file," = ");
put(file," rco : "); put(file,rcond,2,3,3);
put_line(file," ==");
end Write_Statistics_and_Condition;
procedure Winding_Numbers
( file : in file_type; sols : in out Solution_List;
proj : in boolean; target : in Complex_Number;
mwc : in natural ) is
-- DESCRIPTION :
-- Computes the winding number for the given list of solutions.
-- REQUIRED :
-- The homotopy is already defined and stored in the package homotopy.
sa : Solu_Info_Array(1..Length_Of(sols)) := Deep_Create(sols);
pen : Pred_Pars := Continuation_Parameters.Create_End_Game;
cen : Corr_Pars := Continuation_Parameters.Create_End_Game;
tol : constant double_float := 10.0**(-10);
epslop : constant double_float := 10.0**(-6);
wc : natural;
sum,allsum : Standard_Complex_Vectors.Vector(sa(sa'first).sol.v'range);
procedure CCont2 is
new Circular_Single_Conditioned_Reporting_Continue
(Max_Norm,Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
begin
for i in sa'range loop
if AbsVal(Create(1.0) - sa(i).sol.t) > epslop
then
CCont2(file,sa(i),target,tol,epslop,wc,mwc,sum,allsum,false,pen,cen);
sa(i).sol.m := wc;
sa(i).sol.v := allsum;
sa(i).sol.t := target;
Write_Statistics_and_Condition
(file,i,sa(i).nstep,sa(i).nfail,sa(i).niter,sa(i).nsyst,sa(i).rcond);
put(file,sa(i).sol.all);
end if;
end loop;
Clear(sols);
sols := Shallow_Create(sa);
end Winding_Numbers;
begin
-- READING MAXIMUM WINDING NUMBER :
new_line;
put("Give the maximum winding number : "); Read_Natural(max_wc);
-- READING THE START SYSTEM :
new_line;
put_line("Reading the name of the file for start system.");
Read_Name_and_Open_File(infile); get(infile,q);
new_line;
put_line(file,"THE START SYSTEM :");
new_line(file); put(file,q); new_line(file);
-- CONSTRUCTING THE HOMOTOPY AND TUNING OF PARAMETERS :
Copy(p,pp);
Driver_for_Homotopy_Construction(file,pp,q,sols,target);
proj := (Number_of_Unknowns(q(q'first)) > q'last);
Driver_for_Continuation_Parameters(file);
new_line;
put("Do you want intermediate output during continuation ? (y/n) ");
Ask_Yes_or_No(ans);
if ans = 'y'
then Driver_for_Process_io(file,oc);
end if;
-- COMPUTATION OF WINDING NUMBERS :
new_line;
put_line("No more input desired. Computing winding numbers ...");
put_line("The program can now run in the background.");
new_line;
tstart(timer);
Winding_Numbers(file,sols,proj,target,max_wc);
tstop(timer);
new_line(file);
print_times(file,timer,"computation of winding numbers");
new_line(file);
end Driver_for_Winding_Numbers;