Annotation of OpenXM_contrib/PHC/Ada/Continuation/drivers_for_path_directions.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 Standard_Floating_Numbers; use Standard_Floating_Numbers;
5: with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
6: with Numbers_io; use Numbers_io;
7: with Standard_Complex_Vectors;
8: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
9: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
10: with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
11: with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
12: with Homotopy;
13: with Projective_Transformations; use Projective_Transformations;
14: with Continuation_Data; use Continuation_Data;
15: with Continuation_Parameters; use Continuation_Parameters;
16: with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
17:
18: package body Drivers_for_Path_Directions is
19:
20: -- AUXILIARIES TO INSTANTIATE :
21:
22: function Scale ( v : Vector ) return Vector is
23:
24: -- DESCRIPTION :
25: -- Divides the vector by its largest component.
26:
27: res : Vector(v'range);
28: tol : constant double_float := 10.0**(-12);
29: ind : integer := v'first;
30: max : double_float := abs(v(ind));
31:
32: begin
33: for i in v'range loop
34: if abs(v(i)) > max
35: then max := abs(v(i));
36: ind := i;
37: end if;
38: end loop;
39: if max > tol
40: then for i in v'range loop
41: res(i) := v(i)/max;
42: end loop;
43: end if;
44: return res;
45: end Scale;
46:
47: function Toric_Evaluator ( x : Standard_Complex_Vectors.Vector;
48: t : Complex_Number )
49: return Standard_Complex_Vectors.Vector is
50: begin
51: return Homotopy.Eval(x,t);
52: end Toric_Evaluator;
53:
54: function Toric_Differentiator
55: ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
56: return Standard_Complex_Vectors.Vector is
57: begin
58: return Homotopy.Eval(x,t);
59: end Toric_Differentiator;
60:
61: function Toric_Differentiator
62: ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
63: return Standard_Complex_Matrices.Matrix is
64: begin
65: return Homotopy.Diff(x,t);
66: end Toric_Differentiator;
67:
68: -- TARGET ROUTINES :
69:
70: procedure Init_Path_Directions
71: ( n,nv : in natural; v : in out Link_to_VecVec;
72: errv : in out Link_to_Vector ) is
73:
74: begin
75: v := new VecVec(1..nv);
76: for i in v'range loop
77: v(i) := new Vector'(1..n => 0.0);
78: end loop;
79: errv := new Vector'(1..nv => 1.0);
80: end Init_Path_Directions;
81:
82: procedure Toric_Continue
83: ( file : in file_type; sols : in out Solution_List;
84: proj,report : in boolean; v : in out VecVec;
85: errv : in out Vector; target : in Complex_Number ) is
86:
87: -- DESCRIPTION :
88: -- Performs the continuation with online toric compactifications.
89:
90: timer : timing_widget;
91: h : constant Poly_Sys := Homotopy.Eval(target); --Homotopy.Homotopy_System;
92: n : constant natural := h'length;
93: hh : Poly_Sys(h'range);
94:
95: procedure Sil_Cont is
96: new Silent_Toric_Continue(Max_Norm,Toric_Evaluator,
97: Toric_Differentiator,Toric_Differentiator);
98: procedure Rep_Cont is
99: new Reporting_Toric_Continue(Max_Norm,Toric_Evaluator,
100: Toric_Differentiator,Toric_Differentiator);
101: begin
102: tstart(timer);
103: if report
104: then Rep_Cont(file,sols,proj,v,errv,target);
105: else Sil_Cont(sols,proj,v,errv,target);
106: end if;
107: tstop(timer);
108: new_line(file);
109: print_times(file,timer,"toric continuation"); new_line(file);
110: end Toric_Continue;
111:
112: procedure Write_Directions
113: ( file : in file_type; v : in VecVec; errv : in Vector ) is
114:
115: procedure Write ( file : in file_type; v : in Vector ) is
116: begin
117: for i in v'range loop
118: put(file,v(i)); new_line(file);
119: end loop;
120: end Write;
121:
122: procedure Write_Direction
123: ( file : in file_type;
124: v : in Vector; error : in double_float; i : natural ) is
125:
126: begin
127: put(file,"Computed direction of path ");
128: put(file,i,1); put_line(file," :"); Write(file,v);
129: put(file,"with error : "); put(file,error); new_line(file);
130: end Write_Direction;
131:
132: begin
133: for i in v'range loop
134: Write_Direction(file,v(i).all,errv(i),i);
135: end loop;
136: end Write_Directions;
137:
138: end Drivers_for_Path_Directions;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>