Annotation of OpenXM_contrib/PHC/Ada/Schubert/pieri_continuation.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
2: with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
3: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
4: with Standard_Random_Numbers; use Standard_Random_Numbers;
5: with Standard_Natural_Vectors;
6: with Standard_Natural_Vectors_io; use Standard_Natural_Vectors_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_VecVecs; use Standard_Complex_VecVecs;
11: with Standard_Natural_Matrices_io; use Standard_Natural_Matrices_io;
12: with Standard_Complex_Matrices_io; use Standard_Complex_Matrices_io;
13: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
14: with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
15: with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
16: with Standard_Complex_Solutions; use Standard_Complex_Solutions;
17: with Homotopy;
18: with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
19: with Standard_Root_Refiners; use Standard_Root_Refiners;
20: with Determinantal_Systems; use Determinantal_Systems;
21: with Plane_Representations; use Plane_Representations;
22: with Curves_into_Grassmannian; use Curves_into_Grassmannian;
23: with Curves_into_Grassmannian_io; use Curves_into_Grassmannian_io;
24:
25: package body Pieri_Continuation is
26:
27: -- AUXILIARY FORMAT CONVERTORS :
28:
29: function Create ( v : Standard_Complex_Vectors.Vector; conpar : natural )
30: return Solution is
31:
32: -- DESCRIPTION :
33: -- Returns a solution that contains a solution with vector v in it.
34: -- The v(conpar) is the continuation parameter.
35:
36: sol : Solution(v'last-1);
37:
38: begin
39: sol.t := v(conpar);
40: sol.m := 1;
41: sol.v(v'first..conpar-1) := v(v'first..conpar-1);
42: sol.v(conpar..sol.v'last) := v(conpar+1..v'last);
43: sol.err := 0.0;
44: sol.rco := 0.0;
45: sol.res := 0.0;
46: return sol;
47: end Create;
48:
49: function Create ( v : Standard_Complex_Vectors.Vector )
50: return Solution_List is
51:
52: -- DESCRIPTION :
53: -- Returns a solution list that contains a solution with vector v in it.
54:
55: -- NOTE : only needed in the original Pieri homotopy algorithm.
56:
57: res : Solution_List;
58:
59: begin
60: Add(res,Create(v,v'last));
61: return res;
62: end Create;
63:
64: function Create ( v : Standard_Complex_VecVecs.VecVec; conpar : natural )
65: return Solution_List is
66:
67: -- DESCRIPTION :
68: -- Returns a solution list that contains the vectors in v.
69:
70: -- REQUIRED : v is not empty.
71: -- The value for the continuation parameter is in v(conpar).
72:
73: res,res_last : Solution_List;
74:
75: begin
76: for i in v'range loop
77: Append(res,res_last,Create(v(i).all,conpar));
78: end loop;
79: return res;
80: end Create;
81:
82: function Convert ( sols : Solution_List )
83: return Standard_Complex_Vectors.Vector is
84:
85: -- DESCRIPTION :
86: -- Converts the solutions in sols to the vector representation.
87:
88: -- NOTE : only needed in the original Pieri homotopy algorithm.
89:
90: ls : Link_to_Solution := Head_Of(sols);
91: res : Standard_Complex_Vectors.Vector(1..ls.n+1);
92:
93: begin
94: res(ls.v'range) := ls.v;
95: res(res'last) := ls.t;
96: return res;
97: end Convert;
98:
99: function Convert ( sol : Solution; conpar : natural )
100: return Standard_Complex_Vectors.Vector is
101:
102: -- DESCRIPTION :
103: -- Converts the solution to the vector representation, where the
104: -- continuation parameter is put into the entry indexed by conpar.
105:
106: res : Standard_Complex_Vectors.Vector(1..sol.n+1);
107:
108: begin
109: res(sol.v'first..conpar-1) := sol.v(sol.v'first..conpar-1);
110: res(conpar) := sol.t;
111: res(conpar+1..sol.n+1) := sol.v(conpar..sol.v'last);
112: return res;
113: end Convert;
114:
115: function Convert ( sols : Solution_List; conpar : natural )
116: return Standard_Complex_VecVecs.VecVec is
117:
118: -- DESCRIPTION :
119: -- Converts the solutions in sols to the vector representation.
120:
121: res : Standard_Complex_VecVecs.VecVec(1..Length_Of(sols));
122: tmp : Solution_List := sols;
123:
124: begin
125: for i in res'range loop
126: declare
127: ls : Link_to_Solution := Head_Of(tmp);
128: v : Standard_Complex_Vectors.Vector(ls.v'first..ls.v'last+1);
129: begin
130: v := Convert(ls.all,conpar);
131: res(i) := new Standard_Complex_Vectors.Vector'(v);
132: end;
133: tmp := Tail_Of(tmp);
134: end loop;
135: return res;
136: end Convert;
137:
138: function Square ( n : natural; p : Poly_Sys ) return Poly_Sys is
139:
140: -- DESCRIPTION :
141: -- Returns a n-by-n system, by adding random linear combinations of
142: -- the polynomials p(i), i > n, to the first n polynomials.
143:
144: res : Poly_Sys(1..n);
145: acc : Poly;
146:
147: begin
148: for i in res'range loop
149: Copy(p(i),res(i));
150: for j in n+1..p'last loop
151: acc := Random1*p(j);
152: Add(res(i),acc);
153: Clear(acc);
154: end loop;
155: end loop;
156: return res;
157: end Square;
158:
159: procedure Refine_Roots
160: ( file : in file_type; p : in Poly_Sys;
161: sols : in out Solution_List; report : in boolean ) is
162:
163: epsxa : constant double_float := 10.0**(-8);
164: epsfa : constant double_float := 10.0**(-8);
165: tolsing : constant double_float := 10.0**(-8);
166: max : constant natural := 3;
167: numit : natural := 0;
168:
169: begin
170: if report
171: then
172: Reporting_Root_Refiner(file,p,sols,epsxa,epsfa,tolsing,numit,max,false);
173: else
174: Silent_Root_Refiner(p,sols,epsxa,epsfa,tolsing,numit,max);
175: end if;
176: end Refine_Roots;
177:
178: procedure Determinantal_Polynomial_Continuation
179: ( file : in file_type; lochom : in Poly_Sys; -- out added
180: locsol : in out Standard_Complex_Vectors.Vector;
181: report,outlog : in boolean ) is
182:
183: -- DESCRIPTION :
184: -- Given a homotopy and solution with as last variable the continuation
185: -- parameter the polynomial continuation will launched.
186: -- This homotopy uses the determinantal expansions defined in the paper.
187:
188: -- NOTE :
189: -- This routine is only used in the original Pieri homotopy algorithm.
190:
191: sols : Solution_List := Create(locsol);
192: thesys,squsys : Poly_Sys(locsol'first..locsol'last-1);
193: square_case : boolean;
194: -- ran : Complex_Number;
195:
196: begin
197: -- for i in lochom'range loop -- added to deal with real input
198: -- ran := Random1;
199: -- Mul(lochom(i),ran);
200: -- end loop;
201: if lochom'length + 1 = locsol'length
202: then square_case := true;
203: Homotopy.Create(lochom,locsol'last);
204: else square_case := false;
205: squsys := Square(locsol'last-1,lochom);
206: if outlog
207: then put_line(file,"The homotopy as square system : ");
208: put_line(file,squsys);
209: end if;
210: Homotopy.Create(squsys,locsol'last);
211: end if;
212: declare
213: procedure Sil_Cont is
214: new Silent_Continue(Max_Norm,
215: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
216: procedure Rep_Cont is
217: new Reporting_Continue(Max_Norm,
218: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
219: begin
220: if report
221: then Rep_Cont(file,sols,false,Create(1.0));
222: else Sil_Cont(sols,false,Create(1.0));
223: end if;
224: end;
225: if square_case
226: then thesys := Eval(lochom,Create(1.0),locsol'last);
227: Refine_Roots(file,thesys,sols,report);
228: else thesys := Eval(squsys,Create(1.0),locsol'last);
229: Refine_Roots(file,thesys,sols,report);
230: Clear(squsys);
231: end if;
232: Clear(thesys);
233: locsol := Convert(sols);
234: Clear(sols);
235: Homotopy.Clear;
236: end Determinantal_Polynomial_Continuation;
237:
238: procedure Determinantal_Polynomial_Continuation
239: ( file : in file_type; lochom : in Poly_Sys;
240: conpar : in natural;
241: locsols : in out Standard_Complex_VecVecs.VecVec;
242: report,outlog : in boolean ) is
243:
244: -- DESCRIPTION :
245: -- Given a homotopy and solution with as continuation parameter
246: -- the variable indexed by conpar, the continuation will be launched.
247: -- This homotopy uses the determinantal expansions defined in the paper.
248:
249: sols : Solution_List := Create(locsols,conpar);
250: firstsol : constant Standard_Complex_Vectors.Vector
251: := locsols(locsols'first).all;
252: thesys,squsys : Poly_Sys(firstsol'first..firstsol'last-1);
253: square_case : boolean;
254:
255: begin
256: if lochom'length + 1 = firstsol'length
257: then square_case := true;
258: Homotopy.Create(lochom,conpar);
259: else square_case := false;
260: squsys := Square(firstsol'last-1,lochom);
261: if outlog
262: then put_line(file,"The homotopy as square system : ");
263: put_line(file,squsys);
264: end if;
265: Homotopy.Create(squsys,conpar);
266: end if;
267: declare
268: procedure Sil_Cont is
269: new Silent_Continue(Max_Norm,
270: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
271: procedure Rep_Cont is
272: new Reporting_Continue(Max_Norm,
273: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
274: begin
275: if report
276: then Rep_Cont(file,sols,false,Create(1.0));
277: else Sil_Cont(sols,false,Create(1.0));
278: end if;
279: end;
280: if square_case
281: then thesys := Eval(lochom,Create(1.0),conpar);
282: Refine_Roots(file,thesys,sols,report);
283: else thesys := Eval(squsys,Create(1.0),conpar);
284: Refine_Roots(file,thesys,sols,report);
285: Clear(squsys);
286: end if;
287: Clear(thesys);
288: locsols := Convert(sols,conpar);
289: Clear(sols);
290: Homotopy.Clear;
291: end Determinantal_Polynomial_Continuation;
292:
293: -- TARGET ROUTINES :
294:
295: procedure Trace_Paths ( file : in file_type; homsys : in Poly_Sys;
296: locmap : in Standard_Natural_Matrices.Matrix;
297: report,outlog : in boolean;
298: plane : in out Standard_Complex_Matrices.Matrix ) is
299:
300: -- NOTE :
301: -- This routine is only used in the original Pieri homotopy algorithm.
302:
303: plavec : constant Standard_Complex_Vectors.Vector := Vector_Rep(plane);
304: solvec : Standard_Complex_Vectors.Vector(1..plavec'last+1);
305: lochom : Poly_Sys(homsys'range) := Localize(locmap,homsys);
306: plaloc : constant Standard_Complex_Matrices.Matrix
307: := Localize(locmap,plane);
308: solloc : constant Standard_Complex_Vectors.Vector
309: := Vector_Rep(locmap,plaloc);
310: locsol : Standard_Complex_Vectors.Vector(1..solloc'last+1);
311: evahom : Standard_Complex_Vectors.Vector(homsys'range);
312: evaloc : Standard_Complex_Vectors.Vector(lochom'range);
313:
314: begin
315: if outlog
316: then put_line(file,"The localization pattern :"); put(file,locmap);
317: put_line(file,"The localized homotopy : "); put_line(file,lochom);
318: end if;
319: -- put_line(file,"The plane in local coordinates : "); put(file,plaloc,2);
320: solvec(plavec'range) := plavec;
321: solvec(solvec'last) := Create(0.0);
322: -- put_line(file,"The solution vector : "); put_line(file,solvec);
323: evahom := Eval(homsys,solvec);
324: -- put_line(file,"Evaluation of solution vector at homotopy for t=0 :");
325: -- put_line(file,evahom);
326: put(file,"Residual of start solution :");
327: put(file,Max_Norm(evahom),3); put_line(file,".");
328: -- put_line(file,"The localized vector representation of the plane at t=0 :");
329: -- put_line(file,solloc);
330: locsol(solloc'range) := solloc;
331: locsol(locsol'last) := Create(0.0);
332: evaloc := Eval(lochom,locsol);
333: put(file,"Residual of localized plane at t=0 :");
334: put(file,Max_Norm(evaloc),3); put_line(file,".");
335: -- Linear_Polynomial_Continuation(file,lochom,locsol,report,outlog);
336: Determinantal_Polynomial_Continuation(file,lochom,locsol,report,outlog);
337: -- put_line(file,"The localized vector representation of the plane at t=1 :");
338: -- put_line(file,locsol);
339: evaloc := Eval(lochom,locsol);
340: -- put_line(file,"Evaluation of localized vector at homotopy for t=1 :");
341: -- put_line(file,evaloc);
342: put(file,"Residual of localized plane at t=1 :");
343: put(file,Max_Norm(evaloc),3); put_line(file,".");
344: plane := Matrix_Rep(locmap,locsol(1..locsol'last-1));
345: -- put_line(file,"The solution plane at t=1 :");
346: -- put(file,plane,2);
347: solvec(plavec'range) := Vector_Rep(plane);
348: solvec(solvec'last) := Create(1.0);
349: evahom := Eval(homsys,solvec);
350: -- put_line(file,"Evaluation of solution vector at homotopy for t=1 :");
351: -- put_line(file,evahom);
352: put(file,"Residual of target solution :");
353: put(file,Max_Norm(evahom),3); put_line(file,".");
354: Clear(lochom);
355: end Trace_Paths;
356:
357: procedure Trace_Paths
358: ( file : in file_type; homsys : in Poly_Sys;
359: locmap : in Standard_Natural_Matrices.Matrix;
360: report,outlog : in boolean;
361: planes : in Standard_Complex_VecMats.VecMat ) is
362:
363: lochom : Poly_Sys(homsys'range) := Localize(locmap,homsys);
364: locsols : Standard_Complex_VecVecs.VecVec(planes'range);
365: evaloc : Standard_Complex_Vectors.Vector(lochom'range);
366: evahom : Standard_Complex_Vectors.Vector(homsys'range);
367: lastvar : natural;
368:
369: begin
370: if outlog
371: then put_line(file,"The localized homotopy :");
372: put_line(file,lochom);
373: end if;
374: for i in planes'range loop -- solution planes into vectors
375: declare
376: plaloc : constant Standard_Complex_Matrices.Matrix
377: := Localize(locmap,planes(i).all);
378: solloc : constant Standard_Complex_Vectors.Vector
379: := Vector_Rep(locmap,plaloc);
380: locsol : Standard_Complex_Vectors.Vector(1..solloc'last+1);
381: begin
382: locsol(solloc'range) := solloc;
383: locsol(locsol'last) := Create(0.0);
384: locsols(i) := new Standard_Complex_Vectors.Vector'(locsol);
385: if outlog
386: then evaloc := Eval(lochom,locsol);
387: put(file,"Residual of localized plane at start :");
388: put(file,Max_Norm(evaloc),3); put_line(file,".");
389: end if;
390: end;
391: end loop;
392: lastvar := locsols(locsols'first)'last;
393: Determinantal_Polynomial_Continuation
394: (file,lochom,lastvar,locsols,report,outlog);
395: for i in locsols'range loop -- solution vectors into planes
396: planes(i).all := Matrix_Rep(locmap,locsols(i)(1..locsols(i)'last-1));
397: if outlog
398: then
399: evaloc := Eval(lochom,locsols(i).all);
400: put(file,"Residual of localized plane at target :");
401: put(file,Max_Norm(evaloc),3); put_line(file,".");
402: declare
403: plavec : constant Standard_Complex_Vectors.Vector
404: := Vector_Rep(planes(i).all);
405: solvec : Standard_Complex_Vectors.Vector
406: (plavec'first..plavec'last+1);
407: begin
408: solvec(plavec'range) := plavec;
409: solvec(solvec'last) := Create(1.0);
410: evahom := Eval(homsys,solvec);
411: put(file,"Residual of target solution :");
412: put(file,Max_Norm(evahom),3); put_line(file,".");
413: end;
414: end if;
415: end loop;
416: Clear(lochom); Clear(locsols);
417: end Trace_Paths;
418:
419: procedure Quantum_Trace_Paths
420: ( file : in file_type; m,p,q : in natural; nd : in Node;
421: -- (m,p,q) are needed for the symbol table
422: homsys : in Poly_Sys; conpar,s_mode : in natural;
423: locmap : in Standard_Natural_Matrices.Matrix;
424: report,outlog : in boolean;
425: planes : in Standard_Complex_VecMats.VecMat ) is
426:
427: lochom : Poly_Sys(homsys'range)
428: := Column_Localize(nd.top,nd.bottom,locmap,homsys);
429: locsols : Standard_Complex_VecVecs.VecVec(planes'range);
430: evaloc : Standard_Complex_Vectors.Vector(lochom'range);
431: evahom : Standard_Complex_Vectors.Vector(homsys'range);
432: addmix : natural;
433:
434: begin
435: if outlog
436: then put_line(file,"The localization map : "); put(file,locmap);
437: if nd.tp = mixed
438: then Two_Set_up_Symbol_Table(m,p,q,nd.top,nd.bottom);
439: else One_Set_up_Symbol_Table(m,p,q,nd.top,nd.bottom);
440: end if;
441: Reduce_Symbols(nd.top,nd.bottom,locmap);
442: put_line(file,"The localized homotopy : "); put_line(file,lochom);
443: end if;
444: if nd.tp = mixed
445: then addmix := 1;
446: else addmix := 0;
447: end if;
448: for i in planes'range loop -- solution planes into vectors
449: declare
450: plaloc : constant Standard_Complex_Matrices.Matrix
451: := Localize(locmap,planes(i).all);
452: solloc : constant Standard_Complex_Vectors.Vector
453: := Column_Vector_Rep(locmap,plaloc);
454: locsol : Standard_Complex_Vectors.Vector(1..solloc'last+2+addmix);
455: begin
456: locsol(solloc'range) := solloc;
457: if s_mode = 0
458: then locsol(locsol'last-1) := Create(0.0); -- start value for s is 0
459: else locsol(locsol'last-1) := Create(1.0); -- start value for s is 1
460: end if;
461: locsol(locsol'last) := Create(0.0); -- start value for t is 0
462: if addmix = 1
463: then locsol(locsol'last-2) := Create(1.0); -- additional s-value
464: end if;
465: locsols(i) := new Standard_Complex_Vectors.Vector'(locsol);
466: if outlog
467: then evaloc := Eval(lochom,locsol);
468: put(file,"Residual of localized plane at start :");
469: put(file,Max_Norm(evaloc),3); put_line(file,".");
470: end if;
471: end;
472: end loop;
473: Determinantal_Polynomial_Continuation
474: (file,lochom,conpar,locsols,report,outlog);
475: for i in locsols'range loop -- solution vectors into planes
476: planes(i).all
477: := Column_Matrix_Rep(locmap,locsols(i)(1..locsols(i)'last-2-addmix));
478: if outlog
479: then evaloc := Eval(lochom,locsols(i).all);
480: put(file,"Residual of localized plane at target :");
481: put(file,Max_Norm(evaloc),3); put_line(file,".");
482: put_line(file,"The computed plane : ");
483: put(file,planes(i).all,2);
484: end if;
485: declare
486: plavec : constant Standard_Complex_Vectors.Vector
487: := Column_Vector_Rep(nd.top,nd.bottom,planes(i).all);
488: solvec : Standard_Complex_Vectors.Vector
489: (plavec'first..plavec'last+2+addmix);
490: begin
491: solvec(plavec'range) := plavec;
492: solvec(solvec'last) := Create(1.0); -- target value for t
493: solvec(solvec'last-1) := locsols(i)(locsols(i)'last-1);
494: if addmix = 1
495: then solvec(solvec'last-2) := locsols(i)(locsols(i)'last-2);
496: end if;
497: if outlog
498: then evahom := Eval(homsys,solvec);
499: put(file,"Residual of target solution :");
500: put(file,Max_Norm(evahom),3); put_line(file,".");
501: end if;
502: end;
503: end loop;
504: Clear(lochom); Clear(locsols);
505: end Quantum_Trace_Paths;
506:
507: end Pieri_Continuation;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to pieri_continuation.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