File: [local] / OpenXM_contrib / PHC / Ada / Schubert / solve_pieri_leaves.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:32 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_Complex_Vectors_io; use Standard_Complex_Vectors_io;
with Standard_Natural_Vectors;
with Standard_Complex_Vectors;
with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
function Solve_Pieri_Leaves
( file : file_type; b1,b2 : Bracket; m : Matrix ) return Matrix is
res : Matrix(m'range(1),b1'range);
n : constant natural := m'last(1);
dim : constant natural := m'last(2)-1; -- dimension of linear system
equ : Bracket(1..dim); -- equations : res(equ(i),*) = 0
sys : Matrix(1..dim,1..dim);
rhs : Standard_Complex_Vectors.Vector(1..dim);
gen : Standard_Complex_Vectors.Vector(1..n);
procedure Initialize_Solution_Plane is
-- DESCRIPTION :
-- Fills in the patterns of the resulting p-plane and sets up the
-- equations that the resulting plane has to satisfy.
-- These equations define the space C.
cnt : natural := 0;
zero_row : boolean;
begin
for i in 1..n loop
zero_row := true;
for j in b1'range loop
if ((i < b2(j)) or (i > n+1-b1(b1'last+1-j)))
then null;
else zero_row := false;
end if;
end loop;
if zero_row
then cnt := cnt+1;
equ(cnt) := i;
end if;
end loop;
end Initialize_Solution_Plane;
procedure Solve_Linear_System is
-- DESCRIPTION :
-- Solves the linear system that determines the solution plane.
-- The solution to the linear system gives the coefficients in the
-- linear combination of the columns of the given m-plane that
-- determines the generator for the line of intersection with the
-- space C and the given m-plane.
ipvt : Standard_Natural_Vectors.Vector(1..dim);
info : natural;
-- orgsys : Standard_Complex_Matrices.Matrix(1..dim,1..dim);
-- orgrhs,resi : Standard_Complex_Vectors.Vector(1..dim);
use Standard_Complex_Vectors,Standard_Complex_Matrices;
begin
for j in 1..dim loop
for i in equ'range loop
sys(i,j) := m(equ(i),j);
end loop;
end loop;
-- orgsys := sys;
for i in 1..dim loop
rhs(i) := m(equ(i),m'last(2));
end loop;
-- orgrhs := rhs;
lufac(sys,dim,ipvt,info);
lusolve(sys,dim,ipvt,rhs);
for i in 1..n loop -- linear combination of the columns of m-plane
gen(i) := -m(i,m'last(2));
for j in rhs'range loop
gen(i) := gen(i) + rhs(j)*m(i,j);
end loop;
end loop;
-- put_line(file,"The solution to the linear system : ");
-- put_line(file,rhs); new_line(file);
-- resi := orgrhs - orgsys*rhs;
-- put_line(file,"The residual vector of the linear system : ");
-- put(file,resi,2); new_line(file);
-- put_line(file,"The generator of the intersection : ");
-- put_line(file,gen); new_line(file);
end Solve_Linear_System;
procedure Determine_Solution_Plane is
-- DESCRIPTION :
-- Determines the solution plane from the generator of the intersection
-- of the space C with the given m-plane.
begin
for i in 1..n loop
for j in b1'range loop
if ((i < b2(j)) or (i > n+1-b1(b1'last+1-j)))
then res(i,j) := Create(0.0);
else res(i,j) := gen(i);
end if;
end loop;
end loop;
end Determine_Solution_Plane;
begin
Initialize_Solution_Plane;
Solve_Linear_System;
Determine_Solution_Plane;
return res;
end Solve_Pieri_Leaves;