Annotation of OpenXM_contrib/PHC/Ada/Schubert/specialization_of_planes.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 2: with Standard_Random_Numbers; use Standard_Random_Numbers;
! 3: with Standard_Natural_Vectors;
! 4: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 5:
! 6: package body Specialization_of_Planes is
! 7:
! 8: function Random_Upper_Triangular
! 9: ( n : natural ) return Standard_Complex_Matrices.Matrix is
! 10:
! 11: res : Standard_Complex_Matrices.Matrix(1..n,1..n);
! 12:
! 13: begin
! 14: for j in 1..n loop -- assign values to jth column
! 15: for i in 1..n-j loop
! 16: res(i,j) := Random1; -- randoms above anti-diagonal
! 17: end loop;
! 18: res(n-j+1,j) := Create(1.0); -- 1 = anti-diagonal element
! 19: for i in n-j+2..n loop
! 20: res(i,j) := Create(0.0); -- zeros under anti-diagonal
! 21: end loop;
! 22: end loop;
! 23: return res;
! 24: end Random_Upper_Triangular;
! 25:
! 26: function Random_Lower_Triangular
! 27: ( n : natural ) return Standard_Complex_Matrices.Matrix is
! 28:
! 29: res : Standard_Complex_Matrices.Matrix(1..n,1..n);
! 30:
! 31: begin
! 32: for j in 1..n loop -- assign values to jth column
! 33: for i in 1..(j-1) loop
! 34: res(i,j) := Create(0.0); -- zeros above diagonal
! 35: end loop;
! 36: res(j,j) := Create(1.0); -- 1 = diagonal element
! 37: for i in (j+1)..n loop
! 38: res(i,j) := Random1; -- randoms under diagonal
! 39: end loop;
! 40: end loop;
! 41: return res;
! 42: end Random_Lower_Triangular;
! 43:
! 44: function U_Matrix ( F : Standard_Complex_Matrices.Matrix; b : Bracket )
! 45: return Standard_Complex_Matrices.Matrix is
! 46:
! 47: m : constant natural := F'length(1) - b'length;
! 48: res : Standard_Complex_Matrices.Matrix(F'range(1),1..m);
! 49: rvf : Standard_Complex_Matrices.Matrix(F'range(1),F'range(2)) := F;
! 50: rng : constant natural := F'length(2) - b(b'last);
! 51: tmp : Complex_Number;
! 52: ind : natural := 1;
! 53: cnt : natural := 0;
! 54:
! 55: begin
! 56: for j in 1..(rng/2) loop -- reverse last columns
! 57: for i in F'range(1) loop
! 58: tmp := rvf(i,F'last(2)-j+1);
! 59: rvf(i,F'last(2)-j+1) := rvf(i,F'last(2)-rng+j);
! 60: rvf(i,F'last(2)-rng+j) := tmp;
! 61: end loop;
! 62: end loop;
! 63: for j in F'range(2) loop -- remove columns indexed by b
! 64: if ((ind <= b'last) and then (j = b(ind)))
! 65: then ind := ind+1;
! 66: else cnt := cnt+1;
! 67: for i in F'range(1) loop
! 68: res(i,cnt) := rvf(i,j);
! 69: end loop;
! 70: end if;
! 71: end loop;
! 72: return res;
! 73: end U_Matrix;
! 74:
! 75: function Special_Plane ( m : natural; b : Bracket )
! 76: return Standard_Complex_Matrices.Matrix is
! 77:
! 78: p : constant natural := b'length;
! 79: n : constant natural := m+p;
! 80: res : Standard_Complex_Matrices.Matrix(1..n,1..m);
! 81: row,col : natural;
! 82:
! 83: begin
! 84: row := 1; col := 0;
! 85: for i in res'range(1) loop
! 86: for j in res'range(2) loop
! 87: res(i,j) := Create(0.0);
! 88: end loop;
! 89: if ((row <= p) and then (b(row) = i))
! 90: then row := row + 1;
! 91: else col := col + 1;
! 92: res(i,col) := Create(1.0);
! 93: end if;
! 94: end loop;
! 95: return res;
! 96: end Special_Plane;
! 97:
! 98: function Special_Bottom_Plane ( m : natural; b : Bracket )
! 99: return Standard_Complex_Matrices.Matrix is
! 100:
! 101: p : constant natural := b'length;
! 102: n : constant natural := m+p;
! 103: res : Standard_Complex_Matrices.Matrix(1..n,1..m);
! 104: row,col : natural;
! 105:
! 106: begin
! 107: row := 1; col := 0;
! 108: for i in res'range(1) loop
! 109: if ((row <= p) and then (b(row) = i))
! 110: then row := row + 1;
! 111: else col := col + 1;
! 112: for k in 1..i-1 loop -- randoms above the diagonal
! 113: res(k,col) := Random1;
! 114: end loop;
! 115: res(i,col) := Create(1.0);
! 116: for k in i+1..n loop -- zeros below the diagonal
! 117: res(k,col) := Create(0.0);
! 118: end loop;
! 119: end if;
! 120: end loop;
! 121: return res;
! 122: end Special_Bottom_Plane;
! 123:
! 124: function Special_Top_Plane ( m : natural; b : Bracket )
! 125: return Standard_Complex_Matrices.Matrix is
! 126:
! 127: p : constant natural := b'length;
! 128: n : constant natural := m+p;
! 129: res : Standard_Complex_Matrices.Matrix(1..n,1..m);
! 130: row,col : natural;
! 131:
! 132: begin
! 133: row := 1; col := 0;
! 134: for i in res'range(1) loop
! 135: if ((row <= p) and then (b(row) = i))
! 136: then row := row + 1;
! 137: else col := col + 1;
! 138: for k in 1..i-1 loop -- zeros above the diagonal
! 139: res(k,col) := Create(0.0);
! 140: end loop;
! 141: res(i,col) := Create(1.0);
! 142: for k in i+1..n loop -- randoms below the diagonal
! 143: res(k,col) := Random1;
! 144: end loop;
! 145: end if;
! 146: end loop;
! 147: return res;
! 148: end Special_Top_Plane;
! 149:
! 150: function Special_Plane
! 151: ( n,m,k : natural; b : Bracket;
! 152: special : in Standard_Complex_Matrices.Matrix )
! 153: return Standard_Complex_Matrices.Matrix is
! 154:
! 155: res : Standard_Complex_Matrices.Matrix(1..n,1..m+1-k);
! 156: ran : Complex_Number;
! 157:
! 158: begin
! 159: for i in res'range(1) loop -- initialize
! 160: for j in res'range(2) loop
! 161: res(i,j) := Create(0.0);
! 162: end loop;
! 163: end loop;
! 164: for j in res'range(2) loop -- build j-th column
! 165: for k in b'range loop
! 166: ran := Random1; -- random for column k of mat
! 167: for i in special'range(1) loop
! 168: res(i,j) := res(i,j) + ran*special(i,b(k));
! 169: end loop;
! 170: end loop;
! 171: end loop;
! 172: return res;
! 173: end Special_Plane;
! 174:
! 175: function Special_Bottom_Plane ( n,m,k : natural; b : Bracket )
! 176: return Standard_Complex_Matrices.Matrix is
! 177:
! 178: mat : Standard_Complex_Matrices.Matrix(1..n,b'range);
! 179:
! 180: begin
! 181: for j in b'range loop -- j-th column of matrix
! 182: for i in 1..b(j) loop
! 183: mat(i,j) := Random1; -- random numbers above row b(j)
! 184: end loop;
! 185: for i in b(j)+1..n loop
! 186: mat(i,j) := Create(0.0); -- zeros below row b(j)
! 187: end loop;
! 188: end loop;
! 189: return Special_Plane(n,m,k,b,mat);
! 190: end Special_Bottom_Plane;
! 191:
! 192: function Special_Top_Plane ( n,m,k : natural; b : Bracket )
! 193: return Standard_Complex_Matrices.Matrix is
! 194:
! 195: mat : Standard_Complex_Matrices.Matrix(1..n,b'range);
! 196:
! 197: begin
! 198: for j in b'range loop -- j-th column of matrix
! 199: for i in 1..b(j)-1 loop
! 200: mat(i,j) := Create(0.0); -- zeros below row b(j)
! 201: end loop;
! 202: for i in b(j)..n loop
! 203: mat(i,j) := Random1; -- random numbers below row b(j)
! 204: end loop;
! 205: end loop;
! 206: return Special_Plane(n,m,k,b,mat);
! 207: end Special_Top_Plane;
! 208:
! 209: function Homotopy ( dim : natural; start,target : Complex_Number )
! 210: return Poly is
! 211:
! 212: -- DESCRIPTION :
! 213: -- Returns the polynomial start*(1-t) + target*t, where t is the
! 214: -- is the last variable of index dim.
! 215: -- This procedure is an auxiliary to building the moving U-matrices.
! 216:
! 217: res : Poly;
! 218: t : Term;
! 219: tdg : Degrees := new Standard_Natural_Vectors.Vector'(1..dim => 0);
! 220:
! 221: begin
! 222: t.cf := start;
! 223: t.dg := tdg;
! 224: res := Create(t); -- res = start
! 225: tdg(tdg'last) := 1; -- introduce t
! 226: t.dg := tdg;
! 227: Sub(res,t); -- res = (1-t)*start
! 228: t.cf := target;
! 229: Add(res,t); -- res = (1-t)*start + t*target
! 230: Clear(tdg);
! 231: return res;
! 232: end Homotopy;
! 233:
! 234: function Constant_Poly ( dim : natural; c : Complex_Number ) return Poly is
! 235:
! 236: -- DESCRIPTION :
! 237: -- Returns the constant c represented as polynomial with as many
! 238: -- variables as the given dimension.
! 239:
! 240: res : Poly;
! 241: t : Term;
! 242: tdg : Degrees := new Standard_Natural_Vectors.Vector'(1..dim => 0);
! 243:
! 244: begin
! 245: t.cf := c;
! 246: t.dg := tdg;
! 247: res := Create(t);
! 248: return res;
! 249: end Constant_Poly;
! 250:
! 251: function Moving_U_Matrix
! 252: ( n : natural; U,L : Standard_Complex_Matrices.Matrix )
! 253: return Standard_Complex_Poly_Matrices.Matrix is
! 254:
! 255: res : Standard_Complex_Poly_Matrices.Matrix(L'range(1),L'range(2));
! 256: t : Term;
! 257:
! 258: begin
! 259: for i in res'range(1) loop
! 260: for j in res'range(2) loop
! 261: res(i,j) := Homotopy(n,U(i,j),L(i,j));
! 262: end loop;
! 263: end loop;
! 264: return res;
! 265: end Moving_U_Matrix;
! 266:
! 267: function Moving_U_Matrix
! 268: ( U : Standard_Complex_Matrices.Matrix;
! 269: i,r : natural; b : bracket )
! 270: return Standard_Complex_Poly_Matrices.Matrix is
! 271:
! 272: p : constant natural := b'last;
! 273: m : constant natural := U'length(1) - p;
! 274: dim : constant natural := (m+p)*p+1;
! 275: res : Standard_Complex_Poly_Matrices.Matrix(U'range(1),1..m+1-r);
! 276:
! 277: begin
! 278: for j in res'range(2) loop
! 279: if j+i-1 < b(b'last) - b'last
! 280: then for k in res'range(1) loop
! 281: res(k,j) := Homotopy(dim,U(k,j+i),U(k,j+i-1));
! 282: end loop;
! 283: elsif j+i-1 = b(b'last) - b'last
! 284: then for k in res'range(1) loop
! 285: res(k,j) := Homotopy(dim,U(k,m+1+i-r),U(k,j+i-1));
! 286: end loop;
! 287: else for k in res'range(1) loop
! 288: res(k,j) := Constant_Poly(dim,U(k,j+i-1));
! 289: end loop;
! 290: end if;
! 291: end loop;
! 292: return res;
! 293: end Moving_U_Matrix;
! 294:
! 295: function Slice
! 296: ( M : Standard_Complex_Poly_Matrices.Matrix;
! 297: ind : natural ) return Standard_Complex_Poly_Matrices.Matrix is
! 298:
! 299: -- DESCRIPTION :
! 300: -- Returns the columns of M up to the given index.
! 301:
! 302: res : Standard_Complex_Poly_Matrices.Matrix(M'range(1),M'first(2)..ind);
! 303:
! 304: begin
! 305: for j in res'range(2) loop
! 306: for i in res'range(1) loop
! 307: Copy(M(i,j),res(i,j));
! 308: end loop;
! 309: end loop;
! 310: return res;
! 311: end Slice;
! 312:
! 313: function Lower_Section
! 314: ( M : Standard_Complex_Poly_Matrices.Matrix;
! 315: row : natural ) return Standard_Complex_Poly_Matrices.Matrix is
! 316:
! 317: res : Standard_Complex_Poly_Matrices.Matrix(M'range(1),M'range(2));
! 318: cnt : natural := M'first(2)-1;
! 319: add : boolean;
! 320:
! 321: begin
! 322: for j in M'range(2) loop
! 323: for i in (row+1)..M'last(1) loop
! 324: add := (M(i,j) = Null_Poly);
! 325: exit when not add;
! 326: end loop;
! 327: if add
! 328: then cnt := cnt+1;
! 329: for i in M'range(1) loop
! 330: res(i,cnt) := M(i,j);
! 331: end loop;
! 332: end if;
! 333: end loop;
! 334: return Slice(res,cnt);
! 335: end Lower_Section;
! 336:
! 337: function Upper_Section
! 338: ( M : Standard_Complex_Poly_Matrices.Matrix;
! 339: row : natural ) return Standard_Complex_Poly_Matrices.Matrix is
! 340:
! 341: res : Standard_Complex_Poly_Matrices.Matrix(M'range(1),M'range(2));
! 342: cnt : natural := M'first(2)-1;
! 343: add : boolean;
! 344:
! 345: begin
! 346: for j in M'range(2) loop
! 347: for i in M'first(1)..(row-1) loop
! 348: add := (M(i,j) = Null_Poly);
! 349: exit when not add;
! 350: end loop;
! 351: if add
! 352: then cnt := cnt+1;
! 353: for i in M'range(1) loop
! 354: res(i,cnt) := M(i,j);
! 355: end loop;
! 356: end if;
! 357: end loop;
! 358: return Slice(res,cnt);
! 359: end Upper_Section;
! 360:
! 361: end Specialization_of_Planes;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to specialization_of_planes.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