Annotation of OpenXM_contrib/PHC/Ada/Continuation/black_polynomial_continuations.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_Complex_Numbers; use Standard_Complex_Numbers;
7: with Standard_Complex_Numbers_io; use Standard_Complex_Numbers_io;
8: with Standard_Random_Numbers; use Standard_Random_Numbers;
9: with 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 Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
14: with Standard_Complex_Solutions_io; use Standard_Complex_Solutions_io;
15: with Scaling; use Scaling;
16: with Projective_Transformations; use Projective_Transformations;
17: with Continuation_Parameters;
18: with Continuation_Parameters_io;
19: with Homotopy,Process_io; use Process_io;
20: with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
21: with Standard_Root_Refiners; use Standard_Root_Refiners;
22: with Scanners_for_Continuation; use Scanners_for_Continuation;
23:
24: package body Black_Polynomial_Continuations is
25:
26: procedure Scan_Input
27: ( infile,outfile : in file_type; p,q : in out Link_to_Poly_Sys;
28: sols : in out Solution_List; arti : out boolean ) is
29:
30: -- DESCRIPTION :
31: -- Scans the input file for target system and, if the homotopy is
32: -- artificial (in that case arti = true, otherwise arti = false),
33: -- for a start system. In both cases, start solutions are required.
34:
35: found,artificial : boolean;
36:
37: begin
38: get(infile,p);
39: put(outfile,p'last,p.all);
40: artificial := (Number_of_Unknowns(p(p'first)) = p'last);
41: if artificial
42: then Scan_and_Skip(infile,"START SYSTEM",found);
43: if found
44: then get(infile,q);
45: new_line(outfile);
46: put_line(outfile,"THE START SYSTEM : ");
47: new_line(outfile);
48: put_line(outfile,q.all);
49: end if;
50: end if;
51: Scan_and_Skip(infile,"SOLUTIONS",found);
52: if found
53: then get(infile,sols);
54: new_line(outfile);
55: put_line(outfile,"THE START SOLUTIONS : ");
56: new_line(outfile);
57: put(outfile,Length_Of(sols),Head_Of(sols).n,sols);
58: new_line(outfile);
59: end if;
60: arti := artificial;
61: end Scan_Input;
62:
63: procedure Set_Homotopy_Parameters
64: ( file : in file_type; k : in out natural;
65: a,t : in out Complex_Number; prt : in out boolean ) is
66:
67: -- DESCRIPTION :
68: -- Sets the default values for the homotopy parameters.
69:
70: begin
71: k := 2;
72: a := Random1;
73: t := Create(1.0);
74: prt := false;
75: new_line(file);
76: put_line(file,"HOMOTOPY PARAMETERS :");
77: put(file," k : "); put(file,k,2); new_line(file);
78: put(file," a : "); put(file,a); new_line(file);
79: put(file," t : "); put(file,t); new_line(file);
80: if prt
81: then put_line(file," projective transformation");
82: else put_line(file," no projective transformation");
83: end if;
84: end Set_Homotopy_Parameters;
85:
86: procedure Tune_Continuation_Parameters ( outfile : in file_type ) is
87:
88: -- DESCRIPTION :
89: -- Scans the input file for continuation parameters and the
90: -- output parameter.
91:
92: begin
93: Continuation_Parameters.Tune(2);
94: new_line(outfile);
95: put_line(outfile,"****************** CURRENT CONTINUATION PARAMETERS "
96: & "*****************");
97: Continuation_Parameters_io.put(outfile);
98: put_line(outfile,"***************************************************"
99: & "*****************");
100: Process_io.Set_Output_Code(nil);
101: end Tune_Continuation_Parameters;
102:
103: procedure Black_Box_Refine
104: ( outfile : in file_type; p : in Poly_Sys;
105: artificial : in boolean; target : in Complex_Number;
106: k : in natural; sols,refsols : in out Solution_List ) is
107:
108: -- DESCRIPTION :
109: -- Refines the roots in sols w.r.t. the system p.
110:
111: epsxa,epsfa : constant double_float := 10.0**(-8);
112: tolsing : constant double_float := 10.0**(-8);
113: len,nb : natural := 0;
114:
115: begin
116: if Length_Of(sols) > 0
117: then
118: if artificial
119: then if not Is_Null(sols) and then Head_Of(sols).n > p'last
120: then Affine_Transformation(sols);
121: end if;
122: Reporting_Root_Refiner
123: (outfile,p,sols,refsols,epsxa,epsfa,tolsing,nb,5,false);
124: else declare
125: pt : Poly_Sys(p'range);
126: begin
127: pt := Eval(p,target,k);
128: Reporting_Root_Refiner
129: (outfile,pt,sols,refsols,epsxa,epsfa,tolsing,nb,5,false);
130: Clear(pt);
131: end;
132: end if;
133: end if;
134: end Black_Box_Refine;
135:
136: procedure Black_Box_Polynomial_Continuation
137: ( infile,outfile : in file_type; pocotime : out duration ) is
138:
139: p,q,sp : Link_to_Poly_Sys;
140: sols,refsols : Solution_List;
141: timer : timing_widget;
142: k : natural := 0;
143: a,target : Complex_Number;
144: proj,artificial : boolean;
145: rcond : double_float;
146: scalecoeff : Standard_Complex_Vectors.Link_to_Vector;
147:
148: procedure Cont is
149: new Reporting_Continue(Max_Norm,
150: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
151:
152: begin
153: Scan_Input(infile,outfile,p,q,sols,artificial);
154: scalecoeff := new Standard_Complex_Vectors.Vector(1..2*p'length);
155: sp := new Poly_Sys(p'range);
156: Copy(p.all,sp.all);
157: Scale(sp.all,2,false,rcond,scalecoeff.all);
158: Set_Homotopy_Parameters(outfile,k,a,target,proj);
159: if artificial
160: then Homotopy.Create(sp.all,q.all,k,a);
161: else Homotopy.Create(sp.all,k);
162: target := a;
163: end if;
164: Tune_Continuation_Parameters(outfile);
165: new_line(outfile);
166: put_line(outfile,"THE SCALED SOLUTIONS :");
167: new_line(outfile);
168: tstart(timer);
169: Cont(outfile,sols,proj,target);
170: tstop(timer);
171: new_line(outfile);
172: print_times(outfile,timer,"continuation");
173: pocotime := Elapsed_User_Time(timer);
174: Scale(2,scalecoeff.all,sols);
175: Clear(sp);
176: Black_Box_Refine(outfile,p.all,artificial,target,k,sols,refsols);
177: end Black_Box_Polynomial_Continuation;
178:
179: procedure Black_Box_Polynomial_Continuation
180: ( outfile : in file_type;
181: p,q : in Poly_Sys; sols : in out Solution_List;
182: pocotime : out duration ) is
183:
184: refsols : Solution_List;
185: timer : timing_widget;
186: k : natural := 0;
187: a,target : Complex_Number := Create(0.0);
188: proj : boolean;
189: rcond : double_float;
190: scalecoeff : Standard_Complex_Vectors.Vector(1..2*p'length);
191: sp : Poly_Sys(p'range);
192:
193: procedure Cont is
194: new Reporting_Continue(Max_Norm,
195: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
196:
197: begin
198: Set_Homotopy_Parameters(outfile,k,a,target,proj);
199: Copy(p,sp);
200: Scale(sp,2,false,rcond,scalecoeff);
201: Homotopy.Create(sp,q,k,a);
202: Tune_Continuation_Parameters(outfile);
203: new_line(outfile);
204: put_line(outfile,"THE SCALED SOLUTIONS :");
205: new_line(outfile);
206: tstart(timer);
207: Cont(outfile,sols,proj,target);
208: tstop(timer);
209: new_line(outfile);
210: print_times(outfile,timer,"continuation");
211: pocotime := Elapsed_User_Time(timer);
212: Scale(2,scalecoeff,sols);
213: Clear(sp);
214: Black_Box_Refine(outfile,p,true,target,k,sols,refsols);
215: sols := refsols;
216: end Black_Box_Polynomial_Continuation;
217:
218: end Black_Polynomial_Continuations;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to black_polynomial_continuations.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Continuation