File: [local] / OpenXM_contrib / PHC / Ada / Schubert / specialization_of_planes.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:33 2000 UTC (23 years, 8 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
with Standard_Complex_Numbers; use Standard_Complex_Numbers;
with Standard_Random_Numbers; use Standard_Random_Numbers;
with Standard_Natural_Vectors;
with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
package body Specialization_of_Planes is
function Random_Upper_Triangular
( n : natural ) return Standard_Complex_Matrices.Matrix is
res : Standard_Complex_Matrices.Matrix(1..n,1..n);
begin
for j in 1..n loop -- assign values to jth column
for i in 1..n-j loop
res(i,j) := Random1; -- randoms above anti-diagonal
end loop;
res(n-j+1,j) := Create(1.0); -- 1 = anti-diagonal element
for i in n-j+2..n loop
res(i,j) := Create(0.0); -- zeros under anti-diagonal
end loop;
end loop;
return res;
end Random_Upper_Triangular;
function Random_Lower_Triangular
( n : natural ) return Standard_Complex_Matrices.Matrix is
res : Standard_Complex_Matrices.Matrix(1..n,1..n);
begin
for j in 1..n loop -- assign values to jth column
for i in 1..(j-1) loop
res(i,j) := Create(0.0); -- zeros above diagonal
end loop;
res(j,j) := Create(1.0); -- 1 = diagonal element
for i in (j+1)..n loop
res(i,j) := Random1; -- randoms under diagonal
end loop;
end loop;
return res;
end Random_Lower_Triangular;
function U_Matrix ( F : Standard_Complex_Matrices.Matrix; b : Bracket )
return Standard_Complex_Matrices.Matrix is
m : constant natural := F'length(1) - b'length;
res : Standard_Complex_Matrices.Matrix(F'range(1),1..m);
rvf : Standard_Complex_Matrices.Matrix(F'range(1),F'range(2)) := F;
rng : constant natural := F'length(2) - b(b'last);
tmp : Complex_Number;
ind : natural := 1;
cnt : natural := 0;
begin
for j in 1..(rng/2) loop -- reverse last columns
for i in F'range(1) loop
tmp := rvf(i,F'last(2)-j+1);
rvf(i,F'last(2)-j+1) := rvf(i,F'last(2)-rng+j);
rvf(i,F'last(2)-rng+j) := tmp;
end loop;
end loop;
for j in F'range(2) loop -- remove columns indexed by b
if ((ind <= b'last) and then (j = b(ind)))
then ind := ind+1;
else cnt := cnt+1;
for i in F'range(1) loop
res(i,cnt) := rvf(i,j);
end loop;
end if;
end loop;
return res;
end U_Matrix;
function Special_Plane ( m : natural; b : Bracket )
return Standard_Complex_Matrices.Matrix is
p : constant natural := b'length;
n : constant natural := m+p;
res : Standard_Complex_Matrices.Matrix(1..n,1..m);
row,col : natural;
begin
row := 1; col := 0;
for i in res'range(1) loop
for j in res'range(2) loop
res(i,j) := Create(0.0);
end loop;
if ((row <= p) and then (b(row) = i))
then row := row + 1;
else col := col + 1;
res(i,col) := Create(1.0);
end if;
end loop;
return res;
end Special_Plane;
function Special_Bottom_Plane ( m : natural; b : Bracket )
return Standard_Complex_Matrices.Matrix is
p : constant natural := b'length;
n : constant natural := m+p;
res : Standard_Complex_Matrices.Matrix(1..n,1..m);
row,col : natural;
begin
row := 1; col := 0;
for i in res'range(1) loop
if ((row <= p) and then (b(row) = i))
then row := row + 1;
else col := col + 1;
for k in 1..i-1 loop -- randoms above the diagonal
res(k,col) := Random1;
end loop;
res(i,col) := Create(1.0);
for k in i+1..n loop -- zeros below the diagonal
res(k,col) := Create(0.0);
end loop;
end if;
end loop;
return res;
end Special_Bottom_Plane;
function Special_Top_Plane ( m : natural; b : Bracket )
return Standard_Complex_Matrices.Matrix is
p : constant natural := b'length;
n : constant natural := m+p;
res : Standard_Complex_Matrices.Matrix(1..n,1..m);
row,col : natural;
begin
row := 1; col := 0;
for i in res'range(1) loop
if ((row <= p) and then (b(row) = i))
then row := row + 1;
else col := col + 1;
for k in 1..i-1 loop -- zeros above the diagonal
res(k,col) := Create(0.0);
end loop;
res(i,col) := Create(1.0);
for k in i+1..n loop -- randoms below the diagonal
res(k,col) := Random1;
end loop;
end if;
end loop;
return res;
end Special_Top_Plane;
function Special_Plane
( n,m,k : natural; b : Bracket;
special : in Standard_Complex_Matrices.Matrix )
return Standard_Complex_Matrices.Matrix is
res : Standard_Complex_Matrices.Matrix(1..n,1..m+1-k);
ran : Complex_Number;
begin
for i in res'range(1) loop -- initialize
for j in res'range(2) loop
res(i,j) := Create(0.0);
end loop;
end loop;
for j in res'range(2) loop -- build j-th column
for k in b'range loop
ran := Random1; -- random for column k of mat
for i in special'range(1) loop
res(i,j) := res(i,j) + ran*special(i,b(k));
end loop;
end loop;
end loop;
return res;
end Special_Plane;
function Special_Bottom_Plane ( n,m,k : natural; b : Bracket )
return Standard_Complex_Matrices.Matrix is
mat : Standard_Complex_Matrices.Matrix(1..n,b'range);
begin
for j in b'range loop -- j-th column of matrix
for i in 1..b(j) loop
mat(i,j) := Random1; -- random numbers above row b(j)
end loop;
for i in b(j)+1..n loop
mat(i,j) := Create(0.0); -- zeros below row b(j)
end loop;
end loop;
return Special_Plane(n,m,k,b,mat);
end Special_Bottom_Plane;
function Special_Top_Plane ( n,m,k : natural; b : Bracket )
return Standard_Complex_Matrices.Matrix is
mat : Standard_Complex_Matrices.Matrix(1..n,b'range);
begin
for j in b'range loop -- j-th column of matrix
for i in 1..b(j)-1 loop
mat(i,j) := Create(0.0); -- zeros below row b(j)
end loop;
for i in b(j)..n loop
mat(i,j) := Random1; -- random numbers below row b(j)
end loop;
end loop;
return Special_Plane(n,m,k,b,mat);
end Special_Top_Plane;
function Homotopy ( dim : natural; start,target : Complex_Number )
return Poly is
-- DESCRIPTION :
-- Returns the polynomial start*(1-t) + target*t, where t is the
-- is the last variable of index dim.
-- This procedure is an auxiliary to building the moving U-matrices.
res : Poly;
t : Term;
tdg : Degrees := new Standard_Natural_Vectors.Vector'(1..dim => 0);
begin
t.cf := start;
t.dg := tdg;
res := Create(t); -- res = start
tdg(tdg'last) := 1; -- introduce t
t.dg := tdg;
Sub(res,t); -- res = (1-t)*start
t.cf := target;
Add(res,t); -- res = (1-t)*start + t*target
Clear(tdg);
return res;
end Homotopy;
function Constant_Poly ( dim : natural; c : Complex_Number ) return Poly is
-- DESCRIPTION :
-- Returns the constant c represented as polynomial with as many
-- variables as the given dimension.
res : Poly;
t : Term;
tdg : Degrees := new Standard_Natural_Vectors.Vector'(1..dim => 0);
begin
t.cf := c;
t.dg := tdg;
res := Create(t);
return res;
end Constant_Poly;
function Moving_U_Matrix
( n : natural; U,L : Standard_Complex_Matrices.Matrix )
return Standard_Complex_Poly_Matrices.Matrix is
res : Standard_Complex_Poly_Matrices.Matrix(L'range(1),L'range(2));
t : Term;
begin
for i in res'range(1) loop
for j in res'range(2) loop
res(i,j) := Homotopy(n,U(i,j),L(i,j));
end loop;
end loop;
return res;
end Moving_U_Matrix;
function Moving_U_Matrix
( U : Standard_Complex_Matrices.Matrix;
i,r : natural; b : bracket )
return Standard_Complex_Poly_Matrices.Matrix is
p : constant natural := b'last;
m : constant natural := U'length(1) - p;
dim : constant natural := (m+p)*p+1;
res : Standard_Complex_Poly_Matrices.Matrix(U'range(1),1..m+1-r);
begin
for j in res'range(2) loop
if j+i-1 < b(b'last) - b'last
then for k in res'range(1) loop
res(k,j) := Homotopy(dim,U(k,j+i),U(k,j+i-1));
end loop;
elsif j+i-1 = b(b'last) - b'last
then for k in res'range(1) loop
res(k,j) := Homotopy(dim,U(k,m+1+i-r),U(k,j+i-1));
end loop;
else for k in res'range(1) loop
res(k,j) := Constant_Poly(dim,U(k,j+i-1));
end loop;
end if;
end loop;
return res;
end Moving_U_Matrix;
function Slice
( M : Standard_Complex_Poly_Matrices.Matrix;
ind : natural ) return Standard_Complex_Poly_Matrices.Matrix is
-- DESCRIPTION :
-- Returns the columns of M up to the given index.
res : Standard_Complex_Poly_Matrices.Matrix(M'range(1),M'first(2)..ind);
begin
for j in res'range(2) loop
for i in res'range(1) loop
Copy(M(i,j),res(i,j));
end loop;
end loop;
return res;
end Slice;
function Lower_Section
( M : Standard_Complex_Poly_Matrices.Matrix;
row : natural ) return Standard_Complex_Poly_Matrices.Matrix is
res : Standard_Complex_Poly_Matrices.Matrix(M'range(1),M'range(2));
cnt : natural := M'first(2)-1;
add : boolean;
begin
for j in M'range(2) loop
for i in (row+1)..M'last(1) loop
add := (M(i,j) = Null_Poly);
exit when not add;
end loop;
if add
then cnt := cnt+1;
for i in M'range(1) loop
res(i,cnt) := M(i,j);
end loop;
end if;
end loop;
return Slice(res,cnt);
end Lower_Section;
function Upper_Section
( M : Standard_Complex_Poly_Matrices.Matrix;
row : natural ) return Standard_Complex_Poly_Matrices.Matrix is
res : Standard_Complex_Poly_Matrices.Matrix(M'range(1),M'range(2));
cnt : natural := M'first(2)-1;
add : boolean;
begin
for j in M'range(2) loop
for i in M'first(1)..(row-1) loop
add := (M(i,j) = Null_Poly);
exit when not add;
end loop;
if add
then cnt := cnt+1;
for i in M'range(1) loop
res(i,cnt) := M(i,j);
end loop;
end if;
end loop;
return Slice(res,cnt);
end Upper_Section;
end Specialization_of_Planes;