File: [local] / OpenXM_contrib / PHC / Ada / Schubert / pieri_continuation.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_Floating_Numbers; use Standard_Floating_Numbers;
with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
with Standard_Complex_Numbers; use Standard_Complex_Numbers;
with Standard_Random_Numbers; use Standard_Random_Numbers;
with Standard_Natural_Vectors;
with Standard_Natural_Vectors_io; use Standard_Natural_Vectors_io;
with Standard_Complex_Vectors;
with Standard_Complex_Vectors_io; use Standard_Complex_Vectors_io;
with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
with Standard_Complex_VecVecs; use Standard_Complex_VecVecs;
with Standard_Natural_Matrices_io; use Standard_Natural_Matrices_io;
with Standard_Complex_Matrices_io; use Standard_Complex_Matrices_io;
with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
with Standard_Complex_Solutions; use Standard_Complex_Solutions;
with Homotopy;
with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
with Standard_Root_Refiners; use Standard_Root_Refiners;
with Determinantal_Systems; use Determinantal_Systems;
with Plane_Representations; use Plane_Representations;
with Curves_into_Grassmannian; use Curves_into_Grassmannian;
with Curves_into_Grassmannian_io; use Curves_into_Grassmannian_io;
package body Pieri_Continuation is
-- AUXILIARY FORMAT CONVERTORS :
function Create ( v : Standard_Complex_Vectors.Vector; conpar : natural )
return Solution is
-- DESCRIPTION :
-- Returns a solution that contains a solution with vector v in it.
-- The v(conpar) is the continuation parameter.
sol : Solution(v'last-1);
begin
sol.t := v(conpar);
sol.m := 1;
sol.v(v'first..conpar-1) := v(v'first..conpar-1);
sol.v(conpar..sol.v'last) := v(conpar+1..v'last);
sol.err := 0.0;
sol.rco := 0.0;
sol.res := 0.0;
return sol;
end Create;
function Create ( v : Standard_Complex_Vectors.Vector )
return Solution_List is
-- DESCRIPTION :
-- Returns a solution list that contains a solution with vector v in it.
-- NOTE : only needed in the original Pieri homotopy algorithm.
res : Solution_List;
begin
Add(res,Create(v,v'last));
return res;
end Create;
function Create ( v : Standard_Complex_VecVecs.VecVec; conpar : natural )
return Solution_List is
-- DESCRIPTION :
-- Returns a solution list that contains the vectors in v.
-- REQUIRED : v is not empty.
-- The value for the continuation parameter is in v(conpar).
res,res_last : Solution_List;
begin
for i in v'range loop
Append(res,res_last,Create(v(i).all,conpar));
end loop;
return res;
end Create;
function Convert ( sols : Solution_List )
return Standard_Complex_Vectors.Vector is
-- DESCRIPTION :
-- Converts the solutions in sols to the vector representation.
-- NOTE : only needed in the original Pieri homotopy algorithm.
ls : Link_to_Solution := Head_Of(sols);
res : Standard_Complex_Vectors.Vector(1..ls.n+1);
begin
res(ls.v'range) := ls.v;
res(res'last) := ls.t;
return res;
end Convert;
function Convert ( sol : Solution; conpar : natural )
return Standard_Complex_Vectors.Vector is
-- DESCRIPTION :
-- Converts the solution to the vector representation, where the
-- continuation parameter is put into the entry indexed by conpar.
res : Standard_Complex_Vectors.Vector(1..sol.n+1);
begin
res(sol.v'first..conpar-1) := sol.v(sol.v'first..conpar-1);
res(conpar) := sol.t;
res(conpar+1..sol.n+1) := sol.v(conpar..sol.v'last);
return res;
end Convert;
function Convert ( sols : Solution_List; conpar : natural )
return Standard_Complex_VecVecs.VecVec is
-- DESCRIPTION :
-- Converts the solutions in sols to the vector representation.
res : Standard_Complex_VecVecs.VecVec(1..Length_Of(sols));
tmp : Solution_List := sols;
begin
for i in res'range loop
declare
ls : Link_to_Solution := Head_Of(tmp);
v : Standard_Complex_Vectors.Vector(ls.v'first..ls.v'last+1);
begin
v := Convert(ls.all,conpar);
res(i) := new Standard_Complex_Vectors.Vector'(v);
end;
tmp := Tail_Of(tmp);
end loop;
return res;
end Convert;
function Square ( n : natural; p : Poly_Sys ) return Poly_Sys is
-- DESCRIPTION :
-- Returns a n-by-n system, by adding random linear combinations of
-- the polynomials p(i), i > n, to the first n polynomials.
res : Poly_Sys(1..n);
acc : Poly;
begin
for i in res'range loop
Copy(p(i),res(i));
for j in n+1..p'last loop
acc := Random1*p(j);
Add(res(i),acc);
Clear(acc);
end loop;
end loop;
return res;
end Square;
procedure Refine_Roots
( file : in file_type; p : in Poly_Sys;
sols : in out Solution_List; report : in boolean ) is
epsxa : constant double_float := 10.0**(-8);
epsfa : constant double_float := 10.0**(-8);
tolsing : constant double_float := 10.0**(-8);
max : constant natural := 3;
numit : natural := 0;
begin
if report
then
Reporting_Root_Refiner(file,p,sols,epsxa,epsfa,tolsing,numit,max,false);
else
Silent_Root_Refiner(p,sols,epsxa,epsfa,tolsing,numit,max);
end if;
end Refine_Roots;
procedure Determinantal_Polynomial_Continuation
( file : in file_type; lochom : in Poly_Sys; -- out added
locsol : in out Standard_Complex_Vectors.Vector;
report,outlog : in boolean ) is
-- DESCRIPTION :
-- Given a homotopy and solution with as last variable the continuation
-- parameter the polynomial continuation will launched.
-- This homotopy uses the determinantal expansions defined in the paper.
-- NOTE :
-- This routine is only used in the original Pieri homotopy algorithm.
sols : Solution_List := Create(locsol);
thesys,squsys : Poly_Sys(locsol'first..locsol'last-1);
square_case : boolean;
-- ran : Complex_Number;
begin
-- for i in lochom'range loop -- added to deal with real input
-- ran := Random1;
-- Mul(lochom(i),ran);
-- end loop;
if lochom'length + 1 = locsol'length
then square_case := true;
Homotopy.Create(lochom,locsol'last);
else square_case := false;
squsys := Square(locsol'last-1,lochom);
if outlog
then put_line(file,"The homotopy as square system : ");
put_line(file,squsys);
end if;
Homotopy.Create(squsys,locsol'last);
end if;
declare
procedure Sil_Cont is
new Silent_Continue(Max_Norm,
Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
procedure Rep_Cont is
new Reporting_Continue(Max_Norm,
Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
begin
if report
then Rep_Cont(file,sols,false,Create(1.0));
else Sil_Cont(sols,false,Create(1.0));
end if;
end;
if square_case
then thesys := Eval(lochom,Create(1.0),locsol'last);
Refine_Roots(file,thesys,sols,report);
else thesys := Eval(squsys,Create(1.0),locsol'last);
Refine_Roots(file,thesys,sols,report);
Clear(squsys);
end if;
Clear(thesys);
locsol := Convert(sols);
Clear(sols);
Homotopy.Clear;
end Determinantal_Polynomial_Continuation;
procedure Determinantal_Polynomial_Continuation
( file : in file_type; lochom : in Poly_Sys;
conpar : in natural;
locsols : in out Standard_Complex_VecVecs.VecVec;
report,outlog : in boolean ) is
-- DESCRIPTION :
-- Given a homotopy and solution with as continuation parameter
-- the variable indexed by conpar, the continuation will be launched.
-- This homotopy uses the determinantal expansions defined in the paper.
sols : Solution_List := Create(locsols,conpar);
firstsol : constant Standard_Complex_Vectors.Vector
:= locsols(locsols'first).all;
thesys,squsys : Poly_Sys(firstsol'first..firstsol'last-1);
square_case : boolean;
begin
if lochom'length + 1 = firstsol'length
then square_case := true;
Homotopy.Create(lochom,conpar);
else square_case := false;
squsys := Square(firstsol'last-1,lochom);
if outlog
then put_line(file,"The homotopy as square system : ");
put_line(file,squsys);
end if;
Homotopy.Create(squsys,conpar);
end if;
declare
procedure Sil_Cont is
new Silent_Continue(Max_Norm,
Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
procedure Rep_Cont is
new Reporting_Continue(Max_Norm,
Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
begin
if report
then Rep_Cont(file,sols,false,Create(1.0));
else Sil_Cont(sols,false,Create(1.0));
end if;
end;
if square_case
then thesys := Eval(lochom,Create(1.0),conpar);
Refine_Roots(file,thesys,sols,report);
else thesys := Eval(squsys,Create(1.0),conpar);
Refine_Roots(file,thesys,sols,report);
Clear(squsys);
end if;
Clear(thesys);
locsols := Convert(sols,conpar);
Clear(sols);
Homotopy.Clear;
end Determinantal_Polynomial_Continuation;
-- TARGET ROUTINES :
procedure Trace_Paths ( file : in file_type; homsys : in Poly_Sys;
locmap : in Standard_Natural_Matrices.Matrix;
report,outlog : in boolean;
plane : in out Standard_Complex_Matrices.Matrix ) is
-- NOTE :
-- This routine is only used in the original Pieri homotopy algorithm.
plavec : constant Standard_Complex_Vectors.Vector := Vector_Rep(plane);
solvec : Standard_Complex_Vectors.Vector(1..plavec'last+1);
lochom : Poly_Sys(homsys'range) := Localize(locmap,homsys);
plaloc : constant Standard_Complex_Matrices.Matrix
:= Localize(locmap,plane);
solloc : constant Standard_Complex_Vectors.Vector
:= Vector_Rep(locmap,plaloc);
locsol : Standard_Complex_Vectors.Vector(1..solloc'last+1);
evahom : Standard_Complex_Vectors.Vector(homsys'range);
evaloc : Standard_Complex_Vectors.Vector(lochom'range);
begin
if outlog
then put_line(file,"The localization pattern :"); put(file,locmap);
put_line(file,"The localized homotopy : "); put_line(file,lochom);
end if;
-- put_line(file,"The plane in local coordinates : "); put(file,plaloc,2);
solvec(plavec'range) := plavec;
solvec(solvec'last) := Create(0.0);
-- put_line(file,"The solution vector : "); put_line(file,solvec);
evahom := Eval(homsys,solvec);
-- put_line(file,"Evaluation of solution vector at homotopy for t=0 :");
-- put_line(file,evahom);
put(file,"Residual of start solution :");
put(file,Max_Norm(evahom),3); put_line(file,".");
-- put_line(file,"The localized vector representation of the plane at t=0 :");
-- put_line(file,solloc);
locsol(solloc'range) := solloc;
locsol(locsol'last) := Create(0.0);
evaloc := Eval(lochom,locsol);
put(file,"Residual of localized plane at t=0 :");
put(file,Max_Norm(evaloc),3); put_line(file,".");
-- Linear_Polynomial_Continuation(file,lochom,locsol,report,outlog);
Determinantal_Polynomial_Continuation(file,lochom,locsol,report,outlog);
-- put_line(file,"The localized vector representation of the plane at t=1 :");
-- put_line(file,locsol);
evaloc := Eval(lochom,locsol);
-- put_line(file,"Evaluation of localized vector at homotopy for t=1 :");
-- put_line(file,evaloc);
put(file,"Residual of localized plane at t=1 :");
put(file,Max_Norm(evaloc),3); put_line(file,".");
plane := Matrix_Rep(locmap,locsol(1..locsol'last-1));
-- put_line(file,"The solution plane at t=1 :");
-- put(file,plane,2);
solvec(plavec'range) := Vector_Rep(plane);
solvec(solvec'last) := Create(1.0);
evahom := Eval(homsys,solvec);
-- put_line(file,"Evaluation of solution vector at homotopy for t=1 :");
-- put_line(file,evahom);
put(file,"Residual of target solution :");
put(file,Max_Norm(evahom),3); put_line(file,".");
Clear(lochom);
end Trace_Paths;
procedure Trace_Paths
( file : in file_type; homsys : in Poly_Sys;
locmap : in Standard_Natural_Matrices.Matrix;
report,outlog : in boolean;
planes : in Standard_Complex_VecMats.VecMat ) is
lochom : Poly_Sys(homsys'range) := Localize(locmap,homsys);
locsols : Standard_Complex_VecVecs.VecVec(planes'range);
evaloc : Standard_Complex_Vectors.Vector(lochom'range);
evahom : Standard_Complex_Vectors.Vector(homsys'range);
lastvar : natural;
begin
if outlog
then put_line(file,"The localized homotopy :");
put_line(file,lochom);
end if;
for i in planes'range loop -- solution planes into vectors
declare
plaloc : constant Standard_Complex_Matrices.Matrix
:= Localize(locmap,planes(i).all);
solloc : constant Standard_Complex_Vectors.Vector
:= Vector_Rep(locmap,plaloc);
locsol : Standard_Complex_Vectors.Vector(1..solloc'last+1);
begin
locsol(solloc'range) := solloc;
locsol(locsol'last) := Create(0.0);
locsols(i) := new Standard_Complex_Vectors.Vector'(locsol);
if outlog
then evaloc := Eval(lochom,locsol);
put(file,"Residual of localized plane at start :");
put(file,Max_Norm(evaloc),3); put_line(file,".");
end if;
end;
end loop;
lastvar := locsols(locsols'first)'last;
Determinantal_Polynomial_Continuation
(file,lochom,lastvar,locsols,report,outlog);
for i in locsols'range loop -- solution vectors into planes
planes(i).all := Matrix_Rep(locmap,locsols(i)(1..locsols(i)'last-1));
if outlog
then
evaloc := Eval(lochom,locsols(i).all);
put(file,"Residual of localized plane at target :");
put(file,Max_Norm(evaloc),3); put_line(file,".");
declare
plavec : constant Standard_Complex_Vectors.Vector
:= Vector_Rep(planes(i).all);
solvec : Standard_Complex_Vectors.Vector
(plavec'first..plavec'last+1);
begin
solvec(plavec'range) := plavec;
solvec(solvec'last) := Create(1.0);
evahom := Eval(homsys,solvec);
put(file,"Residual of target solution :");
put(file,Max_Norm(evahom),3); put_line(file,".");
end;
end if;
end loop;
Clear(lochom); Clear(locsols);
end Trace_Paths;
procedure Quantum_Trace_Paths
( file : in file_type; m,p,q : in natural; nd : in Node;
-- (m,p,q) are needed for the symbol table
homsys : in Poly_Sys; conpar,s_mode : in natural;
locmap : in Standard_Natural_Matrices.Matrix;
report,outlog : in boolean;
planes : in Standard_Complex_VecMats.VecMat ) is
lochom : Poly_Sys(homsys'range)
:= Column_Localize(nd.top,nd.bottom,locmap,homsys);
locsols : Standard_Complex_VecVecs.VecVec(planes'range);
evaloc : Standard_Complex_Vectors.Vector(lochom'range);
evahom : Standard_Complex_Vectors.Vector(homsys'range);
addmix : natural;
begin
if outlog
then put_line(file,"The localization map : "); put(file,locmap);
if nd.tp = mixed
then Two_Set_up_Symbol_Table(m,p,q,nd.top,nd.bottom);
else One_Set_up_Symbol_Table(m,p,q,nd.top,nd.bottom);
end if;
Reduce_Symbols(nd.top,nd.bottom,locmap);
put_line(file,"The localized homotopy : "); put_line(file,lochom);
end if;
if nd.tp = mixed
then addmix := 1;
else addmix := 0;
end if;
for i in planes'range loop -- solution planes into vectors
declare
plaloc : constant Standard_Complex_Matrices.Matrix
:= Localize(locmap,planes(i).all);
solloc : constant Standard_Complex_Vectors.Vector
:= Column_Vector_Rep(locmap,plaloc);
locsol : Standard_Complex_Vectors.Vector(1..solloc'last+2+addmix);
begin
locsol(solloc'range) := solloc;
if s_mode = 0
then locsol(locsol'last-1) := Create(0.0); -- start value for s is 0
else locsol(locsol'last-1) := Create(1.0); -- start value for s is 1
end if;
locsol(locsol'last) := Create(0.0); -- start value for t is 0
if addmix = 1
then locsol(locsol'last-2) := Create(1.0); -- additional s-value
end if;
locsols(i) := new Standard_Complex_Vectors.Vector'(locsol);
if outlog
then evaloc := Eval(lochom,locsol);
put(file,"Residual of localized plane at start :");
put(file,Max_Norm(evaloc),3); put_line(file,".");
end if;
end;
end loop;
Determinantal_Polynomial_Continuation
(file,lochom,conpar,locsols,report,outlog);
for i in locsols'range loop -- solution vectors into planes
planes(i).all
:= Column_Matrix_Rep(locmap,locsols(i)(1..locsols(i)'last-2-addmix));
if outlog
then evaloc := Eval(lochom,locsols(i).all);
put(file,"Residual of localized plane at target :");
put(file,Max_Norm(evaloc),3); put_line(file,".");
put_line(file,"The computed plane : ");
put(file,planes(i).all,2);
end if;
declare
plavec : constant Standard_Complex_Vectors.Vector
:= Column_Vector_Rep(nd.top,nd.bottom,planes(i).all);
solvec : Standard_Complex_Vectors.Vector
(plavec'first..plavec'last+2+addmix);
begin
solvec(plavec'range) := plavec;
solvec(solvec'last) := Create(1.0); -- target value for t
solvec(solvec'last-1) := locsols(i)(locsols(i)'last-1);
if addmix = 1
then solvec(solvec'last-2) := locsols(i)(locsols(i)'last-2);
end if;
if outlog
then evahom := Eval(homsys,solvec);
put(file,"Residual of target solution :");
put(file,Max_Norm(evahom),3); put_line(file,".");
end if;
end;
end loop;
Clear(lochom); Clear(locsols);
end Quantum_Trace_Paths;
end Pieri_Continuation;