Annotation of OpenXM_contrib/PHC/Ada/Schubert/specialization_of_planes.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
2: with Standard_Random_Numbers; use Standard_Random_Numbers;
3: with Standard_Natural_Vectors;
4: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
5:
6: package body Specialization_of_Planes is
7:
8: function Random_Upper_Triangular
9: ( n : natural ) return Standard_Complex_Matrices.Matrix is
10:
11: res : Standard_Complex_Matrices.Matrix(1..n,1..n);
12:
13: begin
14: for j in 1..n loop -- assign values to jth column
15: for i in 1..n-j loop
16: res(i,j) := Random1; -- randoms above anti-diagonal
17: end loop;
18: res(n-j+1,j) := Create(1.0); -- 1 = anti-diagonal element
19: for i in n-j+2..n loop
20: res(i,j) := Create(0.0); -- zeros under anti-diagonal
21: end loop;
22: end loop;
23: return res;
24: end Random_Upper_Triangular;
25:
26: function Random_Lower_Triangular
27: ( n : natural ) return Standard_Complex_Matrices.Matrix is
28:
29: res : Standard_Complex_Matrices.Matrix(1..n,1..n);
30:
31: begin
32: for j in 1..n loop -- assign values to jth column
33: for i in 1..(j-1) loop
34: res(i,j) := Create(0.0); -- zeros above diagonal
35: end loop;
36: res(j,j) := Create(1.0); -- 1 = diagonal element
37: for i in (j+1)..n loop
38: res(i,j) := Random1; -- randoms under diagonal
39: end loop;
40: end loop;
41: return res;
42: end Random_Lower_Triangular;
43:
44: function U_Matrix ( F : Standard_Complex_Matrices.Matrix; b : Bracket )
45: return Standard_Complex_Matrices.Matrix is
46:
47: m : constant natural := F'length(1) - b'length;
48: res : Standard_Complex_Matrices.Matrix(F'range(1),1..m);
49: rvf : Standard_Complex_Matrices.Matrix(F'range(1),F'range(2)) := F;
50: rng : constant natural := F'length(2) - b(b'last);
51: tmp : Complex_Number;
52: ind : natural := 1;
53: cnt : natural := 0;
54:
55: begin
56: for j in 1..(rng/2) loop -- reverse last columns
57: for i in F'range(1) loop
58: tmp := rvf(i,F'last(2)-j+1);
59: rvf(i,F'last(2)-j+1) := rvf(i,F'last(2)-rng+j);
60: rvf(i,F'last(2)-rng+j) := tmp;
61: end loop;
62: end loop;
63: for j in F'range(2) loop -- remove columns indexed by b
64: if ((ind <= b'last) and then (j = b(ind)))
65: then ind := ind+1;
66: else cnt := cnt+1;
67: for i in F'range(1) loop
68: res(i,cnt) := rvf(i,j);
69: end loop;
70: end if;
71: end loop;
72: return res;
73: end U_Matrix;
74:
75: function Special_Plane ( m : natural; b : Bracket )
76: return Standard_Complex_Matrices.Matrix is
77:
78: p : constant natural := b'length;
79: n : constant natural := m+p;
80: res : Standard_Complex_Matrices.Matrix(1..n,1..m);
81: row,col : natural;
82:
83: begin
84: row := 1; col := 0;
85: for i in res'range(1) loop
86: for j in res'range(2) loop
87: res(i,j) := Create(0.0);
88: end loop;
89: if ((row <= p) and then (b(row) = i))
90: then row := row + 1;
91: else col := col + 1;
92: res(i,col) := Create(1.0);
93: end if;
94: end loop;
95: return res;
96: end Special_Plane;
97:
98: function Special_Bottom_Plane ( m : natural; b : Bracket )
99: return Standard_Complex_Matrices.Matrix is
100:
101: p : constant natural := b'length;
102: n : constant natural := m+p;
103: res : Standard_Complex_Matrices.Matrix(1..n,1..m);
104: row,col : natural;
105:
106: begin
107: row := 1; col := 0;
108: for i in res'range(1) loop
109: if ((row <= p) and then (b(row) = i))
110: then row := row + 1;
111: else col := col + 1;
112: for k in 1..i-1 loop -- randoms above the diagonal
113: res(k,col) := Random1;
114: end loop;
115: res(i,col) := Create(1.0);
116: for k in i+1..n loop -- zeros below the diagonal
117: res(k,col) := Create(0.0);
118: end loop;
119: end if;
120: end loop;
121: return res;
122: end Special_Bottom_Plane;
123:
124: function Special_Top_Plane ( m : natural; b : Bracket )
125: return Standard_Complex_Matrices.Matrix is
126:
127: p : constant natural := b'length;
128: n : constant natural := m+p;
129: res : Standard_Complex_Matrices.Matrix(1..n,1..m);
130: row,col : natural;
131:
132: begin
133: row := 1; col := 0;
134: for i in res'range(1) loop
135: if ((row <= p) and then (b(row) = i))
136: then row := row + 1;
137: else col := col + 1;
138: for k in 1..i-1 loop -- zeros above the diagonal
139: res(k,col) := Create(0.0);
140: end loop;
141: res(i,col) := Create(1.0);
142: for k in i+1..n loop -- randoms below the diagonal
143: res(k,col) := Random1;
144: end loop;
145: end if;
146: end loop;
147: return res;
148: end Special_Top_Plane;
149:
150: function Special_Plane
151: ( n,m,k : natural; b : Bracket;
152: special : in Standard_Complex_Matrices.Matrix )
153: return Standard_Complex_Matrices.Matrix is
154:
155: res : Standard_Complex_Matrices.Matrix(1..n,1..m+1-k);
156: ran : Complex_Number;
157:
158: begin
159: for i in res'range(1) loop -- initialize
160: for j in res'range(2) loop
161: res(i,j) := Create(0.0);
162: end loop;
163: end loop;
164: for j in res'range(2) loop -- build j-th column
165: for k in b'range loop
166: ran := Random1; -- random for column k of mat
167: for i in special'range(1) loop
168: res(i,j) := res(i,j) + ran*special(i,b(k));
169: end loop;
170: end loop;
171: end loop;
172: return res;
173: end Special_Plane;
174:
175: function Special_Bottom_Plane ( n,m,k : natural; b : Bracket )
176: return Standard_Complex_Matrices.Matrix is
177:
178: mat : Standard_Complex_Matrices.Matrix(1..n,b'range);
179:
180: begin
181: for j in b'range loop -- j-th column of matrix
182: for i in 1..b(j) loop
183: mat(i,j) := Random1; -- random numbers above row b(j)
184: end loop;
185: for i in b(j)+1..n loop
186: mat(i,j) := Create(0.0); -- zeros below row b(j)
187: end loop;
188: end loop;
189: return Special_Plane(n,m,k,b,mat);
190: end Special_Bottom_Plane;
191:
192: function Special_Top_Plane ( n,m,k : natural; b : Bracket )
193: return Standard_Complex_Matrices.Matrix is
194:
195: mat : Standard_Complex_Matrices.Matrix(1..n,b'range);
196:
197: begin
198: for j in b'range loop -- j-th column of matrix
199: for i in 1..b(j)-1 loop
200: mat(i,j) := Create(0.0); -- zeros below row b(j)
201: end loop;
202: for i in b(j)..n loop
203: mat(i,j) := Random1; -- random numbers below row b(j)
204: end loop;
205: end loop;
206: return Special_Plane(n,m,k,b,mat);
207: end Special_Top_Plane;
208:
209: function Homotopy ( dim : natural; start,target : Complex_Number )
210: return Poly is
211:
212: -- DESCRIPTION :
213: -- Returns the polynomial start*(1-t) + target*t, where t is the
214: -- is the last variable of index dim.
215: -- This procedure is an auxiliary to building the moving U-matrices.
216:
217: res : Poly;
218: t : Term;
219: tdg : Degrees := new Standard_Natural_Vectors.Vector'(1..dim => 0);
220:
221: begin
222: t.cf := start;
223: t.dg := tdg;
224: res := Create(t); -- res = start
225: tdg(tdg'last) := 1; -- introduce t
226: t.dg := tdg;
227: Sub(res,t); -- res = (1-t)*start
228: t.cf := target;
229: Add(res,t); -- res = (1-t)*start + t*target
230: Clear(tdg);
231: return res;
232: end Homotopy;
233:
234: function Constant_Poly ( dim : natural; c : Complex_Number ) return Poly is
235:
236: -- DESCRIPTION :
237: -- Returns the constant c represented as polynomial with as many
238: -- variables as the given dimension.
239:
240: res : Poly;
241: t : Term;
242: tdg : Degrees := new Standard_Natural_Vectors.Vector'(1..dim => 0);
243:
244: begin
245: t.cf := c;
246: t.dg := tdg;
247: res := Create(t);
248: return res;
249: end Constant_Poly;
250:
251: function Moving_U_Matrix
252: ( n : natural; U,L : Standard_Complex_Matrices.Matrix )
253: return Standard_Complex_Poly_Matrices.Matrix is
254:
255: res : Standard_Complex_Poly_Matrices.Matrix(L'range(1),L'range(2));
256: t : Term;
257:
258: begin
259: for i in res'range(1) loop
260: for j in res'range(2) loop
261: res(i,j) := Homotopy(n,U(i,j),L(i,j));
262: end loop;
263: end loop;
264: return res;
265: end Moving_U_Matrix;
266:
267: function Moving_U_Matrix
268: ( U : Standard_Complex_Matrices.Matrix;
269: i,r : natural; b : bracket )
270: return Standard_Complex_Poly_Matrices.Matrix is
271:
272: p : constant natural := b'last;
273: m : constant natural := U'length(1) - p;
274: dim : constant natural := (m+p)*p+1;
275: res : Standard_Complex_Poly_Matrices.Matrix(U'range(1),1..m+1-r);
276:
277: begin
278: for j in res'range(2) loop
279: if j+i-1 < b(b'last) - b'last
280: then for k in res'range(1) loop
281: res(k,j) := Homotopy(dim,U(k,j+i),U(k,j+i-1));
282: end loop;
283: elsif j+i-1 = b(b'last) - b'last
284: then for k in res'range(1) loop
285: res(k,j) := Homotopy(dim,U(k,m+1+i-r),U(k,j+i-1));
286: end loop;
287: else for k in res'range(1) loop
288: res(k,j) := Constant_Poly(dim,U(k,j+i-1));
289: end loop;
290: end if;
291: end loop;
292: return res;
293: end Moving_U_Matrix;
294:
295: function Slice
296: ( M : Standard_Complex_Poly_Matrices.Matrix;
297: ind : natural ) return Standard_Complex_Poly_Matrices.Matrix is
298:
299: -- DESCRIPTION :
300: -- Returns the columns of M up to the given index.
301:
302: res : Standard_Complex_Poly_Matrices.Matrix(M'range(1),M'first(2)..ind);
303:
304: begin
305: for j in res'range(2) loop
306: for i in res'range(1) loop
307: Copy(M(i,j),res(i,j));
308: end loop;
309: end loop;
310: return res;
311: end Slice;
312:
313: function Lower_Section
314: ( M : Standard_Complex_Poly_Matrices.Matrix;
315: row : natural ) return Standard_Complex_Poly_Matrices.Matrix is
316:
317: res : Standard_Complex_Poly_Matrices.Matrix(M'range(1),M'range(2));
318: cnt : natural := M'first(2)-1;
319: add : boolean;
320:
321: begin
322: for j in M'range(2) loop
323: for i in (row+1)..M'last(1) loop
324: add := (M(i,j) = Null_Poly);
325: exit when not add;
326: end loop;
327: if add
328: then cnt := cnt+1;
329: for i in M'range(1) loop
330: res(i,cnt) := M(i,j);
331: end loop;
332: end if;
333: end loop;
334: return Slice(res,cnt);
335: end Lower_Section;
336:
337: function Upper_Section
338: ( M : Standard_Complex_Poly_Matrices.Matrix;
339: row : natural ) return Standard_Complex_Poly_Matrices.Matrix is
340:
341: res : Standard_Complex_Poly_Matrices.Matrix(M'range(1),M'range(2));
342: cnt : natural := M'first(2)-1;
343: add : boolean;
344:
345: begin
346: for j in M'range(2) loop
347: for i in M'first(1)..(row-1) loop
348: add := (M(i,j) = Null_Poly);
349: exit when not add;
350: end loop;
351: if add
352: then cnt := cnt+1;
353: for i in M'range(1) loop
354: res(i,cnt) := M(i,j);
355: end loop;
356: end if;
357: end loop;
358: return Slice(res,cnt);
359: end Upper_Section;
360:
361: end Specialization_of_Planes;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to specialization_of_planes.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Schubert