Annotation of OpenXM_contrib/PHC/Ada/Schubert/driver_for_pieri_homotopies.adb, Revision 1.1.1.1
1.1 maekawa 1: with text_io,integer_io; use text_io,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_Natural_Vectors;
6: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
7: with Standard_Complex_Vectors;
8: with Standard_Complex_Vectors_io; use Standard_Complex_Vectors_io;
9: with Standard_Natural_Matrices;
10: with Standard_Natural_Matrices_io; use Standard_Natural_Matrices_io;
11: with Standard_Floating_Matrices;
12: with Standard_Complex_Matrices;
13: with Standard_Complex_Matrices_io; use Standard_Complex_Matrices_io;
14: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
15: with Standard_Random_Numbers; use Standard_Random_Numbers;
16: with Standard_Random_Vectors; use Standard_Random_Vectors;
17: with Standard_Matrix_Inversion; use Standard_Matrix_Inversion;
18: with Standard_Complex_VecMats; use Standard_Complex_VecMats;
19: with Symbol_Table; use Symbol_Table;
20: with Matrix_Indeterminates;
21: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
22: with Standard_Complex_Poly_Matrices;
23: with Standard_Complex_Poly_Matrices_io; use Standard_Complex_Poly_Matrices_io;
24: with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
25: with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
26: with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
27: with Homotopy;
28: with Standard_Complex_Solutions; use Standard_Complex_Solutions;
29: with Standard_Complex_Solutions_io; use Standard_Complex_Solutions_io;
30: with Standard_Root_Refiners; use Standard_Root_Refiners;
31: with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
32: with Drivers_for_Poly_Continuation; use Drivers_for_Poly_Continuation;
33: with Brackets,Brackets_io; use Brackets,Brackets_io;
34: with Bracket_Monomials,Bracket_Systems; use Bracket_Monomials,Bracket_Systems;
35: with Specialization_of_Planes; use Specialization_of_Planes;
36: with Symbolic_Minor_Equations; use Symbolic_Minor_Equations;
37: with Numeric_Minor_Equations; use Numeric_Minor_Equations;
38: with Determinantal_Systems; use Determinantal_Systems;
39: with Plane_Representations; use Plane_Representations;
40: with Localization_Posets; use Localization_Posets;
41: with Localization_Posets_io; use Localization_Posets_io;
42: with Deformation_Posets; use Deformation_Posets;
43: with Drivers_for_Input_Planes; use Drivers_for_Input_Planes;
44:
45: procedure Driver_for_Pieri_Homotopies
46: ( file : in file_type; n,d : in natural ) is
47:
48: -- AUXILIARIES IN THE SOLVERS :
49:
50: procedure Add_t_Symbol is
51:
52: -- DESCRIPTION :
53: -- Adds the symbol for the continuation parameter t to the symbol table.
54:
55: tsb : Symbol;
56:
57: begin
58: Symbol_Table.Enlarge(1);
59: tsb(1) := 't';
60: for i in 2..tsb'last loop
61: tsb(i) := ' ';
62: end loop;
63: Symbol_Table.Add(tsb);
64: end Add_t_Symbol;
65:
66: procedure Set_Parameters ( file : in file_type; report : out boolean ) is
67:
68: -- DESCRIPTION :
69: -- Interactive determination of the continuation and output parameters.
70:
71: oc : natural;
72:
73: begin
74: new_line;
75: Driver_for_Continuation_Parameters(file);
76: new_line;
77: Driver_for_Process_io(file,oc);
78: report := not (oc = 0);
79: new_line;
80: put_line("No more input expected. See output file for results...");
81: new_line;
82: new_line(file);
83: end Set_Parameters;
84:
85: function First_Standard_Plane
86: ( m,p : natural ) return Standard_Complex_Matrices.Matrix is
87:
88: -- DESCRIPTION :
89: -- Returns the m-plane spanned by the first m standard basis vectors.
90:
91: res : Standard_Complex_Matrices.Matrix(1..m+p,1..m);
92:
93: begin
94: for j in 1..m loop -- assign j-th column
95: for i in 1..m+p loop
96: res(i,j) := Create(0.0); -- initialize with zeros
97: end loop;
98: res(j,j) := Create(1.0); -- j-th column = j-th basis vector
99: end loop;
100: return res;
101: end First_Standard_Plane;
102:
103: function Last_Standard_Plane
104: ( m,p : natural ) return Standard_Complex_Matrices.Matrix is
105:
106: -- DESCRIPTION :
107: -- Returns the m-plane spanned by the last m standard basis vectors.
108:
109: res : Standard_Complex_Matrices.Matrix(1..m+p,1..m);
110:
111: begin
112: for j in 1..m loop -- assign j-th column
113: for i in 1..m+p loop
114: res(i,j) := Create(0.0); -- initialize with zeros
115: end loop;
116: res(p+j,j) := Create(1.0); -- j-th vector = (p+j)-th basis vector
117: end loop;
118: return res;
119: end Last_Standard_Plane;
120:
121: procedure Basis_Change ( n : in natural; vm : in VecMat; transvm : out VecMat;
122: trans : out Standard_Complex_Matrices.Matrix ) is
123:
124: -- DESCRIPTION :
125: -- Changes basis so that the last two planes in vm are spanned by
126: -- the first and last standard basis vectors respectively.
127:
128: wrk : Standard_Complex_Matrices.Matrix(1..n,1..n);
129: -- penuplane : constant Standard_Complex_Matrices.Matrix
130: -- := vm(vm'last-1).all;
131: frstplane : constant Standard_Complex_Matrices.Matrix := vm(vm'first).all;
132: lastplane : constant Standard_Complex_Matrices.Matrix := vm(vm'last).all;
133: colind : natural := 0;
134: use Standard_Complex_Matrices;
135: ranflt : double_float;
136:
137: begin
138: for j in frstplane'range(2) loop -- fill up the work matrix
139: colind := colind + 1;
140: for i in frstplane'range(1) loop
141: wrk(i,colind) := frstplane(i,j);
142: end loop;
143: end loop;
144: for j in colind+1..(n-lastplane'length(2)) loop -- random spacing
145: for i in 1..n loop
146: ranflt := Random;
147: wrk(i,j) := Create(ranflt);
148: end loop;
149: end loop;
150: colind := n+1;
151: for j in reverse lastplane'range(2) loop
152: colind := colind - 1;
153: for i in lastplane'range(1) loop
154: wrk(i,colind) := lastplane(i,j);
155: end loop;
156: end loop;
157: trans := Inverse(wrk); -- transformation = inverse
158: for i in vm'first+1..vm'last-1 loop
159: transvm(i-1) := new Standard_Complex_Matrices.Matrix'
160: (trans*vm(i).all);
161: end loop;
162: transvm(transvm'last-1)
163: := new Standard_Complex_Matrices.Matrix'(trans*vm(vm'first).all);
164: transvm(transvm'last)
165: := new Standard_Complex_Matrices.Matrix'(trans*vm(vm'last).all);
166: end Basis_Change;
167:
168: function Solution_Plane ( locmap : Standard_Natural_Matrices.Matrix;
169: mat : Standard_Complex_Matrices.Matrix )
170: return Solution is
171:
172: -- DESCRIPTION :
173: -- Returns the representation of the solution plane as a vector.
174:
175: solloc : constant Standard_Complex_Vectors.Vector
176: := Vector_Rep(locmap,Localize(locmap,mat));
177: sol : Solution(solloc'length);
178:
179: begin
180: sol.m := 1;
181: sol.t := Create(0.0);
182: sol.res := 0.0;
183: sol.err := 0.0;
184: sol.rco := 0.0;
185: sol.v := solloc;
186: return sol;
187: end Solution_Plane;
188:
189: function Solution_Planes ( locmap : Standard_Natural_Matrices.Matrix;
190: vm : VecMat ) return Solution_List is
191:
192: -- DESCRIPTION :
193: -- Returns the representation of the vector of planes as a solution list.
194:
195: res,res_last : Solution_List;
196:
197: begin
198: for i in vm'range loop
199: Append(res,res_last,Solution_Plane(locmap,vm(i).all));
200: end loop;
201: return res;
202: end Solution_Planes;
203:
204: function Create_Hypersurface_System
205: ( n : natural; locmap : Standard_Natural_Matrices.Matrix;
206: xpm : Standard_Complex_Poly_Matrices.Matrix; planes : VecMat )
207: return Poly_Sys is
208:
209: -- DESCRIPTION :
210: -- Returns the polynomial system that collects the hypersurface
211: -- intersection conditions for meeting the given m-planes.
212: -- The system is localized according to the given localization map.
213:
214: wrksys : Poly_Sys(1..n) := Polynomial_Equations(planes(1..n),xpm);
215: res : Poly_Sys(1..n) := Localize(locmap,wrksys);
216:
217: begin
218: Clear(wrksys);
219: return res;
220: end Create_Hypersurface_System;
221:
222: function Create_General_System
223: ( locmap : Standard_Natural_Matrices.Matrix;
224: xpm : Standard_Complex_Poly_Matrices.Matrix; planes : VecMat )
225: return Poly_Sys is
226:
227: -- DESCRIPTION :
228: -- Returns the polynomial system that collects the general
229: -- intersection conditions for meeting the given (m+1-k(i))-planes.
230: -- The system is localized according to the given localization map.
231:
232: wrksys : constant Poly_Sys := Polynomial_Equations(planes,xpm);
233: res : constant Poly_Sys := Localize(locmap,wrksys);
234: tmp : Poly;
235:
236: begin
237: for i in wrksys'range loop
238: tmp := wrksys(i);
239: Clear(tmp);
240: end loop;
241: return res;
242: end Create_General_System;
243:
244: procedure Evaluate_Roots ( file : in file_type;
245: p : in Poly_sys; sols : in out Solution_List ) is
246:
247: -- DESCRIPTION :
248: -- Evaluates the roots at the system, good for overdetermined systems.
249:
250: tmp : Solution_List := sols;
251:
252: begin
253: for i in 1..Length_Of(sols) loop
254: put(file,"Solution "); put(file,i,1);
255: put_line(file," evaluated at the system : ");
256: put_line(file,Eval(p,Head_Of(tmp).v));
257: tmp := Tail_Of(tmp);
258: end loop;
259: end Evaluate_Roots;
260:
261: function Square ( n : natural; p : Poly_Sys ) return Poly_Sys is
262:
263: -- DESCRIPTION :
264: -- Returns a n-by-n system, by adding random linear combinations of
265: -- the polynomials p(i), i > n, to the first n polynomials.
266:
267: res : Poly_Sys(1..n);
268: acc : Poly;
269:
270: begin
271: for i in res'range loop
272: Copy(p(i),res(i));
273: for j in n+1..p'last loop
274: acc := Random1*p(j);
275: Add(res(i),acc);
276: Clear(acc);
277: end loop;
278: end loop;
279: return res;
280: end Square;
281:
282: procedure Refine_Roots ( file : in file_type;
283: p : in Poly_Sys; sols : in out Solution_List ) is
284:
285: -- DESCRIPTION :
286: -- Calls the root refiner, for square polynomial systems only.
287:
288: epsxa : constant double_float := 10.0**(-8);
289: epsfa : constant double_float := 10.0**(-8);
290: tolsing : constant double_float := 10.0**(-8);
291: max : constant natural := 3;
292: numit : natural := 0;
293:
294: begin
295: Reporting_Root_Refiner(file,p,sols,epsxa,epsfa,tolsing,numit,max,false);
296: end Refine_Roots;
297:
298: function General_Homotopy
299: ( xpm : Standard_Complex_Poly_Matrices.Matrix;
300: locmap : Standard_Natural_Matrices.Matrix;
301: start,target : VecMat ) return Poly_Sys is
302:
303: -- DESCRIPTION :
304: -- Returns a general homotopy between start and target planes.
305: -- The continuation parameter is inside the determinants.
306:
307: -- REQUIRED : dimensions of start and target matrices must correspond.
308:
309: res : Link_to_Poly_Sys;
310: n : constant natural := xpm'length(1);
311: p : constant natural := xpm'length(2);
312: m : constant natural := n-p;
313: nva : constant natural := n*p + 1;
314:
315: begin
316: for i in start'range loop
317: declare
318: moving : Standard_Complex_Poly_Matrices.Matrix
319: (start(i)'range(1),start(i)'range(2))
320: := Moving_U_Matrix(nva,start(i).all,target(i).all);
321: kd : constant natural := p + start(i)'length(2);
322: bm : Bracket_Monomial := Maximal_Minors(n,kd);
323: bs : Bracket_System(0..Number_of_Brackets(bm))
324: := Minor_Equations(kd,kd-p,bm);
325: sys : Poly_Sys(1..bs'last)
326: := Lifted_Expanded_Minors(moving,xpm,bs);
327: begin
328: Concat(res,sys);
329: Standard_Complex_Poly_Matrices.Clear(moving);
330: end;
331: end loop;
332: declare
333: locres : constant Poly_Sys := Localize(locmap,res.all);
334: begin
335: Clear(res);
336: return locres;
337: end;
338: end General_Homotopy;
339:
340: procedure Continuation ( file : in file_type; sols : in out Solution_List;
341: report : in boolean ) is
342:
343: -- DESCRIPTION :
344: -- Calls the path trackers starting at the given solutions.
345:
346: -- REQUIRED : Homotopy is properly created.
347:
348: procedure Sil_Cont is
349: new Silent_Continue(Max_Norm,
350: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
351: procedure Rep_Cont is
352: new Reporting_Continue(Max_Norm,
353: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
354:
355: begin
356: Set_Continuation_Parameter(sols,Create(0.0));
357: if report
358: then Rep_Cont(file,sols,false,Create(1.0));
359: else Sil_Cont(sols,false,Create(1.0));
360: end if;
361: end Continuation;
362:
363: procedure General_Path_Following
364: ( file : in file_type; target,homsys : in Poly_Sys;
365: sols : in out Solution_List; report : in boolean ) is
366:
367: -- DESCRIPTION :
368: -- Given a homotopy with start solutions, the path trackers will
369: -- trace the paths defined by this homotopy.
370:
371: -- REQUIRED : not Is_Null(sols).
372:
373: timer : Timing_Widget;
374: neq : constant natural := homsys'last; -- number of equations
375: dim : constant natural := Head_Of(sols).n; -- actual dimension
376: squhom : Poly_Sys(1..dim) := Square(dim,homsys);
377:
378: begin
379: tstart(timer);
380: Homotopy.Create(squhom,dim+1);
381: Continuation(file,sols,report);
382: tstop(timer);
383: new_line(file);
384: print_times(file,timer,"Determinantal Cheater's homotopy to target system");
385: new_line(file);
386: Evaluate_Roots(file,target,sols);
387: end General_Path_Following;
388:
389: procedure Path_Following_with_Cheater
390: ( file : in file_type; start,target : in Poly_Sys;
391: sols : in out Solution_List; report : in boolean ) is
392:
393: -- DESCRIPTION :
394: -- Calls the standard continuation routines to solve a specific
395: -- target system, starting at the solutions of the start system.
396: -- This is the usual linear cheater between start and target system,
397: -- although it will take nonsquare inputs, it should only be used
398: -- for hypersurface intersection conditions.
399:
400: -- REQUIRED : not Is_Null(sols).
401:
402: timer : Timing_Widget;
403: n : constant natural := target'last;
404: a : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(1,n);
405: b : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(1,n);
406: dimsol : constant natural := Head_Of(sols).n;
407: squ_target,squ_start : Poly_Sys(1..dimsol);
408:
409: begin
410: tstart(timer);
411: if dimsol = n
412: then Homotopy.Create(target,start,2,a,b,true); -- linear cheater
413: else squ_target := Square(dimsol,target);
414: squ_start := Square(dimsol,start);
415: Homotopy.Create(squ_target,squ_start,2,a,b,true);
416: end if;
417: Continuation(file,sols,report);
418: tstop(timer);
419: new_line(file);
420: print_times(file,timer,"Linear Cheater's homotopy to target system");
421: if dimsol = n
422: then Refine_Roots(file,target,sols);
423: else Refine_Roots(file,squ_target,sols);
424: Evaluate_Roots(file,target,sols);
425: end if;
426: end Path_Following_with_Cheater;
427:
428: procedure Solve_Hypersurface_Target_System
429: ( file : in file_type; m,p : in natural;
430: start_planes,target_planes : in VecMat;
431: index_poset : in Array_of_Array_of_Nodes;
432: deform_poset : in Array_of_Array_of_VecMats;
433: target_level : in natural; report : in boolean ) is
434:
435: -- DESCRIPTION :
436: -- This procedure tests the output of the deformation poset,
437: -- creating polynomial representations of the intersection conditions
438: -- and solution lists representing the solution planes.
439:
440: -- ON ENTRY :
441: -- file to write intermediate output on;
442: -- m dimension of the input planes;
443: -- p dimension of the solution planes;
444: -- start_planes input m-planes in general position;
445: -- target_planes specific input m-planes;
446: -- index_poset indexed localization poset;
447: -- deform_poset poset with the solution p-planes;
448: -- target_level indicates lowest node in deform_poset that is filled;
449: -- report switch for intermediate output during continuation.
450:
451: dim : constant natural := m*p;
452: xpm : Standard_Complex_Poly_Matrices.Matrix(1..m+p,1..p)
453: := Localization_Pattern(m+p,index_poset(dim)(1).top,
454: index_poset(dim)(1).bottom);
455: solplanes : constant VecMat := deform_poset(target_level)(1).all;
456: locmap : Standard_Natural_Matrices.Matrix(1..m+p,1..p)
457: := Standard_Coordinate_Frame(xpm,solplanes(1).all);
458: locsys : Poly_Sys(start_planes'range)
459: := Create_Hypersurface_System(dim,locmap,xpm,start_planes);
460: target : Poly_Sys(target_planes'range)
461: := Create_Hypersurface_System(dim,locmap,xpm,target_planes);
462: sols : Solution_List := Solution_Planes(locmap,solplanes);
463:
464: begin
465: Matrix_Indeterminates.Reduce_Symbols(locmap);
466: put_line(file,"THE GENERIC SYSTEM : ");
467: put_line(file,locsys);
468: new_line(file);
469: put_line(file,"The localization map :"); put(file,locmap);
470: Refine_Roots(file,locsys,sols);
471: new_line(file);
472: put_line(file,"THE TARGET SYSTEM : ");
473: put_line(file,target);
474: new_line(file);
475: Path_Following_with_Cheater(file,locsys,target,sols,report);
476: end Solve_Hypersurface_Target_System;
477:
478: procedure Solve_General_Target_System
479: ( file : in file_type; m,p : in natural; k : in Bracket;
480: start_planes,target_planes : in VecMat;
481: index_poset : in Array_of_Array_of_Nodes;
482: deform_poset : in Array_of_Array_of_VecMats;
483: target_level : in natural; report : in boolean ) is
484:
485: -- DESCRIPTION :
486: -- This procedure tests the output of the deformation poset,
487: -- creating polynomial representations of the intersection conditions
488: -- and solution lists representing the solution planes.
489:
490: -- ON ENTRY :
491: -- file to write intermediate output on;
492: -- m dimension of the input planes;
493: -- p dimension of the solution planes;
494: -- k co-dimension conditions;
495: -- start_planes input (m+1-k(i))-planes in general position;
496: -- target_planes specific input (m+1-k(i))-planes;
497: -- index_poset indexed localization poset;
498: -- deform_poset poset with the solution p-planes;
499: -- target_level indicates lowest node in deform_poset that is filled;
500: -- report switch for intermediate output during continuation.
501:
502: dim : constant natural := m*p;
503: xpm : Standard_Complex_Poly_Matrices.Matrix(1..m+p,1..p)
504: := Localization_Pattern(m+p,index_poset(dim)(1).top,
505: index_poset(dim)(1).bottom);
506: solplanes : constant VecMat := deform_poset(target_level)(1).all;
507: locmap : Standard_Natural_Matrices.Matrix(1..m+p,1..p)
508: := Standard_Coordinate_Frame(xpm,solplanes(1).all);
509: homtpy : constant Poly_Sys
510: := General_Homotopy(xpm,locmap,start_planes,target_planes);
511: locsys : constant Poly_Sys
512: := Create_General_System(locmap,xpm,start_planes);
513: target : constant Poly_Sys
514: := Create_General_System(locmap,xpm,target_planes);
515: sols : Solution_List := Solution_Planes(locmap,solplanes);
516:
517: begin
518: Matrix_Indeterminates.Reduce_Symbols(locmap);
519: put_line(file,"THE GENERIC SYSTEM : ");
520: put_line(file,locsys);
521: new_line(file);
522: put_line(file,"The localization map :"); put(file,locmap);
523: Evaluate_Roots(file,locsys,sols);
524: new_line(file);
525: put_line(file,"THE TARGET SYSTEM : ");
526: put_line(file,target);
527: new_line(file);
528: General_Path_Following(file,target,homtpy,sols,report);
529: end Solve_General_Target_System;
530:
531: procedure Write_Poset_Times
532: ( file : in file_type; timer : in Timing_Widget;
533: npaths : in Standard_Natural_Vectors.Vector;
534: timings : in Duration_Array ) is
535:
536: -- DESCRIPTION :
537: -- Writes a overview of #paths and timings spent during deformations
538: -- along the poset structure.
539:
540: begin
541: new_line(file);
542: put_line(file,"--------------------------------------");
543: put_line(file,"| TIMING INFORMATION OVERVIEW |");
544: put_line(file,"--------------------------------------");
545: put_line(file,"| n | #paths | user cpu time |");
546: put_line(file,"--------------------------------------");
547: for i in npaths'range loop
548: if npaths(i) /= 0
549: then put(file,"|"); put(file,i,4); put(file," |");
550: put(file,npaths(i),7); put(file," | ");
551: print_hms(file,timings(i)); put(file," |"); new_line(file);
552: end if;
553: end loop;
554: put_line(file,"--------------------------------------");
555: put(file,"| total |");
556: put(file,Standard_Natural_Vectors.Sum(npaths),7);
557: put(file," | ");
558: print_hms(file,Elapsed_User_Time(timer));
559: put(file," |"); new_line(file);
560: put_line(file,"--------------------------------------");
561: end Write_Poset_Times;
562:
563: procedure Solve_Hypersurface_Deformation_Poset
564: ( file : in file_type; m,p : in natural;
565: level_poset : in Array_of_Nodes;
566: index_poset : in Array_of_Array_of_Nodes ) is
567:
568: -- DESCRIPTION :
569: -- Creates a deformation poset and applies the Solve operator.
570:
571: n : constant natural := m+p;
572: deform_poset : Array_of_Array_of_VecMats(index_poset'range)
573: := Create(index_poset);
574: planes : VecMat(1..m*p) := Random_Complex_Planes(m,p);
575: target_planes : VecMat(1..m*p);
576: ans : character;
577: report,outlog : boolean;
578: timer : Timing_Widget;
579: target_level : natural := m*p;
580: nbp : natural := 0;
581: npaths : Standard_Natural_Vectors.Vector(1..target_level)
582: := (1..target_level => 0);
583: timings : Duration_Array(1..target_level) := (1..target_level => 0.0);
584:
585: begin
586: put_line("The size of the deformation poset : ");
587: put_line(file,"The size of the deformation poset : ");
588: put_roco(index_poset);
589: put_roco(file,index_poset);
590: new_line;
591: put("Give target level <= "); put(target_level,1);
592: put(" = root level by default : "); get(target_level);
593: for i in 1..target_level loop
594: nbp := nbp + Row_Root_Count_Sum(level_poset,i);
595: end loop;
596: put("The number of paths : "); put(nbp,1); new_line;
597: put(file,"The number of paths : "); put(file,nbp,1); new_line(file);
598: new_line;
599: put("Do you want to have the homotopies on file ? (y/n) ");
600: Ask_Yes_or_No(ans);
601: outlog := (ans = 'y');
602: if target_level < m*p
603: then planes(m*p).all := Last_Standard_Plane(m,p);
604: if target_level < m*p-1
605: then planes(m*p-1).all := First_Standard_Plane(m,p);
606: end if;
607: end if;
608: Driver_for_Input_Planes(file,m,p,target_planes);
609: Matrix_Indeterminates.Initialize_Symbols(m+p,p);
610: Add_t_Symbol;
611: Set_Parameters(file,report);
612: tstart(timer);
613: for i in index_poset(target_level)'range loop
614: declare
615: root : Node := index_poset(target_level)(i).all;
616: begin
617: Solve(file,m+p,deform_poset,root,planes(1..target_level),
618: report,outlog,npaths,timings);
619: end;
620: end loop;
621: tstop(timer);
622: new_line(file);
623: print_times(file,timer,"Solving along the deformation poset");
624: Write_Poset_Times(file,timer,
625: npaths(1..target_level),timings(1..target_level));
626: if deform_poset(target_level)(1) /= null
627: then new_line(file);
628: Solve_Hypersurface_Target_System
629: (file,m,p,planes,target_planes,index_poset,deform_poset,
630: target_level,report);
631: end if;
632: end Solve_Hypersurface_Deformation_Poset;
633:
634: procedure Solve_General_Deformation_Poset
635: ( file : in file_type; m,p : in natural; k : in Bracket;
636: index_poset : in Array_of_Array_of_Nodes ) is
637:
638: -- DESCRIPTION :
639: -- Applies the solver to general intersection conditions.
640:
641: deform_poset : Array_of_Array_of_VecMats(index_poset'range)
642: := Create(index_poset);
643: planes : VecMat(k'range) := Random_Complex_Planes(m,p,k);
644: target_planes : VecMat(k'range);
645: ans : character;
646: report,outlog : boolean;
647: timer : Timing_Widget;
648: target_level : natural := m*p;
649: npaths : Standard_Natural_Vectors.Vector(1..target_level)
650: := (1..target_level => 0);
651: timings : Duration_Array(1..target_level) := (1..target_level => 0.0);
652:
653: begin
654: put_line("The size of the deformation poset : ");
655: put_line(file,"The size of the deformation poset : ");
656: put_roco(index_poset);
657: put_roco(file,index_poset);
658: new_line;
659: put("Give target level <= "); put(target_level,1);
660: put(" = root level by default : "); get(target_level);
661: new_line;
662: put("Do you want to have the homotopies on file ? (y/n) ");
663: Ask_Yes_or_No(ans);
664: outlog := (ans = 'y');
665: Driver_for_Input_Planes(file,m,p,k,target_planes);
666: Matrix_Indeterminates.Initialize_Symbols(m+p,p);
667: Add_t_Symbol;
668: Set_Parameters(file,report);
669: tstart(timer);
670: for i in index_poset(target_level)'range loop
671: declare
672: root : Node := index_poset(target_level)(i).all;
673: begin
674: Solve(file,m+p,k,deform_poset,root,planes,report,outlog,
675: npaths,timings);
676: end;
677: end loop;
678: tstop(timer);
679: new_line(file);
680: print_times(file,timer,"Solving along the deformation poset");
681: Write_Poset_Times(file,timer,
682: npaths(1..target_level),timings(1..target_level));
683: if deform_poset(target_level)(1) /= null
684: then new_line(file);
685: Solve_General_Target_System
686: (file,m,p,k,planes,target_planes,index_poset,deform_poset,
687: target_level,report);
688: end if;
689: end Solve_General_Deformation_Poset;
690:
691: -- HYPERSURFACE SOLVERS :
692:
693: procedure Create_Top_Hypersurface_Poset
694: ( file : in file_type; m,p : in natural ) is
695:
696: -- DESCRIPTION :
697: -- Create the poset by incrementing only top pivots.
698:
699: timer : Timing_Widget;
700: root : Node(p) := Trivial_Root(m,p);
701: lnkroot : Link_to_Node := new Node'(root);
702: level_poset : Array_of_Nodes(0..m*p);
703: index_poset : Array_of_Array_of_Nodes(0..m*p);
704:
705: begin
706: tstart(timer);
707: Top_Create(lnkroot,m+p);
708: put_line("The poset created from the top : ");
709: put_line(file,"The poset created from the top : ");
710: level_poset := Create_Leveled_Poset(lnkroot);
711: Count_Roots(level_poset);
712: index_poset := Create_Indexed_Poset(level_poset);
713: put(index_poset);
714: put(file,index_poset);
715: Solve_Hypersurface_Deformation_Poset(file,m,p,level_poset,index_poset);
716: tstop(timer);
717: new_line(file);
718: print_times(file,timer,
719: "Total time for Hypersurface Pieri Homotopy Algorithm");
720: end Create_Top_Hypersurface_Poset;
721:
722: procedure Create_Bottom_Hypersurface_Poset
723: ( file : in file_type; m,p : in natural ) is
724:
725: -- DESCRIPTION :
726: -- Create the poset by decrementing only bottom pivots.
727:
728: timer : Timing_Widget;
729: root : Node(p) := Trivial_Root(m,p);
730: lnkroot : Link_to_Node := new Node'(root);
731: level_poset : Array_of_Nodes(0..m*p);
732: index_poset : Array_of_Array_of_Nodes(0..m*p);
733:
734: begin
735: tstart(timer);
736: Bottom_Create(lnkroot);
737: put_line("The poset created from the bottom : ");
738: put_line(file,"The poset created from the bottom : ");
739: level_poset := Create_Leveled_Poset(lnkroot);
740: Count_Roots(level_poset);
741: index_poset := Create_Indexed_Poset(level_poset);
742: put(index_poset);
743: put(file,index_poset);
744: Solve_Hypersurface_Deformation_Poset(file,m,p,level_poset,index_poset);
745: tstop(timer);
746: new_line(file);
747: print_times(file,timer,
748: "Total time for Hypersurface Pieri Homotopy Algorithm");
749: end Create_Bottom_Hypersurface_Poset;
750:
751: procedure Create_Mixed_Hypersurface_Poset
752: ( file : in file_type; m,p : in natural ) is
753:
754: -- DESCRIPTION :
755: -- Create the poset by incrementing top and decrementing bottom pivots.
756:
757: timer : Timing_Widget;
758: root : Node(p) := Trivial_Root(m,p);
759: lnkroot : Link_to_Node := new Node'(root);
760: level_poset : Array_of_Nodes(0..m*p);
761: index_poset : Array_of_Array_of_Nodes(0..m*p);
762:
763: begin
764: tstart(timer);
765: Top_Bottom_Create(lnkroot,m+p);
766: put_line("The poset created in a mixed fashion : ");
767: put_line(file,"The poset created in a mixed fashion : ");
768: level_poset := Create_Leveled_Poset(lnkroot);
769: Count_Roots(level_poset);
770: index_poset := Create_Indexed_Poset(level_poset);
771: put(index_poset);
772: put(file,index_poset);
773: Solve_Hypersurface_Deformation_Poset(file,m,p,level_poset,index_poset);
774: tstop(timer);
775: new_line(file);
776: print_times(file,timer,
777: "Total time for Hypersurface Pieri Homotopy Algorithm");
778: end Create_Mixed_Hypersurface_Poset;
779:
780: -- GENERAL SOLVERS :
781:
782: procedure Create_Top_General_Poset
783: ( file : in file_type; m,p : in natural ) is
784:
785: -- DESCRIPTION :
786: -- Creates a poset for counting general subspace intersections,
787: -- by consistently incrementing the top pivots.
788:
789: timer : Timing_Widget;
790: root : Node(p) := Trivial_Root(m,p);
791: lnkroot : Link_to_Node := new Node'(root);
792: codim : constant Bracket := Read_Codimensions(m,p,0);
793: level_poset : Array_of_Nodes(0..m*p);
794: index_poset : Array_of_Array_of_Nodes(0..m*p);
795:
796: begin
797: put(file," k = "); put(file,codim); new_line(file);
798: tstart(timer);
799: Top_Create(lnkroot,codim,m+p);
800: put_line("The poset created from the top : ");
801: put_line(file,"The poset created from the top : ");
802: level_poset := Create_Leveled_Poset(lnkroot);
803: Count_Roots(level_poset);
804: index_poset := Create_Indexed_Poset(level_poset);
805: put(index_poset);
806: put(file,index_poset);
807: Solve_General_Deformation_Poset(file,m,p,codim,index_poset);
808: tstop(timer);
809: new_line(file);
810: print_times(file,timer,
811: "Total time for General Pieri Homotopy Algorithm");
812: end Create_Top_General_Poset;
813:
814: procedure Create_Bottom_General_Poset
815: ( file : in file_type; m,p : in natural ) is
816:
817: -- DESCRIPTION :
818: -- Creates a poset for counting general subspace intersections,
819: -- by consistently incrementing the top pivots.
820:
821: timer : Timing_Widget;
822: root : Node(p) := Trivial_Root(m,p);
823: lnkroot : Link_to_Node := new Node'(root);
824: codim : constant Bracket := Read_Codimensions(m,p,0);
825: level_poset : Array_of_Nodes(0..m*p);
826: index_poset : Array_of_Array_of_Nodes(0..m*p);
827:
828: begin
829: put(file," k = "); put(file,codim); new_line(file);
830: tstart(timer);
831: Bottom_Create(lnkroot,codim);
832: put_line("The poset created from the bottom : ");
833: put_line(file,"The poset created from the bottom : ");
834: level_poset := Create_Leveled_Poset(lnkroot);
835: Count_Roots(level_poset);
836: index_poset := Create_Indexed_Poset(level_poset);
837: put(index_poset);
838: put(file,index_poset);
839: Solve_General_Deformation_Poset(file,m,p,codim,index_poset);
840: tstop(timer);
841: new_line(file);
842: print_times(file,timer,
843: "Total time for General Pieri Homotopy Algorithm");
844: end Create_Bottom_General_Poset;
845:
846: procedure Create_Mixed_General_Poset
847: ( file : in file_type; m,p : in natural ) is
848:
849: -- DESCRIPTION :
850: -- Creates a poset for counting general subspace intersections,
851: -- by incrementing the top and decrementing the bottom pivots.
852:
853: timer : Timing_Widget;
854: root : Node(p) := Trivial_Root(m,p);
855: lnkroot : Link_to_Node := new Node'(root);
856: codim : constant Bracket := Read_Codimensions(m,p,0);
857: level_poset : Array_of_Nodes(0..m*p);
858: index_poset : Array_of_Array_of_Nodes(0..m*p);
859:
860: begin
861: put(file," k = "); put(file,codim); new_line(file);
862: tstart(timer);
863: Top_Bottom_Create(lnkroot,codim,m+p);
864: put_line("The poset created in a mixed fashion : ");
865: put_line(file,"The poset created in a mixed fashion : ");
866: level_poset := Create_Leveled_Poset(lnkroot);
867: Count_Roots(level_poset);
868: index_poset := Create_Indexed_Poset(level_poset);
869: put(index_poset);
870: put(file,index_poset);
871: Solve_General_Deformation_Poset(file,m,p,codim,index_poset);
872: tstop(timer);
873: new_line(file);
874: print_times(file,timer,
875: "Total time for General Pieri Homotopy Algorithm");
876: end Create_Mixed_General_Poset;
877:
878: procedure Main is
879:
880: p : constant natural := d;
881: m : constant natural := n-d;
882: ans : character;
883:
884: begin
885: new_line;
886: put_line("MENU for deforming p-planes in (m+p)-space, co-dimensions k : ");
887: put_line(" 1. k_i = 1 consistently incrementing the top pivots.");
888: put_line(" 2. consistently decrementing the bottom pivots.");
889: put_line(" 3. mixed top-bottom sequence for poset creation.");
890: put_line(" 4. k_i >= 1 consistently incrementing the top pivots.");
891: put_line(" 5. consistently decrementing the bottom pivots.");
892: put_line(" 6. mixed top-bottom sequence for poset creation.");
893: put("Type 1, 2, 3, 4, 5, or 6 to choose : ");
894: Ask_Alternative(ans,"123456");
895: new_line;
896: case ans is
897: when '1' => new_line(file); Create_Top_Hypersurface_Poset(file,m,p);
898: when '2' => new_line(file); Create_Bottom_Hypersurface_Poset(file,m,p);
899: when '3' => new_line(file); Create_Mixed_Hypersurface_Poset(file,m,p);
900: when '4' => Create_Top_General_Poset(file,m,p);
901: when '5' => Create_Bottom_General_Poset(file,m,p);
902: when '6' => Create_Mixed_General_Poset(file,m,p);
903: when others => put_line("Option not recognized. Please try again.");
904: end case;
905: end Main;
906:
907: begin
908: Main;
909: end Driver_for_Pieri_Homotopies;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to driver_for_pieri_homotopies.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