Annotation of OpenXM_contrib/PHC/Ada/Schubert/determinantal_systems.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 2: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 3: with Standard_Natural_Vectors;
! 4: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 5: with Bracket_Monomials; use Bracket_Monomials;
! 6: with Bracket_Systems; use Bracket_Systems;
! 7: with Evaluated_Minors; use Evaluated_Minors;
! 8: with Symbolic_Minor_Equations; use Symbolic_Minor_Equations;
! 9: with Numeric_Minor_Equations; use Numeric_Minor_Equations;
! 10:
! 11: package body Determinantal_Systems is
! 12:
! 13: -- AUXILIARY TO EVALUATION OF MAXIMAL MINORS :
! 14:
! 15: function Number_of_Maximal_Minors ( nrows,ncols : natural ) return natural is
! 16:
! 17: -- DESCRIPTION :
! 18: -- Returns the number of maximal minors of a matrix with nrows and ncols.
! 19:
! 20: -- REQUIRED : nrows > ncols.
! 21:
! 22: begin
! 23: if ncols = 1
! 24: then return nrows;
! 25: else return Number_of_Maximal_Minors(nrows-1,ncols-1)*nrows/ncols;
! 26: end if;
! 27: end Number_of_Maximal_Minors;
! 28:
! 29: -- LOCALIZATION MAPS :
! 30:
! 31: function Localize ( locmap : Standard_Natural_Matrices.Matrix;
! 32: t : Term ) return Term is
! 33:
! 34: -- DESCRIPTION :
! 35: -- Applies the localization map to the term, eliminating those xij's
! 36: -- xij for which the corresponding entry in locmap is either 0 or 1.
! 37:
! 38: -- NOTE :
! 39: -- This localization assumes that t.dg(k) = 0 with k for which the
! 40: -- corresponding (i,j) with locmap(i,j) = 0.
! 41:
! 42: res : Term;
! 43: ndg : Standard_Natural_Vectors.Vector(t.dg'range);
! 44: cnt : natural := t.dg'first-1;
! 45: ind : natural := cnt;
! 46:
! 47: begin
! 48: for i in locmap'range(1) loop -- lexicographic order of variables
! 49: for j in locmap'range(2) loop
! 50: ind := ind+1;
! 51: if locmap(i,j) = 2
! 52: then cnt := cnt + 1;
! 53: ndg(cnt) := t.dg(ind);
! 54: end if;
! 55: end loop;
! 56: end loop;
! 57: for i in ind+1..t.dg'last loop -- leave the lifting !
! 58: cnt := cnt+1;
! 59: ndg(cnt) := t.dg(i);
! 60: end loop;
! 61: res.cf := t.cf;
! 62: res.dg := new Standard_Natural_Vectors.Vector'(ndg(1..cnt));
! 63: return res;
! 64: end Localize;
! 65:
! 66: function Localize ( locmap : Standard_Natural_Matrices.Matrix;
! 67: p : Poly ) return Poly is
! 68:
! 69: -- DESCRIPTION :
! 70: -- Applies the localization map to the polynomial, eliminating
! 71: -- those xij's for which locmap(i,j) is either 0 or 1.
! 72:
! 73: res : Poly := Null_Poly;
! 74:
! 75: procedure Localize_Term ( t : in Term; continue : out boolean ) is
! 76:
! 77: lt : Term := Localize(locmap,t);
! 78:
! 79: begin
! 80: Add(res,lt);
! 81: Clear(lt.dg);
! 82: continue := true;
! 83: end Localize_Term;
! 84: procedure Localize_Terms is new Visiting_Iterator(Localize_Term);
! 85:
! 86: begin
! 87: Localize_Terms(p);
! 88: return res;
! 89: end Localize;
! 90:
! 91: -- TARGET ROUTINES :
! 92:
! 93: function Standard_Coordinate_Frame
! 94: ( x : Standard_Complex_Poly_Matrices.Matrix;
! 95: plane : Standard_Complex_Matrices.Matrix )
! 96: return Standard_Natural_Matrices.Matrix is
! 97:
! 98: res : Standard_Natural_Matrices.Matrix(x'range(1),x'range(2));
! 99: tol : constant double_float := 10.0**(-10);
! 100: first : boolean;
! 101:
! 102: begin
! 103: for j in res'range(2) loop
! 104: first := true;
! 105: for i in res'range(1) loop
! 106: if x(i,j) = Null_Poly
! 107: then res(i,j) := 0;
! 108: elsif (first and (AbsVal(plane(i,j)) > tol))
! 109: then res(i,j) := 1; first := false;
! 110: else res(i,j) := 2;
! 111: end if;
! 112: end loop;
! 113: end loop;
! 114: return res;
! 115: end Standard_Coordinate_Frame;
! 116:
! 117: function Maximal_Coordinate_Frame
! 118: ( x : Standard_Complex_Poly_Matrices.Matrix;
! 119: plane : Standard_Complex_Matrices.Matrix )
! 120: return Standard_Natural_Matrices.Matrix is
! 121:
! 122: res : Standard_Natural_Matrices.Matrix(x'range(1),x'range(2));
! 123: max,tmp : double_float;
! 124: ind : natural;
! 125:
! 126: begin
! 127: for j in res'range(2) loop
! 128: for i in res'range(1) loop
! 129: if x(i,j) = Null_Poly
! 130: then res(i,j) := 0;
! 131: else res(i,j) := 2;
! 132: end if;
! 133: end loop;
! 134: max := 0.0; ind := 0;
! 135: for i in res'range(1) loop
! 136: tmp := AbsVal(plane(i,j));
! 137: if tmp > max
! 138: then max := tmp; ind := i;
! 139: end if;
! 140: end loop;
! 141: res(ind,j) := 1;
! 142: end loop;
! 143: return res;
! 144: end Maximal_Coordinate_Frame;
! 145:
! 146: function Localize ( locmap : Standard_Natural_Matrices.Matrix;
! 147: p : Poly_Sys ) return Poly_Sys is
! 148:
! 149: res : Poly_Sys(p'range);
! 150:
! 151: begin
! 152: for i in p'range loop
! 153: res(i) := Localize(locmap,p(i));
! 154: end loop;
! 155: return res;
! 156: end Localize;
! 157:
! 158: -- CONSTRUCT THE POLYNOMIAL EQUATIONS :
! 159:
! 160: function Polynomial_Equations
! 161: ( l : Standard_Complex_Matrices.Matrix;
! 162: x : Standard_Complex_Poly_Matrices.Matrix ) return Poly_Sys is
! 163:
! 164: n : constant natural := x'length(1);
! 165: p : constant natural := x'length(2);
! 166: kd : constant natural := p + l'length(2);
! 167: bm : Bracket_Monomial := Maximal_Minors(n,kd);
! 168: bs : Bracket_System(0..Number_of_Brackets(bm))
! 169: := Minor_Equations(kd,kd-p,bm);
! 170: res : Poly_Sys(bs'first+1..bs'last) := Expanded_Minors(l,x,bs);
! 171:
! 172: begin
! 173: Clear(bm); Clear(bs);
! 174: return res;
! 175: end Polynomial_Equations;
! 176:
! 177: procedure Concat ( l : in out Link_to_Poly_Sys; p : Poly_Sys ) is
! 178: begin
! 179: if l = null
! 180: then declare
! 181: newsys : Poly_Sys(p'range);
! 182: cnt : natural := p'first-1;
! 183: begin
! 184: for i in p'range loop
! 185: if p(i) /= Null_Poly
! 186: then cnt := cnt+1;
! 187: newsys(cnt) := p(i);
! 188: end if;
! 189: end loop;
! 190: l := new Poly_Sys'(newsys(p'first..cnt));
! 191: end;
! 192: else declare
! 193: newsys : Poly_Sys(l'first..l'last+p'length);
! 194: ind : natural := l'last;
! 195: begin
! 196: newsys(l'range) := l(l'range);
! 197: for i in p'range loop
! 198: if p(i) /= Null_Poly
! 199: then ind := ind+1;
! 200: newsys(ind) := p(i);
! 201: end if;
! 202: end loop;
! 203: l := new Poly_Sys'(newsys(l'first..ind));
! 204: end;
! 205: end if;
! 206: end Concat;
! 207:
! 208: function Polynomial_Equations
! 209: ( l : Standard_Complex_VecMats.VecMat;
! 210: x : Standard_Complex_Poly_Matrices.Matrix ) return Poly_Sys is
! 211:
! 212: n : constant natural := x'length(1);
! 213: p : constant natural := x'length(2);
! 214: kd : constant natural := p + l(l'first)'length(2);
! 215: bm : Bracket_Monomial := Maximal_Minors(n,kd);
! 216: bs : Bracket_System(0..Number_of_Brackets(bm))
! 217: := Minor_Equations(kd,kd-p,bm);
! 218: nb : constant natural := l'length;
! 219: res : Poly_Sys(bs'first+1..nb*bs'last);
! 220: cnt : natural := bs'first;
! 221:
! 222: begin
! 223: for i in 1..nb loop
! 224: declare
! 225: sys : constant Poly_Sys := Expanded_Minors(l(i).all,x,bs);
! 226: begin
! 227: for j in sys'range loop
! 228: cnt := cnt + 1;
! 229: res(cnt) := sys(j);
! 230: end loop;
! 231: end;
! 232: end loop;
! 233: Clear(bm); Clear(bs);
! 234: return res;
! 235: end Polynomial_Equations;
! 236:
! 237: -- EVALUATORS AND DIFFERENTIATORS :
! 238:
! 239: function Eval ( l,x : Standard_Complex_Matrices.Matrix )
! 240: return Complex_Number is
! 241:
! 242: wrk : Standard_Complex_Matrices.Matrix(x'range(1),x'range(1));
! 243:
! 244: begin
! 245: for j in l'range(2) loop
! 246: for i in l'range(1) loop
! 247: wrk(i,j) := l(i,j);
! 248: end loop;
! 249: end loop;
! 250: for j in x'range(2) loop
! 251: for i in x'range(1) loop
! 252: wrk(i,l'last(2)+j) := x(i,j);
! 253: end loop;
! 254: end loop;
! 255: return Determinant(wrk);
! 256: end Eval;
! 257:
! 258: function Eval ( l,x : Standard_Complex_Matrices.Matrix ) return Vector is
! 259:
! 260: n : constant natural := l'length(2) + x'length(2);
! 261: dim : constant natural := Number_of_Maximal_Minors(l'length(1),n);
! 262: res : Standard_Complex_Vectors.Vector(1..dim);
! 263: -- := (1..dim => Create(0.0));
! 264: cnt : natural := 0;
! 265: rws : Standard_Natural_Vectors.Vector(1..n);
! 266: wrk : Standard_Complex_Matrices.Matrix(1..n,1..n);
! 267:
! 268: procedure Choose_next_Row ( k,start : in natural ) is
! 269:
! 270: -- DESCRIPTION :
! 271: -- Chooses next k-th row within the range start..l'last(1), if k <= n.
! 272: -- If k > n, then the minor defined by the current row selection
! 273: -- is computed and added to the result.
! 274:
! 275: begin
! 276: if k > n
! 277: then for j in l'range(2) loop -- select rows
! 278: for i in 1..n loop
! 279: wrk(i,j) := l(rws(i),j);
! 280: end loop;
! 281: end loop;
! 282: for j in x'range(2) loop
! 283: for i in 1..n loop
! 284: wrk(i,l'last(2)+j) := x(rws(i),j);
! 285: end loop;
! 286: end loop;
! 287: cnt := cnt + 1;
! 288: res(cnt) := Determinant(wrk);
! 289: else for i in start..l'last(1) loop
! 290: rws(k) := i;
! 291: Choose_next_Row(k+1,i+1);
! 292: end loop;
! 293: end if;
! 294: end Choose_next_Row;
! 295:
! 296: begin
! 297: Choose_next_Row(1,1);
! 298: return res;
! 299: end Eval;
! 300:
! 301: function Minor ( l,x : Standard_Complex_Matrices.Matrix; row,col : natural )
! 302: return Complex_Number is
! 303:
! 304: -- DESCRIPTION :
! 305: -- Returns the minor [l|x] with row and column removed and with
! 306: -- the appropriate sign.
! 307:
! 308: wrk : Standard_Complex_Matrices.Matrix
! 309: (x'first(1)..x'last(1)-1,x'first(1)..x'last(1)-1);
! 310: ind : natural;
! 311:
! 312: begin
! 313: for j in l'range(2) loop
! 314: for i in l'range(1) loop
! 315: if i < row
! 316: then wrk(i,j) := l(i,j);
! 317: elsif i > row
! 318: then wrk(i-1,j) := l(i,j);
! 319: end if;
! 320: end loop;
! 321: end loop;
! 322: for j in x'range(2) loop
! 323: if j < col
! 324: then ind := l'last(2) + j;
! 325: elsif j > col
! 326: then ind := l'last(2) + j - 1;
! 327: end if;
! 328: if j /= col
! 329: then for i in x'range(1) loop
! 330: if i < row
! 331: then wrk(i,ind) := x(i,j);
! 332: elsif i > row
! 333: then wrk(i-1,ind) := x(i,j);
! 334: end if;
! 335: end loop;
! 336: end if;
! 337: end loop;
! 338: if (row + l'last(2) + col) mod 2 = 0
! 339: then return Determinant(wrk);
! 340: else return -Determinant(wrk);
! 341: end if;
! 342: end Minor;
! 343:
! 344: function Diff ( l,x : Standard_Complex_Matrices.Matrix; i : natural )
! 345: return Complex_Number is
! 346:
! 347: p : constant natural := x'length(2);
! 348: row,col : natural;
! 349:
! 350: begin
! 351: row := i/p + 1;
! 352: col := i mod p;
! 353: if col = 0
! 354: then col := p;
! 355: row := row - 1;
! 356: end if;
! 357: return Minor(l,x,row,col);
! 358: end Diff;
! 359:
! 360: function Diff ( l,x : Standard_Complex_Matrices.Matrix;
! 361: locmap : Standard_Natural_Matrices.Matrix; i : natural )
! 362: return Complex_Number is
! 363:
! 364: row,col : natural;
! 365: cnt : natural := 0;
! 366:
! 367: begin
! 368: for k1 in locmap'range(1) loop
! 369: for k2 in locmap'range(2) loop
! 370: if locmap(k1,k2) = 2
! 371: then cnt := cnt+1;
! 372: if cnt = i
! 373: then row := k1;
! 374: col := k2;
! 375: end if;
! 376: end if;
! 377: exit when (cnt = i);
! 378: end loop;
! 379: exit when (cnt = i);
! 380: end loop;
! 381: return Minor(l,x,row,col);
! 382: end Diff;
! 383:
! 384: function Eval ( l : Standard_Complex_VecMats.VecMat;
! 385: x : Standard_Complex_Matrices.Matrix )
! 386: return Standard_Complex_Vectors.Vector is
! 387:
! 388: res : Standard_Complex_Vectors.Vector(l'range);
! 389:
! 390: begin
! 391: for i in l'range loop
! 392: res(i) := Eval(l(i).all,x);
! 393: end loop;
! 394: return res;
! 395: end Eval;
! 396:
! 397: function Diff ( l : Standard_Complex_VecMats.VecMat;
! 398: x : Standard_Complex_Matrices.Matrix )
! 399: return Standard_Complex_Matrices.Matrix is
! 400:
! 401: res : Standard_Complex_Matrices.Matrix(l'range,1..x'last(1)*x'last(2));
! 402:
! 403: begin
! 404: for i in res'range(1) loop
! 405: for j in res'range(2) loop
! 406: res(i,j) := Diff(l(i).all,x,j);
! 407: end loop;
! 408: end loop;
! 409: return res;
! 410: end Diff;
! 411:
! 412: function Old_Diff ( l : Standard_Complex_VecMats.VecMat;
! 413: x : Standard_Complex_Matrices.Matrix; nvars : natural;
! 414: locmap : Standard_Natural_Matrices.Matrix )
! 415: return Standard_Complex_Matrices.Matrix is
! 416:
! 417: res : Standard_Complex_Matrices.Matrix(l'range,1..nvars);
! 418:
! 419: begin
! 420: for i in res'range(1) loop
! 421: for j in res'range(2) loop
! 422: res(i,j) := Diff(l(i).all,x,locmap,j);
! 423: end loop;
! 424: end loop;
! 425: return res;
! 426: end Old_Diff;
! 427:
! 428: function Diff ( l : Standard_Complex_VecMats.VecMat;
! 429: x : Standard_Complex_Matrices.Matrix; nvars : natural;
! 430: locmap : Standard_Natural_Matrices.Matrix )
! 431: return Standard_Complex_Matrices.Matrix is
! 432:
! 433: -- NOTE :
! 434: -- This Diff is organized according to the localization map,
! 435: -- to avoid multiple searches.
! 436:
! 437: res : Standard_Complex_Matrices.Matrix(l'range,1..nvars);
! 438: ind : natural := 0;
! 439:
! 440: begin
! 441: for i in locmap'range loop
! 442: for j in locmap'range loop
! 443: if locmap(i,j) = 2
! 444: then ind := ind+1;
! 445: for k in l'range loop
! 446: res(k,ind) := Minor(l(k).all,x,i,j);
! 447: end loop;
! 448: end if;
! 449: end loop;
! 450: end loop;
! 451: return res;
! 452: end Diff;
! 453:
! 454: end Determinantal_Systems;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to determinantal_systems.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