Annotation of OpenXM_contrib/PHC/Ada/Schubert/pieri_deformations.adb, Revision 1.1.1.1
1.1 maekawa 1: with integer_io; use integer_io;
2: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
3: with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
4: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
5: with Standard_Natural_Matrices;
6: with Standard_Natural_Matrices_io; use Standard_Natural_Matrices_io;
7: with Standard_Complex_Vectors;
8: with Standard_Complex_Vectors_io; use Standard_Complex_Vectors_io;
9: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
10: with Standard_Complex_Matrices_io; use Standard_Complex_Matrices_io;
11: with Standard_Complex_Poly_Matrices;
12: with Standard_Complex_Poly_Matrices_io; use Standard_Complex_Poly_Matrices_io;
13: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
14: with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
15: with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
16: with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
17: with Brackets,Brackets_io; use Brackets,Brackets_io;
18: with Bracket_Monomials; use Bracket_Monomials;
19: with Bracket_Monomials_io; use Bracket_Monomials_io;
20: with Bracket_Polynomials; use Bracket_Polynomials;
21: with Bracket_Polynomials_io; use Bracket_Polynomials_io;
22: with Bracket_Systems; use Bracket_Systems;
23: with Pieri_Trees,Pieri_Trees_io; use Pieri_Trees,Pieri_Trees_io;
24: with Symbolic_Minor_Equations; use Symbolic_Minor_Equations;
25: with Numeric_Minor_Equations; use Numeric_Minor_Equations;
26: with Determinantal_Systems; use Determinantal_Systems;
27: with Solve_Pieri_Leaves;
28: with Plane_Representations; use Plane_Representations;
29: with Specialization_of_Planes; use Specialization_of_Planes;
30: with Pieri_Continuation; use Pieri_Continuation;
31:
32: package body Pieri_Deformations is
33:
34: function Coordinate_Bracket
35: ( nd : Pieri_Node; jmp : natural ) return Bracket is
36:
37: -- DESCRIPTION :
38: -- On return we find the bracket determines the local coordinates.
39: -- This bracket depends whether jmp = nd.node'last or not.
40:
41: begin
42: if Is_Leaf(nd) or jmp < nd.node'last
43: then return nd.node;
44: else return Upper_Jump_Decrease(nd);
45: end if;
46: end Coordinate_Bracket;
47:
48: procedure Test_Solution
49: ( file : in file_type; nd1,nd2 : in Pieri_Node;
50: id : in natural; l1,l2 : in VecMat; l : in Matrix;
51: x : Standard_Complex_Poly_Matrices.Matrix;
52: solpla : in Matrix ) is
53:
54: -- DESCRIPTION :
55: -- Evaluates the system that represents the intersection conditions
56: -- for the planes in l1,l2 and the plane l in the given solution plane.
57:
58: sys : Link_to_Poly_Sys := new Poly_Sys'(Polynomial_Equations(l,x));
59: eva : Standard_Complex_Vectors.Vector(sys'range);
60: fst : natural;
61:
62: begin
63: if Is_Leaf(nd1)
64: then fst := l1'last+1;
65: elsif nd1.i = 0
66: then fst := nd1.c;
67: else fst := nd1.c+1;
68: end if;
69: for i in fst..l1'last loop
70: Concat(sys,Polynomial_Equations(l1(i).all,x));
71: end loop;
72: if Is_Leaf(nd2)
73: then fst := l2'last+1;
74: elsif nd2.i = 0
75: then fst := nd2.c;
76: else fst := nd2.c+1;
77: end if;
78: for i in fst..l2'last loop
79: Concat(sys,Polynomial_Equations(l2(i).all,x));
80: end loop;
81: -- put_line(file,"The polynomial system that has been solved :");
82: -- put_line(file,sys.all);
83: eva := Eval(sys.all,Vector_Rep(solpla));
84: put(file,"Residual at the target system :");
85: put(file,Max_Norm(eva),3);
86: Clear(sys);
87: if (((nd1.i = 0) and (nd1.c = 1))
88: or else ((nd2.i = 0) and (nd2.c = 1)))
89: then put(file," PAIR "); put(file,id,1); put_line(file," DONE.");
90: else put_line(file,".");
91: end if;
92: end Test_Solution;
93:
94: procedure Test_Solutions
95: ( file : in file_type;
96: l1,l2 : in VecMat; l : in Matrix;
97: x : Standard_Complex_Poly_Matrices.Matrix;
98: solpla : in VecMat ) is
99:
100: -- DESCRIPTION :
101: -- Evaluates the system that represents the intersection conditions
102: -- for the planes in l1,l2 and the plane l in the given solution plane.
103:
104: sys : Link_to_Poly_Sys := new Poly_Sys'(Polynomial_Equations(l,x));
105:
106: begin
107: for i in l1'range loop
108: Concat(sys,Polynomial_Equations(l1(i).all,x));
109: end loop;
110: for i in l2'range loop
111: Concat(sys,Polynomial_Equations(l2(i).all,x));
112: end loop;
113: put_line(file,"The polynomial system that has been solved :");
114: put_line(file,sys.all);
115: for i in solpla'range loop
116: exit when (solpla(i) = null);
117: put(file,"Solution plane no. "); put(file,i,1); put_line(file," :");
118: put(file,solpla(i).all,2);
119: put_line(file,"The system evaluated at the plane :");
120: put_line(file,Eval(sys.all,Vector_Rep(solpla(i).all)));
121: end loop;
122: Clear(sys);
123: end Test_Solutions;
124:
125: procedure Start_at_Leaves
126: ( file : in file_type; pnd : in Paired_Nodes;
127: ln : in Matrix; sol : in out Matrix ) is
128:
129: -- DESCRIPTION :
130: -- Computes the start solutions at the leaves.
131:
132: -- xpm : Standard_Complex_Poly_Matrices.Matrix(sol'range(1),sol'range(2))
133: -- := Schubert_Pattern(sol'last(1),pnd.left.node,pnd.right.node);
134: -- sys : Link_to_Poly_Sys := new Poly_Sys'(Polynomial_Equations(ln,xpm));
135:
136: begin
137: sol := Solve_Pieri_Leaves(file,pnd.left.node,pnd.right.node,ln);
138: put_line(file,"The solution plane at the leaves : ");
139: put(file,sol,2);
140: -- put_line(file,"The polynomial equations : "); put(file,sys.all);
141: -- put(file,"The solution evaluated at the polynomial equations : ");
142: -- put(file,Evaluate(sys.all,sol),3); new_line(file);
143: -- Standard_Complex_Poly_Matrices.Clear(xpm);
144: -- Clear(sys);
145: end Start_at_Leaves;
146:
147: function Moving_U_Matrix
148: ( file : in file_type; outlog : in boolean;
149: nd : Pieri_Node; lb : Bracket;
150: f : Standard_Complex_Matrices.Matrix; l : VecMat )
151: return Standard_Complex_Poly_Matrices.Matrix is
152:
153: -- DESCRIPTION :
154: -- Returns the moving U-matrix.
155:
156: -- REQUIRED : nd.c > 0.
157:
158: n : constant natural := f'length(1);
159: p : constant natural := lb'length;
160: nva : constant natural := n*p+1;
161: m : constant natural := n-p;
162: lc : constant natural := l(nd.c)'length(2);
163: r : natural;
164: U : constant Standard_Complex_Matrices.Matrix := U_Matrix(f,lb);
165:
166: begin
167: if outlog
168: then put_line(file,"The U-matrix : "); put(file,U,2);
169: end if;
170: if nd.i = 0
171: then return Moving_U_Matrix(nva,U,l(nd.c).all);
172: else r := m+1 - lc;
173: return Moving_U_Matrix(U,nd.i,r,lb);
174: end if;
175: end Moving_U_Matrix;
176:
177: function Moving_Cycle ( movU,x : Standard_Complex_Poly_Matrices.Matrix )
178: return Poly_Sys is
179:
180: -- DESCRIPTION :
181: -- Returns the moving cycle expaned as polynomial system.
182:
183: n : constant natural := x'length(1);
184: p : constant natural := x'length(2);
185: kd : constant natural := p + movU'length(2);
186: bm : Bracket_Monomial := Maximal_Minors(n,kd);
187: bs : Bracket_System(0..Number_of_Brackets(bm))
188: := Minor_Equations(kd,kd-p,bm);
189: res : constant Poly_Sys := Lifted_Expanded_Minors(movU,x,bs);
190:
191: begin
192: Clear(bm); Clear(bs);
193: return res;
194: end Moving_Cycle;
195:
196: function One_Moving_Cycle
197: ( file : in file_type; nd : Pieri_Node; lb : Bracket;
198: f : Standard_Complex_Matrices.Matrix; l : VecMat;
199: x : Standard_Complex_Poly_Matrices.Matrix; outlog : boolean )
200: return Poly_Sys is
201:
202: -- DESCRIPTION :
203: -- Returns the moving part for one node in case i = 0.
204:
205: -- REQUIRED : nd.c > 0.
206:
207: lc : constant natural := l(nd.c)'length(2);
208: movU : Standard_Complex_Poly_Matrices.Matrix(x'range(1),1..lc)
209: := Moving_U_Matrix(file,outlog,nd,lb,f,l);
210: res : constant Poly_Sys := Moving_Cycle(movU,x);
211:
212: begin
213: if outlog
214: then put_line(file,"The moving U matrix : "); put(file,MovU);
215: put_line(file,"The equations of the moving cycle : "); put(file,res);
216: end if;
217: Standard_Complex_Poly_Matrices.Clear(movU);
218: return res;
219: end One_Moving_Cycle;
220:
221: function Left_Moving_Cycle
222: ( file : file_type; nd : Pieri_Node; lb : Bracket; jmp : natural;
223: f : Standard_Complex_Matrices.Matrix; l : VecMat;
224: x : Standard_Complex_Poly_Matrices.Matrix; outlog : boolean )
225: return Poly_Sys is
226:
227: -- DESCRIPTION :
228: -- Returns the moving part for the case i > 0 and jmp < p,
229: -- for the node from the first Pieri tree.
230:
231: -- REQUIRED : nd.c > 0.
232:
233: n : constant natural := x'length(1);
234: p : constant natural := x'length(2);
235: lc : constant natural := l(nd.c)'length(2);
236: movU : Standard_Complex_Poly_Matrices.Matrix(x'range(1),1..lc)
237: := Moving_U_Matrix(file,outlog,nd,lb,f,l);
238: movUsec : constant Standard_Complex_Poly_Matrices.Matrix
239: := Lower_Section(movU,n+1-lb(jmp));
240: res : constant Poly_Sys := Moving_Cycle(movUsec,x);
241:
242: begin
243: if outlog
244: then put_line(file,"The left moving U matrix : "); put(file,movU);
245: put_line(file,"The left moving cycle : "); put(file,movUsec);
246: put_line(file,"The equations of the moving cycle : "); put(file,res);
247: end if;
248: Standard_Complex_Poly_Matrices.Clear(movU);
249: return res;
250: end Left_Moving_Cycle;
251:
252: function Right_Moving_Cycle
253: ( file : file_type; nd : Pieri_Node; lb : Bracket; jmp : natural;
254: f : Standard_Complex_Matrices.Matrix; l : VecMat;
255: x : Standard_Complex_Poly_Matrices.Matrix; outlog : boolean )
256: return Poly_Sys is
257:
258: -- DESCRIPTION :
259: -- Returns the moving part for the case i > 0 and jmp < p,
260: -- for the node from the second Pieri tree.
261:
262: -- REQUIRED : nd.c > 0.
263:
264: lc : constant natural := l(nd.c)'length(2);
265: movU : Standard_Complex_Poly_Matrices.Matrix(x'range(1),1..lc)
266: := Moving_U_Matrix(file,outlog,nd,lb,f,l);
267: movUsec : constant Standard_Complex_Poly_Matrices.Matrix
268: := Upper_Section(movU,lb(jmp));
269: res : constant Poly_Sys := Moving_Cycle(movUsec,x);
270:
271: begin
272: if outlog
273: then put_line(file,"The right moving U matrix : "); put(file,movU);
274: put_line(file,"The right moving cycle : "); put(file,movUsec);
275: put_line(file,"The equations of the moving cycle :"); put(file,res);
276: end if;
277: Standard_Complex_Poly_Matrices.Clear(movU);
278: return res;
279: end Right_Moving_Cycle;
280:
281: procedure Moving_Cycles
282: ( file : in file_type; pnd : in Paired_Nodes; id : in natural;
283: jmp1,jmp2 : in natural; b1,b2 : in Bracket;
284: f1,f2 : in Standard_Complex_Matrices.Matrix;
285: l1,l2 : in VecMat; ln : in Matrix; report,outlog : in boolean;
286: sol : in out Matrix ) is
287:
288: -- DESCRIPTION :
289: -- Set up the moving cycles for the current pair of nodes down to
290: -- the nodes (b1,b2), jumping along the indices (jmp1,jmp2).
291: -- The other parameters are the same as in Deform_Pair.
292:
293: cb1 : constant Bracket := Coordinate_Bracket(pnd.left.all,jmp1);
294: cb2 : constant Bracket := Coordinate_Bracket(pnd.right.all,jmp2);
295: xpm : Standard_Complex_Poly_Matrices.Matrix(sol'range(1),sol'range(2))
296: := Schubert_Pattern(sol'last(1),cb1,cb2);
297: homotopy : Link_to_Poly_Sys;
298: locmap : constant Standard_Natural_Matrices.Matrix
299: := Standard_Coordinate_Frame(xpm,sol);
300:
301: begin
302: if outlog
303: then
304: put_line(file,"The localization map of the solution at the pair :");
305: put(file,xpm);
306: end if;
307: homotopy := new Poly_Sys'(Polynomial_Equations(ln,xpm));
308: for i in pnd.left.c+1..l1'last loop
309: Concat(homotopy,Polynomial_Equations(l1(i).all,xpm));
310: end loop;
311: for i in pnd.right.c+1..l2'last loop
312: Concat(homotopy,Polynomial_Equations(l2(i).all,xpm));
313: end loop;
314: for i in homotopy'range loop
315: Embed(homotopy(i));
316: end loop;
317: if not Is_Leaf(pnd.left.all) and (pnd.left.c > 0)
318: then
319: if pnd.left.i = 0
320: then Concat(homotopy,
321: One_Moving_Cycle(file,pnd.left.all,b1,f1,l1,xpm,outlog));
322: elsif jmp1 < b1'last
323: then Concat(homotopy,
324: Left_Moving_Cycle
325: (file,pnd.left.all,b1,jmp1,f1,l1,xpm,outlog));
326: end if;
327: end if;
328: if not Is_Leaf(pnd.right.all) and (pnd.right.c > 0)
329: then
330: if pnd.right.i = 0
331: then Concat(homotopy,
332: One_Moving_Cycle(file,pnd.right.all,b2,f2,l2,xpm,outlog));
333: elsif jmp2 < b2'last
334: then
335: Concat(homotopy,
336: Right_Moving_Cycle
337: (file,pnd.right.all,b2,jmp2,f2,l2,xpm,outlog));
338: end if;
339: end if;
340: if outlog
341: then put_line(file,"All equations in the homotopy :");
342: put_line(file,homotopy.all);
343: end if;
344: Trace_Paths(file,homotopy.all,locmap,report,outlog,sol);
345: Test_Solution(file,pnd.left.all,pnd.right.all,id,l1,l2,ln,xpm,sol);
346: Standard_Complex_Poly_Matrices.Clear(xpm);
347: Clear(homotopy);
348: end Moving_Cycles;
349:
350: -- TARGET PROCEDURES :
351:
352: procedure Deform_Pair
353: ( file : in file_type; pnd : in Paired_Nodes; id : in natural;
354: f1,f2 : in Standard_Complex_Matrices.Matrix;
355: l1,l2 : in VecMat; ln : in Matrix; report,outlog : in boolean;
356: sol : in out Matrix ) is
357:
358: down : Paired_Nodes;
359: jmp1 : constant natural := Jump(pnd.left.all);
360: jmp2 : constant natural := Jump(pnd.right.all);
361: lb1 : constant Bracket := Lower_Jump_Decrease(pnd.left.all);
362: lb2 : constant Bracket := Lower_Jump_Decrease(pnd.right.all);
363:
364: begin
365: put(file,"JUMP @(");
366: put(file,pnd.left.all); put(file,",");
367: put(file,pnd.right.all); put(file,")"); put(file," = (");
368: put(file,jmp1,1); put(file,","); put(file,jmp2,1); put(file,") -> (");
369: put(file,lb1); put(file,","); put(file,lb2); put_line(file,").");
370: if (Is_Leaf(pnd.left.all) and Is_Leaf(pnd.right.all))
371: then Start_at_Leaves(file,pnd,ln,sol);
372: else Moving_Cycles
373: (file,pnd,id,jmp1,jmp2,lb1,lb2,f1,f2,l1,l2,ln,report,outlog,sol);
374: end if;
375: if not At_First_Branch_Point(pnd)
376: then down := Ancestor(pnd);
377: Deform_Pair(file,down,id,f1,f2,l1,l2,ln,report,outlog,sol);
378: end if;
379: -- if pnd.left.h = pnd.right.h
380: -- then if (((pnd.left.c <= 1) and (pnd.right.c <=1))
381: -- and then (((pnd.left.i = 0) and (pnd.left.c = 1))
382: -- or else ((pnd.right.i = 0) and (pnd.right.c = 1))))
383: -- then null;
384: -- else down.left := pnd.left.ancestor;
385: -- down.right := pnd.right.ancestor;
386: -- Deform_Pair(file,down,id,f1,f2,l1,l2,ln,report,outlog,sol);
387: -- end if;
388: -- elsif pnd.left.h > pnd.right.h
389: -- then down.left := pnd.left.ancestor;
390: -- down.right := pnd.right;
391: -- Deform_Pair(file,down,id,f1,f2,l1,l2,ln,report,outlog,sol);
392: -- else down.left := pnd.left;
393: -- down.right := pnd.right.ancestor;
394: -- Deform_Pair(file,down,id,f1,f2,l1,l2,ln,report,outlog,sol);
395: -- end if;
396: end Deform_Pair;
397:
398: procedure Write_Paired_Chains ( file : in file_type;
399: pnd : in Paired_Nodes; ind : in natural ) is
400:
401: -- DESCRIPTION :
402: -- Writes the chains downwards that start at the leaves of pnd.
403: -- The paired nodes have index number ind.
404:
405: begin
406: put(file,"DESCENDING THE PAIRED CHAINS "); put(file,ind,1);
407: put_line(file," :");
408: put(file,pnd.left); new_line(file);
409: put(file,pnd.right); new_line(file);
410: end Write_Paired_Chains;
411:
412: procedure Deform_Pairs
413: ( file : in file_type; n,d : in natural;
414: lp : in List_of_Paired_Nodes; l1,l2 : in VecMat;
415: ln : in Matrix; report,outlog : in boolean;
416: sols : out VecMat ) is
417:
418: tmp : List_of_Paired_Nodes := lp;
419: pnd : Paired_Nodes;
420: f1 : constant Standard_Complex_Matrices.Matrix
421: := Random_Upper_Triangular(n);
422: f2 : constant Standard_Complex_Matrices.Matrix
423: := Random_Lower_Triangular(n);
424: firstpair : Paired_Nodes := Head_Of(lp);
425: lb1 : constant Bracket := Lowest_Jump_Decrease(firstpair.left.all);
426: lb2 : constant Bracket := Lowest_Jump_Decrease(firstpair.right.all);
427: xpm : Standard_Complex_Poly_Matrices.Matrix(ln'range(1),lb1'range)
428: := Schubert_Pattern(ln'last(1),lb1,lb2);
429:
430: begin
431: for i in 1..Length_Of(lp) loop
432: pnd := Head_Of(tmp);
433: Write_Paired_Chains(file,pnd,i);
434: sols(i) := new Standard_Complex_Matrices.Matrix(1..n,1..d);
435: Deform_Pair(file,pnd,i,f1,f2,l1,l2,ln,report,outlog,sols(i).all);
436: tmp := Tail_Of(tmp);
437: end loop;
438: Test_Solutions(file,l1,l2,ln,xpm,sols);
439: end Deform_Pairs;
440:
441: end Pieri_Deformations;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to pieri_deformations.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