![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to global_dynamic_triangulation.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Dynlift|
Return to global_dynamic_triangulation.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Dynlift |
1.1 maekawa 1: with Standard_Random_Numbers; use Standard_Random_Numbers;
2: with Standard_Integer_VecVecs; use Standard_Integer_VecVecs;
3: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
4: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
5: with Vertices;
6:
7: package body Global_Dynamic_Triangulation is
8:
9: function Ordered_Initial_Simplex
10: ( n : in natural; pts : in List ) return List is
11:
12: -- DESCRIPTION :
13: -- The given list pts has at least n+1 points.
14: -- The order of the points should be taken into account when
15: -- constructing the initial cell.
16: -- Returns the list of points which span the initial simplex.
17:
18: -- ALGORITHM :
19: -- Determines the rank of the matrix, defined by the first n+1 points.
20: -- If necessary, more extremal points will be sought.
21:
22: tmp,res,res_last : List;
23: mat : matrix(1..n,1..n);
24: sh,pt : Link_to_Vector;
25: rnk : natural;
26:
27: begin
28: sh := Head_Of(pts);
29: tmp := Tail_Of(pts);
30: for i in 1..n loop
31: pt := Head_Of(tmp);
32: for j in pt'range loop
33: mat(i,j) := pt(j) - sh(j);
34: end loop;
35: tmp := Tail_Of(tmp);
36: end loop;
37: rnk := Rank(mat);
38: if rnk = n
39: then tmp := pts;
40: for i in 1..n+1 loop
41: pt := Head_Of(tmp);
42: Append(res,res_last,pt.all);
43: tmp := Tail_Of(tmp);
44: end loop;
45: end if;
46: return res;
47: end Ordered_Initial_Simplex;
48:
49: procedure Initial_Simplex
50: ( pts : in List; order : in boolean;
51: s : out Simplex; rest : in out List ) is
52:
53: n : natural;
54: res,tmp,rest_last : List;
55: pt : Link_to_Vector;
56:
57: begin
58: if Is_Null(pts)
59: then
60: s := Null_Simplex;
61: else
62: n := Head_Of(pts)'last;
63: if Length_Of(pts) < n+1
64: then
65: s := Null_Simplex;
66: else
67: if order
68: then res := Ordered_Initial_Simplex(n,pts);
69: if Is_Null(res)
70: then res := Vertices.Extremal_Points(n,pts);
71: end if;
72: else res := Vertices.Extremal_Points(n,pts);
73: end if;
74: if Length_Of(res) = n+1
75: then
76: declare
77: points : VecVec(1..n+1);
78: begin
79: tmp := res;
80: for k in points'range loop
81: pt := Head_Of(tmp);
82: points(k) := new vector(1..n+1);
83: points(k)(1..n) := pt.all;
84: points(k)(n+1) := 0;
85: tmp := Tail_Of(tmp);
86: end loop;
87: s := Create(points);
88: end;
89: -- construct the list with the rest of the points :
90: tmp := pts;
91: while not Is_Null(tmp) loop
92: pt := Head_Of(tmp);
93: if not Is_In(res,pt.all)
94: then Append(rest,rest_last,pt.all);
95: end if;
96: tmp := Tail_Of(tmp);
97: end loop;
98: else s := Null_Simplex;
99: end if;
100: Deep_Clear(res);
101: end if;
102: end if;
103: end Initial_Simplex;
104:
105: -- COMPUTING AN EXTREMAL POINT :
106:
107: function Max_Extreme ( l : List; k : natural ) return Link_to_Vector is
108:
109: res : Link_to_Vector := Head_Of(l);
110: mx : integer := res(k);
111: tmp : List := Tail_Of(l);
112: pt : Link_to_Vector;
113:
114: begin
115: while not Is_Null(tmp) loop
116: pt := Head_Of(tmp);
117: if pt(k) > mx
118: then mx := pt(k);
119: res := pt;
120: end if;
121: tmp := Tail_Of(tmp);
122: end loop;
123: return res;
124: end Max_Extreme;
125:
126: function Max_Extreme ( l : List; weights : vector ) return Link_to_Vector is
127:
128: res : Link_to_Vector := Head_Of(l);
129: tmp : List := Tail_Of(l);
130: mxs : integer := res.all*weights;
131: sp : integer;
132: pt : Link_to_Vector;
133:
134: begin
135: while not Is_Null(tmp) loop
136: pt := Head_Of(tmp);
137: sp := pt.all*weights;
138: if sp > mxs
139: then mxs := sp;
140: res := pt;
141: end if;
142: tmp := Tail_Of(tmp);
143: end loop;
144: return res;
145: end Max_Extreme;
146:
147: function Max_Extreme ( l : List; n : natural; low,upp : integer )
148: return Link_to_Vector is
149:
150: w : Vector(1..n);
151:
152: begin
153: for k in w'range loop
154: w(k) := Random(low,upp);
155: end loop;
156: return Max_Extreme(l,w);
157: end Max_Extreme;
158:
159: end Global_Dynamic_Triangulation;