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