File: [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift / binomial_system_solvers.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:28 2000 UTC (23 years, 10 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_Numbers_Polar; use Standard_Complex_Numbers_Polar;
with Standard_Integer_Vectors;
with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
with Transforming_Solutions; use Transforming_Solutions;
package body Binomial_System_Solvers is
-- AN INTERMEDIATE ROUTINE : --
procedure Factorize ( v : in out VecVec; n : in natural;
t : in out Transfo ) is
-- DESCRIPTION :
-- This routines factorizes the binomial system defined by the
-- degrees in the Vector v;
-- ON ENTRY :
-- v defines the degrees of binomial system
-- n the number of unknowns to be eliminated
-- ON RETURN :
-- v the factorized binomial system
-- t the transformations used to factorize
tt : Transfo;
pivot : integer;
tmp : Standard_Integer_Vectors.Link_to_Vector;
begin
t := Id(v'last);
for i in 1..n loop
pivot := 0; -- search for pivot
for j in i..v'last loop
if v(j)(i) /= 0
then pivot := j;
exit;
end if;
end loop;
if pivot /= 0 -- interchange if necessary
then if pivot /= i
then tmp := v(i);
v(i) := v(pivot);
v(pivot) := tmp;
end if;
tt := Rotate(v(i),i);
Apply(tt,v(i));
for j in (i+1)..v'last loop
Apply(tt,v(j));
v(j).all(i) := 0;
end loop;
Mult1(t,tt);
Clear(tt);
end if;
end loop;
end Factorize;
-- AUXILIARY ROUTINES :
procedure Update ( sols1,sols2 : in out Solution_List ) is
tmp : Solution_List := sols2;
ls : Link_to_Solution;
begin
while not Is_Null(tmp) loop
ls := Head_Of(tmp);
declare
nls : Link_to_Solution := new Solution'(ls.all);
begin
Construct(nls,sols1);
end;
tmp := Tail_Of(tmp);
end loop;
Clear(sols2);
end Update;
procedure Solve ( v : in VecVec; cv : in Vector;
i,n : in natural; sols : in out Solution_List;
acc : in out Vector ) is
t : Transfo;
pivot,d : integer;
a : Complex_Number;
c : Vector(cv'range) := cv;
rv : Vector(cv'range);
tmp : Standard_Integer_Vectors.Link_to_Vector;
nv : VecVec(v'range);
temp : array(v'range) of integer;
newsols : Solution_List;
begin
for j in i..v'last loop
nv(j) := new Standard_Integer_Vectors.Vector'(v(j).all);
end loop;
pivot := 0; -- search for pivot
for j in i..v'last loop
if nv(j)(i) /= 0
then pivot := j;
exit;
end if;
end loop;
if pivot /= 0 -- interchange if necessary
then if pivot /= i
then tmp := nv(i);
nv(i) := nv(pivot);
nv(pivot) := tmp;
a := c(i);
c(i) := c(pivot);
c(pivot) := a;
end if;
t := Rotate(nv(i),i);
Apply(t,nv(i));
for j in (i+1)..nv'last loop
Apply(t,nv(j));
temp(j) := nv(j)(i);
nv(j)(i) := 0;
end loop;
d := nv(i)(i); -- compute the solutions
if d < 0
then d := -d;
rv(i) := Create(1.0)/c(i);
else rv(i) := c(i);
end if;
for j in 1..d loop
a := Root(rv(i),d,j);
acc(i) := a;
for k in (i+1)..nv'last loop
rv(k) := c(k)*a**(-temp(k));
end loop;
if i < n
then Solve(nv,rv,i+1,n,newsols,acc);
else declare -- i = n
ls : Link_to_Solution;
s : Solution(n);
begin
s.t := Create(0.0);
s.m := 1;
s.v := acc;
ls := new Solution'(s);
Construct(ls,newsols);
end;
end if;
Transform(t,newsols);
Update(sols,newsols);
end loop;
Clear(t);
end if;
for j in i..nv'last loop
Standard_Integer_Vectors.Clear(nv(j));
end loop;
end Solve;
-- THE SOLVER :
procedure Solve ( v : in VecVec; cv : in Vector;
n : in natural; sols : in out Solution_List ) is
wv : VecVec(v'range); -- workspace
acc : Vector(cv'range); -- accumulator
begin
for i in v'range loop
wv(i) := new Standard_Integer_Vectors.Vector'(v(i).all);
end loop;
Solve(wv,cv,v'first,v'last,sols,acc);
Clear(wv);
end Solve;
-- COMPUTING THE RESIDUALS :
procedure Residuals ( v : in VecVec; cv : in Vector;
n : in natural; sols : in Solution_List;
res : out Vector ) is
nb : natural;
x : Vector(cv'range);
tmp : Solution_List;
function Eval ( v : VecVec; cv,x : Vector ) return Vector is
-- DESCRIPTION :
-- Computes the value of the vector x in the system defined by
-- v and cv.
y : Standard_Complex_Vectors.Vector(x'range);
begin
for i in x'range loop
y(i) := Create(1.0);
for j in v(i)'range loop
y(i) := y(i)*x(j)**v(i)(j);
end loop;
y(i) := y(i) - cv(i);
end loop;
return y;
end Eval;
begin
tmp := sols;
nb := 1;
while not Is_Null(tmp) loop
x := Head_Of(tmp).v;
res(nb) := Create(Max_Norm(Eval(v,cv,x)));
tmp := Tail_Of(tmp);
nb := nb + 1;
end loop;
end Residuals;
end Binomial_System_Solvers;