Annotation of OpenXM_contrib/PHC/Ada/Schubert/pieri_continuation.adb, Revision 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