Return to driver_for_winding_numbers.adb CVS log

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

Annotation of OpenXM_contrib/PHC/Ada/Continuation/driver_for_winding_numbers.adb, Revision 1.1.1.1

1.1       maekawa     1: with integer_io;                         use integer_io;
                      2: with Communications_with_User;           use Communications_with_User;
                      3: with Timing_Package;                     use Timing_Package;
                      4: with File_Scanning;                      use File_Scanning;
                      5: with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
                      6: with Standard_Floating_Numbers_io;       use Standard_Floating_Numbers_io;
                      7: with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
                      8: with Numbers_io;                         use Numbers_io;
                      9: with Standard_Complex_Vectors;           use Standard_Complex_Vectors;
                     10: with Standard_Complex_Norms_Equals;      use Standard_Complex_Norms_Equals;
                     11: with Standard_Complex_Polynomials;       use Standard_Complex_Polynomials;
                     12: with Standard_Complex_Poly_Systems_io;   use Standard_Complex_Poly_Systems_io;
                     13: with Homotopy;
                     14: with Standard_Complex_Solutions_io;      use Standard_Complex_Solutions_io;
                     15: with Continuation_Parameters;            use Continuation_Parameters;
                     16: with Continuation_Data;                  use Continuation_Data;
                     17: with Path_Trackers;                      use Path_Trackers;
                     18: with Drivers_for_Homotopy_Creation;      use Drivers_for_Homotopy_Creation;
                     19: with Drivers_for_Poly_Continuation;      use Drivers_for_Poly_Continuation;
                     20:
                     21: procedure Driver_for_Winding_Numbers
                     22:              ( file : in file_type; p : in Poly_Sys;
                     23:                sols : in out Solution_List ) is
                     24:
                     25:   infile : file_type;
                     26:   timer : timing_widget;
                     27:   pp,q : Poly_Sys(p'range);
                     28:   found,proj : boolean := false;
                     29:   target : Complex_Number;
                     30:   ans : character;
                     31:   oc,max_wc : natural;
                     32:
                     33:   procedure Write_Statistics_and_Condition
                     34:                 ( file : in file_type; i,nstep,nfail,niter,nsyst : in natural;
                     35:                   rcond : in double_float ) is
                     36:
                     37:   -- DESCRIPTION :
                     38:   --   Writes the computing statistics of the ith path on file, followed
                     39:   --   by the estimate for the inverse of the condition number.
                     40:
                     41:   begin
                     42:     put(file,"========================================");
                     43:     put_line(file,"===================================");
                     44:     put(file,"== "); put(file,i,1); put(file," = ");
                     45:     put(file," #step : "); put(file,nstep,3);
                     46:     put(file," #fail : "); put(file,nfail,2);
                     47:     put(file," #iter : "); put(file,niter,3);
                     48:     if nsyst /= niter
                     49:      then put(file," #syst : "); put(file,nsyst,3);
                     50:     end if;
                     51:     put(file," = ");
                     52:     put(file," rco : "); put(file,rcond,2,3,3);
                     53:     put_line(file," ==");
                     54:   end Write_Statistics_and_Condition;
                     55:
                     56:   procedure Winding_Numbers
                     57:                ( file : in file_type; sols : in out Solution_List;
                     58:                  proj : in boolean; target : in Complex_Number;
                     59:                  mwc : in natural ) is
                     60:
                     61:   -- DESCRIPTION :
                     62:   --   Computes the winding number for the given list of solutions.
                     63:
                     64:   -- REQUIRED :
                     65:   --   The homotopy is already defined and stored in the package homotopy.
                     66:
                     67:     sa : Solu_Info_Array(1..Length_Of(sols)) := Deep_Create(sols);
                     68:     pen : Pred_Pars := Continuation_Parameters.Create_End_Game;
                     69:     cen : Corr_Pars := Continuation_Parameters.Create_End_Game;
                     70:     tol : constant double_float := 10.0**(-10);
                     71:     epslop : constant double_float := 10.0**(-6);
                     72:     wc : natural;
                     73:     sum,allsum : Standard_Complex_Vectors.Vector(sa(sa'first).sol.v'range);
                     74:
                     75:     procedure CCont2 is
                     76:       new Circular_Single_Conditioned_Reporting_Continue
                     77:             (Max_Norm,Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
                     78:
                     79:   begin
                     80:     for i in sa'range loop
                     81:       if AbsVal(Create(1.0) - sa(i).sol.t) > epslop
                     82:        then
                     83:          CCont2(file,sa(i),target,tol,epslop,wc,mwc,sum,allsum,false,pen,cen);
                     84:          sa(i).sol.m := wc;
                     85:          sa(i).sol.v := allsum;
                     86:          sa(i).sol.t := target;
                     87:          Write_Statistics_and_Condition
                     88:            (file,i,sa(i).nstep,sa(i).nfail,sa(i).niter,sa(i).nsyst,sa(i).rcond);
                     89:          put(file,sa(i).sol.all);
                     90:       end if;
                     91:     end loop;
                     92:     Clear(sols);
                     93:     sols := Shallow_Create(sa);
                     94:   end Winding_Numbers;
                     95:
                     96: begin
                     97:  -- READING MAXIMUM WINDING NUMBER :
                     98:   new_line;
                     99:   put("Give the maximum winding number : "); Read_Natural(max_wc);
                    100:  -- READING THE START SYSTEM :
                    101:   new_line;
                    102:   put_line("Reading the name of the file for start system.");
                    103:   Read_Name_and_Open_File(infile); get(infile,q);
                    104:   new_line;
                    105:   put_line(file,"THE START SYSTEM :");
                    106:   new_line(file); put(file,q); new_line(file);
                    107:  -- CONSTRUCTING THE HOMOTOPY AND TUNING OF PARAMETERS :
                    108:   Copy(p,pp);
                    109:   Driver_for_Homotopy_Construction(file,pp,q,sols,target);
                    110:   proj := (Number_of_Unknowns(q(q'first)) > q'last);
                    111:   Driver_for_Continuation_Parameters(file);
                    112:   new_line;
                    113:   put("Do you want intermediate output during continuation ? (y/n) ");
                    114:   Ask_Yes_or_No(ans);
                    115:   if ans = 'y'
                    116:    then Driver_for_Process_io(file,oc);
                    117:   end if;
                    118:  -- COMPUTATION OF WINDING NUMBERS :
                    119:   new_line;
                    120:   put_line("No more input desired.  Computing winding numbers ...");
                    121:   put_line("The program can now run in the background.");
                    122:   new_line;
                    123:   tstart(timer);
                    124:   Winding_Numbers(file,sols,proj,target,max_wc);
                    125:   tstop(timer);
                    126:   new_line(file);
                    127:   print_times(file,timer,"computation of winding numbers");
                    128:   new_line(file);
                    129: end Driver_for_Winding_Numbers;