File: [local] / OpenXM_contrib / PHC / Ada / Schubert / deformation_posets.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 unchecked_deallocation;
with integer_io; use integer_io;
with Timing_Package; use Timing_Package;
with Characters_and_Numbers; use Characters_and_Numbers;
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_Complex_Numbers_io; use Standard_Complex_Numbers_io;
with Standard_Complex_Vectors_io; use Standard_Complex_Vectors_io;
with Standard_Natural_Matrices;
with Standard_Complex_Matrices;
with Standard_Complex_Matrices_io; use Standard_Complex_Matrices_io;
with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
with Symbol_Table; use Symbol_Table;
with Standard_Complex_Poly_Functions; use Standard_Complex_Poly_Functions;
with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
with Standard_Complex_Poly_Matrices;
with Standard_Complex_Poly_Matrices_io; use Standard_Complex_Poly_Matrices_io;
with Brackets,Brackets_io; use Brackets,Brackets_io;
with Bracket_Monomials; use Bracket_Monomials;
with Bracket_Polynomials; use Bracket_Polynomials;
with Bracket_Systems; use Bracket_Systems;
with Bracket_Systems_io; use Bracket_Systems_io;
with Plane_Representations; use Plane_Representations;
with Symbolic_Minor_Equations; use Symbolic_Minor_Equations;
with Numeric_Minor_Equations; use Numeric_Minor_Equations;
with Determinantal_Systems; use Determinantal_Systems;
with Specialization_of_Planes; use Specialization_of_Planes;
with Curves_into_Grassmannian; use Curves_into_Grassmannian;
with Curves_into_Grassmannian_io; use Curves_into_Grassmannian_io;
with Pieri_Homotopies; use Pieri_Homotopies;
with Pieri_Continuation; use Pieri_Continuation;
package body Deformation_Posets is
-- BRACKET AUXILIARITIES TO DETERMINE PIVOTS :
function Complement ( n : natural; b : Bracket ) return Bracket is
-- DESCRIPTION :
-- Returns the complement of the bracket b, defined as a bracket
-- of range 1..n-b'length as an ordered subset of {1..n} \ b.
res : Bracket(1..n-b'last);
cnt : natural := 0;
ind : natural := 1;
begin
for i in 1..n loop
if ((ind > b'last) or else (i < b(ind)))
then cnt := cnt+1;
res(cnt) := i;
elsif i = b(ind)
then ind := ind+1;
end if;
end loop;
return res;
end Complement;
function Remove ( b : Bracket; l : natural ) return Bracket is
-- DESCRIPTION :
-- Returns a smaller bracket that does not contain l.
-- REQUIRED : there exists a k: b(k) = l.
res : Bracket(1..b'last-1);
cnt : natural := 0;
begin
for i in b'range loop
if b(i) /= l
then cnt := cnt+1;
res(cnt) := b(i);
end if;
end loop;
return res;
end Remove;
function Is_In ( b : Bracket; l : natural ) return boolean is
-- DESCRIPTION :
-- Returns true if there exists an index k such that b(k) = l.
begin
for k in b'range loop
if b(k) = l
then return true;
end if;
end loop;
return false;
end Is_In;
function Remove ( b1,b2 : Bracket ) return Bracket is
-- DESCRIPTION :
-- Returns b1 minus the first element in b2 that also occurs in b1.
begin
for i in b2'range loop
if Is_In(b1,b2(i))
then return Remove(b1,b2(i));
end if;
end loop;
return b1;
end Remove;
function Remove ( cols,b,subb : Bracket ) return Bracket is
-- DESCRIPTION :
-- The indices in cols correspond to the entries in b.
-- The bracket subb is a sub-bracket of b, with only one entry removed.
-- The indices on return correspond to the entries in subb.
res : Bracket(subb'range);
begin
for i in subb'range loop
if b(i) = subb(i)
then res(i) := cols(i);
else res(i) := cols(i+1);
end if;
end loop;
return res;
end Remove;
-- POSET-ORIENTED PIERI DEFORMATIONS :
function Leaf_Plane ( n : natural; nd : Node )
return Standard_Complex_Matrices.Matrix is
-- DESCRIPTION :
-- Returns the solution plane that corresponds to a leaf of the poset.
res : Standard_Complex_Matrices.Matrix(1..n,nd.top'range);
begin
for i in res'range(1) loop
for j in res'range(2) loop
res(i,j) := Create(0.0);
end loop;
end loop;
for i in nd.top'range loop
res(nd.top(i),i) := Create(1.0);
end loop;
return res;
end Leaf_Plane;
-- DEFORMATIONS :
function Path_Coordinates ( level,label,child : natural )
return string is
begin
return "tracing (User time) at node("
& Convert(level) & ")("
& Convert(label) & ") from child "
& Convert(child);
end Path_Coordinates;
procedure Write_Path_Coordinates
( file : in file_type;
level,label,path,child,childpath : in natural ) is
-- DESCRIPTION :
-- Writes all coordinates from the current path that is to be traced.
begin
put(file,"Tracing at node("); put(file,level,1); put(file,")(");
put(file,label,1); put(file,") path "); put(file,path,1);
put(file," as path "); put(file,childpath,1);
put(file," from child "); put(file,child,1); new_line(file);
end Write_Path_Coordinates;
procedure Write_Path_Coordinates
( file : in file_type;
level,label,child : in natural ) is
-- DESCRIPTION :
-- Writes all coordinates from the current path that is to be traced.
begin
put(file,"Tracing paths at node("); put(file,level,1); put(file,")(");
put(file,label,1);
put(file,") from child "); put(file,child,1); new_line(file);
end Write_Path_Coordinates;
procedure Deform_from_Children
( file : in file_type;
poset : in out Array_of_Array_of_VecMats;
nd : in Node; n,uplevel : in natural;
homotopy : in Poly_Sys; report,outlog : in boolean;
x : in Standard_Complex_Poly_Matrices.Matrix;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Deforms from the i-th non-empty child of nd at uplevel.
-- ON ENTRY :
-- file to write intermediate results on;
-- poset poset of solution p-planes;
-- nd current node in the localization poset;
-- n dimension of the working space;
-- uplevel level where to find the start planes in the poset;
-- homotopy family of moving planes;
-- report indicates whether intermediate output during continuation;
-- outlog flag to write homotopies on file if set to true.
-- x matrix of unknowns according to a localization pattern.
-- ON RETURN :
-- poset updated poset of solution p-planes;
-- npaths updated number of paths traced at that level;
-- timings updated elapsed user timings.
locmap : Standard_Natural_Matrices.Matrix(1..n,1..nd.p);
solcnt : natural := 0;
label : natural;
begin
for i in nd.child_labels'range loop
label := nd.child_labels(i);
if not Empty(poset,uplevel,label) -- child.roco > 0
then
Write_Path_Coordinates(file,nd.level,nd.label,label);
declare
planes : VecMat(poset(uplevel)(label)'range);
timer : Timing_Widget;
begin
tstart(timer);
for i in planes'range loop -- create to avoid sharing
planes(i) := new Standard_Complex_Matrices.Matrix'
(poset(uplevel)(label)(i).all);
end loop;
locmap := Standard_Coordinate_Frame(x,planes(planes'first).all);
Trace_Paths(file,homotopy,locmap,report,outlog,planes);
for i in planes'range loop
solcnt := solcnt+1;
poset(nd.level)(nd.label)(solcnt) := planes(i);
end loop;
tstop(timer);
new_line(file);
print_times(file,timer,Path_Coordinates(nd.level,nd.label,label));
new_line(file);
timings(nd.level) := timings(nd.level) + Elapsed_User_Time(timer);
end;
end if;
end loop;
npaths(nd.level) := npaths(nd.level) + solcnt;
end Deform_from_Children;
procedure Quantum_Deform_from_Children
( file : in file_type;
poset : in out Array_of_Array_of_VecMats;
nd : in Node; n,q,uplevel : in natural;
homotopy : in Poly_Sys; conpar,s_mode : in natural;
report,outlog : in boolean;
x : in Standard_Complex_Poly_Matrices.Matrix;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Deforms from the i-th non-empty child of nd at uplevel.
-- This is the quantum analogue to the hypersurface Pieri.
-- ON ENTRY :
-- file to write intermediate results on;
-- poset poset of solution p-planes;
-- nd current node in the localization poset;
-- n dimension of the working space;
-- q degree of the curve;
-- uplevel level where to find the start planes in the poset;
-- homotopy family of moving planes;
-- conpar number of the continuation parameter;
-- s_mode if = 0, then s = 0, otherwise s = 1 at start;
-- report indicates whether intermediate output during continuation;
-- outlog flag to write homotopies on file if set to true;
-- x symbolic representation of the curve matrix of polynomials.
-- ON RETURN :
-- poset updated poset of solution p-planes;
-- npaths updated number of paths at each level;
-- timings updated CPU user timings for each level.
m : constant natural := n - nd.p;
rws : constant natural := n*(q+1);
locmap : Standard_Natural_Matrices.Matrix(1..rws,1..nd.p);
solcnt : natural := 0;
label : natural;
begin
for i in nd.child_labels'range loop
label := nd.child_labels(i);
if not Empty(poset,uplevel,label) -- child.roco > 0
then
Write_Path_Coordinates(file,nd.level,nd.label,label);
declare
planes : VecMat(poset(uplevel)(label)'range);
timer : Timing_Widget;
begin
tstart(timer);
for i in planes'range loop -- create to avoid sharing
planes(i) := new Standard_Complex_Matrices.Matrix'
(poset(uplevel)(label)(i).all);
end loop;
locmap := Standard_Coordinate_Frame
(m,nd.p,q,nd.top,nd.bottom,planes(planes'first).all);
Quantum_Trace_Paths
(file,m,nd.p,q,nd,homotopy,conpar,s_mode,locmap,
report,outlog,planes);
for i in planes'range loop
solcnt := solcnt+1;
poset(nd.level)(nd.label)(solcnt) := planes(i);
end loop;
tstop(timer);
new_line(file);
print_times(file,timer,Path_Coordinates(nd.level,nd.label,label));
new_line(file);
timings(nd.level) := timings(nd.level) + Elapsed_User_Time(timer);
end;
end if;
end loop;
npaths(nd.level) := npaths(nd.level) + solcnt;
end Quantum_Deform_from_Children;
procedure Hypersurface_Deform
( file : in file_type; n : in natural;
poset : in out Array_of_Array_of_VecMats;
nd : in Node; expbp : in Bracket_Polynomial;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Does the Pieri deformations to the node, same specifications as Solve
-- for the hypersurface case.
-- REQUIRED : nd.level > 0, solutions at children of nd are in poset.
xpm : Standard_Complex_Poly_Matrices.Matrix(1..n,1..nd.p)
:= Localization_Pattern(n,nd.top,nd.bottom);
homsys : Poly_Sys(1..nd.level);
begin
if nd.tp = mixed
then homsys := Two_Hypersurface_Pieri_Homotopy(n,nd,expbp,xpm,planes);
Deform_from_Children
(file,poset,nd,n,nd.level-2,homsys,report,outlog,
xpm,npaths,timings);
else homsys := One_Hypersurface_Pieri_Homotopy(n,nd,expbp,xpm,planes);
Deform_from_Children
(file,poset,nd,n,nd.level-1,homsys,report,outlog,
xpm,npaths,timings);
end if;
Standard_Complex_Poly_Matrices.Clear(xpm);
Clear(homsys);
end Hypersurface_Deform;
procedure One_General_Deform
( file : in file_type; n,ind : in natural;
poset : in out Array_of_Array_of_VecMats; nd : in Node;
start,target : in Standard_Complex_Matrices.Matrix;
planes : in VecMat; bs : in Bracket_System;
report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Does the Pieri deformations to the node, same specifications as Solve
-- for the general case.
-- REQUIRED : nd.level > 0, solutions at children of nd are in poset.
-- ON ENTRY :
-- file to write intermediate output to;
-- n number of rows in the matrices;
-- ind planes(ind) is currently being folded in with this chain;
-- poset contains solution planes at higher levels;
-- nd current node in the localization poset;
-- start start (m+1-k)-plane for pivots;
-- target target (m+1-k)-plane for pivots;
-- planes target planes;
-- bs structure to expand the minors;
-- report switch to determine output during continuation;
-- outlog flag to write homotopies on file if set to true.
-- ON RETURN :
-- poset solution planes at (nd.level)(nd.label) are determined;
-- npaths number of paths followed at each level;
-- timings CPU user time at each level.
xpm : Standard_Complex_Poly_Matrices.Matrix(1..n,1..nd.p)
:= Localization_Pattern(n,nd.top,nd.bottom);
hom : Link_to_Poly_Sys
:= One_General_Pieri_Homotopy(n,ind,nd,bs,start,target,xpm,planes);
begin
Deform_from_Children
(file,poset,nd,n,nd.level-1,hom.all,report,outlog,xpm,npaths,timings);
Standard_Complex_Poly_Matrices.Clear(xpm);
Clear(hom);
end One_General_Deform;
procedure Two_General_Deform
( file : in file_type; n,ind : in natural;
poset : in out Array_of_Array_of_VecMats; nd : in Node;
top_start,top_target,bot_start,bot_target
: in Standard_Complex_Matrices.Matrix;
planes : in VecMat; top_bs,bot_bs : in Bracket_System;
report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Does the Pieri deformations to the node, same specifications as Solve
-- for the general case.
-- REQUIRED : nd.level > 0, solutions at children of nd are in poset.
-- ON ENTRY :
-- file to write intermediate output to;
-- n number of rows in the matrices;
-- ind planes(ind) is currently being folded in with this chain;
-- poset contains solution planes at higher levels;
-- nd current node in the localization poset;
-- top_start start (m+1-k)-plane for top pivots;
-- top_target target (m+1-k)-plane for top pivots;
-- bot_start start (m+1-k)-plane for bottom pivots;
-- bot_target target (m+1-k)-plane for bottom pivots;
-- planes target planes;
-- top_bs structure to expand the minors for top pivots;
-- bot_bs structure to expand the minors for bottom pivots;
-- report switch to determine output during continuation;
-- outlog flag to write homotopies on file if set to true.
-- ON RETURN :
-- poset solution planes at (nd.level)(nd.label) are determined;
-- npaths number of paths traced at each level;
-- timings updated CPU user times for each level.
xpm : Standard_Complex_Poly_Matrices.Matrix(1..n,1..nd.p)
:= Localization_Pattern(n,nd.top,nd.bottom);
homotopy : Link_to_Poly_Sys;
begin
case nd.tp is
when top
=> homotopy := One_General_Pieri_Homotopy
(n,ind,nd,top_bs,top_start,top_target,xpm,planes);
Deform_from_Children
(file,poset,nd,n,nd.level-1,homotopy.all,report,outlog,xpm,
npaths,timings);
when bottom
=> homotopy := One_General_Pieri_Homotopy
(n,ind,nd,bot_bs,bot_start,bot_target,xpm,planes);
Deform_from_Children
(file,poset,nd,n,nd.level-1,homotopy.all,report,outlog,xpm,
npaths,timings);
when mixed
=> homotopy := Two_General_Pieri_Homotopy
(n,ind,nd,top_bs,bot_bs,top_start,top_target,
bot_start,bot_target,xpm,planes);
Deform_from_Children
(file,poset,nd,n,nd.level-2,homotopy.all,report,outlog,xpm,
npaths,timings);
end case;
Standard_Complex_Poly_Matrices.Clear(xpm);
Clear(homotopy);
end Two_General_Deform;
procedure Quantum_Deform
( file : in file_type; n,q : in natural;
poset : in out Array_of_Array_of_VecMats;
nd : in Node; expbp : in Bracket_Polynomial;
planes : in VecMat; s : Standard_Complex_Vectors.Vector;
report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- This is the q-analogue to the Hypersurface Deform.
-- REQUIRED : nd.level > 0, solutions at children of nd are in poset.
m : constant natural := n-nd.p;
xpm : Standard_Complex_Poly_Matrices.Matrix(1..n,1..nd.p)
:= Symbolic_Create(m,nd.p,q,nd.top,nd.bottom);
begin
if outlog
then put(file,"Curve at node(");
put(file,nd.level,1); put(file,")("); put(file,nd.label,1);
put(file,") for pivots (");
put(file,nd.top); put(file,","); put(file,nd.bottom);
put_line(file,") :");
One_Set_up_Symbol_Table(m,nd.p,q,nd.top,nd.bottom);
put(file,xpm);
end if;
if nd.tp = mixed
then
declare
homsys : Poly_Sys(1..nd.level+2);
begin
homsys := Two_Quantum_Pieri_Homotopy(n,nd,expbp,xpm,planes,s);
if outlog
then Two_Set_up_Symbol_Table(m,nd.p,q,nd.top,nd.bottom);
put_line(file,"The homotopy : "); put_line(file,homsys);
end if;
Quantum_Deform_from_Children
(file,poset,nd,n,q,nd.level-2,homsys,nd.level+3,1,
report,outlog,xpm,npaths,timings);
Clear(homsys);
end;
else
declare
homsys : Poly_Sys(1..nd.level+1);
begin
homsys := One_Quantum_Pieri_Homotopy(n,nd,expbp,xpm,planes,s);
if outlog
then put_line(file,"The homotopy : "); put_line(file,homsys);
end if;
Quantum_Deform_from_Children
(file,poset,nd,n,q,nd.level-1,homsys,nd.level+2,1,
report,outlog,xpm,npaths,timings);
Clear(homsys);
end;
end if;
Standard_Complex_Poly_Matrices.Clear(xpm);
end Quantum_Deform;
function Moving_Point_Mode
( l,k : natural; modpiv : Bracket ) return natural is
-- DESCRIPTION :
-- Returns a natural number that indicates the moving of the
-- interpolation point. The value on return means the following
-- when = 0 : s goes from 0 to 1;
-- = 1 : s remains constant at 1;
-- = 2 : s goes from 1 to a target value.
begin
if l = 0
then return 2;
elsif l = k-1 and modpiv(1) > 1
then return 0;
else return 1;
end if;
end Moving_Point_Mode;
procedure One_General_Quantum_Deform
( file : in file_type; n,q,l,k,ind : in natural;
poset : in out Array_of_Array_of_VecMats; nd : in Node;
modpiv : in Bracket;
start,target : in Standard_Complex_Matrices.Matrix;
planes : in VecMat; s : in Standard_Complex_Vectors.Vector;
bs : in Bracket_System; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- This is the quantum analogue to One_General_Deform.
-- REQUIRED : nd.level > 0, solutions at children of nd are in poset.
-- ON ENTRY :
-- file to write intermediate output to;
-- n number of rows in the matrices;
-- q degree of the map;
-- l runs from k-1 downto 0;
-- k co-dimension currently being treated;
-- ind planes(ind) is currently being folded in with this chain;
-- poset contains solution planes at higher levels;
-- nd current node in the localization poset;
-- modpiv bottom or top pivots modulo n;
-- start start (m+1-k)-plane for pivots;
-- target target (m+1-k)-plane for pivots;
-- planes target planes;
-- s interpolation points where the maps meets the planes;
-- bs structure to expand the minors;
-- report switch to determine output during continuation;
-- outlog flag to write homotopies on file if set to true.
-- ON RETURN :
-- poset solution planes at (nd.level)(nd.label) are determined;
-- npaths number of paths followed at each level;
-- timings CPU user time at each level.
m : constant natural := n-nd.p;
xpm : Standard_Complex_Poly_Matrices.Matrix(1..n,1..nd.p)
:= Symbolic_Create(m,nd.p,q,nd.top,nd.bottom);
s_mode : constant natural := Moving_Point_Mode(l,k,modpiv);
hom : Link_to_Poly_Sys
:= One_General_Quantum_Pieri_Homotopy
(n,ind,nd,s_mode,bs,start,target,xpm,planes,s);
begin
if outlog
then put(file,"level l : "); put(file,l,1); put(file," ");
put(file,"codim k : "); put(file,k,1); new_line(file);
put(file,"Curve at node(");
put(file,nd.level,1); put(file,")("); put(file,nd.label,1);
put(file,") for pivots (");
put(file,nd.top); put(file,","); put(file,nd.bottom);
put_line(file,") :");
One_Set_up_Symbol_Table(m,nd.p,q,nd.top,nd.bottom);
put(file,xpm);
put_line(file,"The homotopy : "); put_line(file,hom.all);
end if;
Quantum_Deform_from_Children
(file,poset,nd,n,q,nd.level-1,hom.all,nd.level+2,s_mode,
report,outlog,xpm,npaths,timings);
Standard_Complex_Poly_Matrices.Clear(xpm);
Clear(hom);
end One_General_Quantum_Deform;
-- CREATORS :
function Create ( index_poset : Array_of_Array_of_Nodes )
return Array_of_Array_of_VecMats is
res : Array_of_Array_of_VecMats(index_poset'range);
lnd : Link_to_Node;
begin
for i in index_poset'range loop
if index_poset(i) /= null
then res(i) := new Array_of_VecMats(index_poset(i)'range);
for j in res(i)'range loop
lnd := index_poset(i)(j);
if lnd.roco /= 0
then res(i)(j) := new VecMat(1..lnd.roco);
end if;
end loop;
end if;
end loop;
return res;
end Create;
-- SELECTORS :
function Empty ( poset : Array_of_Array_of_VecMats; level,label : natural )
return boolean is
use Standard_Complex_Matrices;
begin
if poset(level) = null
then return true;
elsif poset(level)(label) = null
then return true;
else declare
lavm : Link_to_VecMat := poset(level)(label);
begin
return (lavm(lavm'first) = null);
end;
end if;
end Empty;
-- ANALOGUES TO THE ROOT COUNTERS :
procedure Recursive_Hypersurface_Solve
( file : in file_type; n : in natural;
nd : in Node; expbp : in Bracket_Polynomial;
poset : in out Array_of_Array_of_VecMats;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- This additional layer is added to avoid the repeated construction
-- of the structure of the equations, that is now in expbp.
begin
if Empty(poset,nd.level,nd.label)
then
if nd.level = 0
then
poset(nd.level)(nd.label)(1)
:= new Standard_Complex_Matrices.Matrix'(Leaf_Plane(n,nd));
else
for i in nd.children'range(1) loop
for j in nd.children'range(2) loop
if nd.children(i,j) /= null
then Recursive_Hypersurface_Solve
(file,n,nd.children(i,j).all,expbp,
poset,planes,report,outlog,npaths,timings);
end if;
end loop;
end loop;
Hypersurface_Deform
(file,n,poset,nd,expbp,planes,report,outlog,npaths,timings);
end if;
end if;
end Recursive_Hypersurface_Solve;
procedure Solve ( file : in file_type; n : in natural;
poset : in out Array_of_Array_of_VecMats;
nd : in Node; planes : in VecMat;
report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
bm : Bracket_Monomial := Maximal_Minors(n,n);
bs : Bracket_System(0..Number_of_Brackets(bm))
:= Minor_Equations(n,n-nd.p,bm);
begin
Recursive_Hypersurface_Solve
(file,n,nd,bs(1),poset,planes,report,outlog,npaths,timings);
Clear(bm); Clear(bs);
end Solve;
procedure One_Solve_along_Chains
( file : in file_type; nd : in Node; n,l,k,ind : in natural;
poset : in out Array_of_Array_of_VecMats;
pivots,columns : in Bracket; bs : in Bracket_System;
special,start,target : in Standard_Complex_Matrices.Matrix;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Applies the general solver along the nodes in all chains that end at
-- the current node. This is the analogue to the hypersurface solver,
-- for use in connection with the grandchildren first recursive solving.
-- This procedure is only called in "One_Solve".
-- ON ENTRY :
-- file to write intermediate results on;
-- nd current node in the localization poset;
-- n working dimension, equation m+p;
-- l runs from 0 to k-1;
-- k current codimension condition;
-- poset structure with all solution p-planes;
-- ind ind-1 planes are already folded in;
-- pivots pivot elements used for the special m-plane;
-- columns which columns of the special m-plane are used;
-- bs Laplace expansion of the polynomial equations;
-- special special m-plane for top pivots;
-- start (m+1-k)-plane used at the start of the deformation;
-- target (m+1-k)-plane used as target;
-- planes sequence of (m+1-k)-planes;
-- report indicates whether intermediate output during continuation;
-- outlog flag to write homotopies on file if set to true.
-- ON RETURN :
-- poset updated structure of all solution p-planes;
-- npaths updated numbers of paths traced at each level;
-- timings updated CPU user times at each level.
m : constant natural := n - nd.p;
new_piv,new_col : Bracket(1..pivots'last-1);
new_start : Standard_Complex_Matrices.Matrix(1..n,start'range(2));
begin
if empty(poset,nd.level,nd.label)
then
if l < k-1
then
for i in nd.children'range(1) loop
for j in nd.children'range(2) loop
if ((nd.children(i,j) /= null)
and then (nd.children(i,j).roco > 0))
then
if nd.children(i,j).tp = top
then new_piv := Remove(pivots,nd.children(i,j).top);
else new_piv := Remove(pivots,nd.children(i,j).bottom);
end if;
new_col := Remove(columns,pivots,new_piv);
new_start := Special_Plane(n,m,k,new_col,special);
One_Solve_along_Chains
(file,nd.children(i,j).all,n,l+1,k,ind,poset,new_piv,
new_col,bs,special,new_start,start,planes,report,outlog,
npaths,timings);
end if;
end loop;
end loop;
end if;
One_General_Deform
(file,n,ind,poset,nd,start,target,planes,bs,report,outlog,
npaths,timings);
end if;
end One_Solve_along_Chains;
procedure One_Quantum_Solve_along_Chains
( file : in file_type; nd : in Node; n,q,l,k,ind : in natural;
poset : in out Array_of_Array_of_VecMats;
pivots,columns : in Bracket; bs : in Bracket_System;
special,start,target : in Standard_Complex_Matrices.Matrix;
planes : in VecMat; s : in Standard_Complex_Vectors.Vector;
report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Applies the general solver along the nodes in all chains that end at
-- the current node. This is the analogue to the hypersurface solver,
-- for use in connection with the grandchildren first recursive solving.
-- This procedure is only called in the q-analogue of "One_Solve".
-- ON ENTRY :
-- file to write intermediate results on;
-- nd current node in the localization poset;
-- n working dimension, equation m+p;
-- q degree of the map;
-- l runs from 0 to k-1;
-- k current codimension condition;
-- poset structure with all solution p-planes;
-- ind ind-1 planes are already folded in;
-- pivots pivot elements used for the special m-plane;
-- columns which columns of the special m-plane are used;
-- bs Laplace expansion of the polynomial equations;
-- special special m-plane for top pivots;
-- start (m+1-k)-plane used at the start of the deformation;
-- target (m+1-k)-plane used as target;
-- planes sequence of (m+1-k)-planes;
-- s interpolation points where the map meets the planes;
-- report indicates whether intermediate output during continuation;
-- outlog flag to write homotopies on file if set to true.
-- ON RETURN :
-- poset updated structure of all solution p-planes;
-- npaths updated numbers of paths traced at each level;
-- timings updated CPU user times at each level.
m : constant natural := n - nd.p;
new_piv,new_col : Bracket(1..pivots'last-1);
new_start : Standard_Complex_Matrices.Matrix(1..n,start'range(2));
mod_piv : Bracket(1..nd.p);
begin
if empty(poset,nd.level,nd.label)
then
if l < k-1
then
for i in nd.children'range(1) loop
for j in nd.children'range(2) loop
if ((nd.children(i,j) /= null)
and then (nd.children(i,j).roco > 0))
then
if nd.children(i,j).tp = top
then mod_piv := Modulo(nd.children(i,j).top,n);
new_piv := Remove(pivots,mod_piv);
put(file,"Top pivots at node : ");
put(file,nd.top);
put(file," child top pivots : ");
put(file,nd.children(i,j).top); new_line(file);
else mod_piv := Modulo(nd.children(i,j).bottom,n);
new_piv := Remove(pivots,mod_piv);
put(file,"Bottom pivots at node : ");
put(file,nd.bottom);
put(file," child bottom pivots : ");
put(file,nd.children(i,j).bottom); new_line(file);
end if;
put(file,"Modular pivots : "); put(file,mod_piv);
put(file," new pivots : "); put(file,new_piv);
new_line(file);
put(file,"Pivot columns : "); put(file,columns);
put(file," new columns : "); put(file,new_piv);
new_line(file);
new_col := Remove(columns,pivots,new_piv);
new_start := Special_Plane(n,m,k,new_col,special);
One_Quantum_Solve_along_Chains
(file,nd.children(i,j).all,n,q,l+1,k,ind,poset,new_piv,
new_col,bs,special,new_start,start,planes,s,report,outlog,
npaths,timings);
end if;
end loop;
end loop;
end if;
One_General_Quantum_Deform
(file,n,q,l,k,ind,poset,nd,mod_piv,start,target,planes,s,bs,
report,outlog,npaths,timings);
end if;
end One_Quantum_Solve_along_Chains;
procedure Solve_along_One_Chain
( file : in file_type; nd : in Node; n,l,k,ind : in natural;
cod : in Bracket; poset : in out Array_of_Array_of_VecMats;
pivots,columns : in Bracket; bs : in Bracket_System;
special,start,target : in Standard_Complex_Matrices.Matrix;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Applies the general solver along the nodes in all chains that end at
-- the current node. This is the analogue to the hypersurface solver,
-- which can be used in a general control structure.
-- ON ENTRY :
-- file to write intermediate results on;
-- nd current node in the localization poset;
-- n working dimension, equation m+p;
-- l runs from 0 to k-1;
-- k current codimension condition;
-- poset structure with all solution p-planes;
-- ind ind-1 planes are already folded in;
-- pivots pivot elements used for the special m-plane;
-- columns which columns of the special m-plane are used;
-- bs Laplace expansion of the polynomial equations;
-- special special m-plane for top pivots;
-- start (m+1-k)-plane used at the start of the deformation;
-- target (m+1-k)-plane used as target;
-- planes sequence of (m+1-k)-planes;
-- report indicates whether intermediate output during continuation;
-- outlog file to write homotopies on if set to true.
-- ON RETURN :
-- poset updated structure of solution p-planes;
-- npaths updated number of paths traced at each level;
-- timings updated CPU user times for each level.
m : constant natural := n - nd.p;
begin
if empty(poset,nd.level,nd.label)
then
if nd.level = 0
then poset(nd.level)(nd.label)(1)
:= new Standard_Complex_Matrices.Matrix'(Leaf_Plane(n,nd));
elsif nd.roco > 0
then
if l = k
then
if cod'last > cod'first
then
declare
kk : constant natural := cod(cod'last-1);
kd : constant natural := n+1-kk;
new_piv,new_col : Bracket(1..m);
new_special : Standard_Complex_Matrices.Matrix(1..n,1..m);
new_target : constant Standard_Complex_Matrices.Matrix
:= planes(cod'last-1).all;
new_start
: Standard_Complex_Matrices.Matrix(1..n,1..m+1-kk);
new_bm : Bracket_Monomial := Maximal_Minors(n,kd);
new_bs : Bracket_System(0..Number_of_Brackets(new_bm))
:= Minor_Equations(kd,kd-nd.p,new_bm);
begin
for i in new_col'range loop
new_col(i) := i;
end loop;
if nd.tp = top
then new_piv := Complement(n,nd.top);
new_special := Special_Top_Plane(m,nd.top);
else new_piv := Complement(n,nd.bottom);
new_special := Special_Bottom_Plane(m,nd.bottom);
end if;
new_start := Special_Plane(n,m,kk,new_col,new_special);
Solve_along_One_Chain
(file,nd,n,0,kk,cod'last-1,cod(cod'first..cod'last-1),
poset,new_piv,new_col,new_bs,new_special,new_start,
new_target,planes,report,outlog,npaths,timings);
Clear(new_bm); Clear(new_bs);
end;
end if;
else
declare
new_piv,new_col : Bracket(1..pivots'last-1);
new_start
: Standard_Complex_Matrices.Matrix(1..n,start'range(2));
begin
for i in nd.children'range(1) loop
for j in nd.children'range(2) loop
if ((nd.children(i,j) /= null)
and then (nd.children(i,j).roco > 0))
then
if nd.children(i,j).tp = top
then new_piv := Remove(pivots,nd.children(i,j).top);
else new_piv
:= Remove(pivots,nd.children(i,j).bottom);
end if;
new_col := Remove(columns,pivots,new_piv);
new_start := Special_Plane(n,m,k,new_col,special);
Solve_along_One_Chain
(file,nd.children(i,j).all,n,l+1,k,ind,cod,poset,
new_piv,new_col,bs,special,new_start,start,
planes,report,outlog,npaths,timings);
end if;
end loop;
end loop;
One_General_Deform
(file,n,ind,poset,nd,start,target,planes,bs,
report,outlog,npaths,timings);
end;
end if;
end if;
end if;
end Solve_along_One_Chain;
procedure Solve_along_Two_Chains
( file : in file_type; nd : in Node;
n,l_top,k_top,l_bot,k_bot,ind : in natural;
cod : in Bracket; poset : in out Array_of_Array_of_VecMats;
top_pivots,top_columns,bot_pivots,bot_columns : in Bracket;
top_bs,bot_bs : in Bracket_System;
top_special,top_start,top_target,bot_special,bot_start,
bot_target : in Standard_Complex_Matrices.Matrix;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array );
-- DESCRIPTION :
-- Applies the general solver along the nodes in all chains that end at
-- the current node. This is the analogue to the hypersurface solver
-- where the type of the node may be anything.
-- ON ENTRY :
-- file to write intermediate results on;
-- nd current node in the localization poset;
-- n working dimension, equation m+p;
-- l_top runs from 0 to k_top-1;
-- k_top co-dimension condition satisfied incrementing top pivots;
-- l_bot runs from 0 to k_bot-1;
-- k_bot co-dimension condition satisfied decrementing bottom pivots;
-- poset structure with all solution p-planes;
-- ind ind-1 planes are already folded in;
-- top_pivots pivot elements used for the special top m-plane;
-- top_columns which columns of the special top m-plane are used;
-- bot_pivots pivot elements used for the special bottom m-plane;
-- bot_columns which columns of the special bottom m-plane are used;
-- top_bs Laplace expansion of the polynomial equations;
-- bot_bs Laplace expansion of the polynomial equations;
-- top_special special m-plane for top pivots;
-- top_start (m+1-k)-plane used at the start of the deformation;
-- top_target (m+1-k)-plane used as target satisfied with top pivots;
-- bot_special special m-plane for top pivots;
-- bot_start (m+1-k)-plane used at the start of the deformation;
-- bot_target (m+1-k)-plane used as target satisfied with bottom pivots;
-- planes sequence of (m+1-k)-planes;
-- report indicates whether intermediate output during continuation;
-- outlog flag to write homotopies on file if set to true.
-- ON RETURN :
-- poset updated structure of solution p-planes;
-- npaths updated number of paths traced at each level;
-- timings updated CPU user timings at each level.
procedure Solve_along_Two_Chains_Deforming_Top_and_Bottom
( file : in file_type; nd : in Node;
n,l_top,k_top,l_bot,k_bot,ind : in natural;
cod : in Bracket; poset : in out Array_of_Array_of_VecMats;
top_pivots,top_columns,bot_pivots,bot_columns : in Bracket;
top_bs,bot_bs : in Bracket_System;
top_special,top_start,top_target,bot_special,bot_start,
bot_target : in Standard_Complex_Matrices.Matrix;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Assumes that k_top > l_top and k_bot > l_bot so that the deformations
-- involve both incrementing top and decrementing bottom pivots.
m : constant natural := n - nd.p;
top_piv,top_col : Bracket(1..top_pivots'last-1);
bot_piv,bot_col : Bracket(1..bot_pivots'last-1);
new_top_start : Standard_Complex_Matrices.Matrix(1..n,top_start'range(2));
new_bot_start : Standard_Complex_Matrices.Matrix(1..n,bot_start'range(2));
begin
for i in nd.children'range(1) loop
for j in nd.children'range(2) loop
if ((nd.children(i,j) /= null)
and then (nd.children(i,j).roco > 0))
then if nd.children(i,j).tp = top or nd.children(i,j).tp = mixed
then top_piv := Remove(top_pivots,nd.children(i,j).top);
top_col := Remove(top_columns,top_pivots,top_piv);
new_top_start
:= Special_Plane(n,m,k_top,top_col,top_special);
end if;
if nd.children(i,j).tp = bottom or nd.children(i,j).tp = mixed
then bot_piv := Remove(bot_pivots,nd.children(i,j).bottom);
bot_col := Remove(bot_columns,bot_pivots,bot_piv);
new_bot_start
:= Special_Plane(n,m,k_bot,bot_col,bot_special);
end if;
Solve_along_Two_Chains
(file,nd.children(i,j).all,
n,l_top+1,k_top,l_bot+1,k_bot,ind,cod,poset,
top_piv,top_col,bot_piv,bot_col,top_bs,bot_bs,
top_special,new_top_start,top_start,
bot_special,new_bot_start,bot_start,
planes,report,outlog,npaths,timings);
end if;
end loop;
end loop;
Two_General_Deform
(file,n,ind,poset,nd,top_start,top_target,bot_start,
bot_target,planes,top_bs,bot_bs,report,outlog,npaths,timings);
end Solve_along_Two_Chains_Deforming_Top_and_Bottom;
procedure Switch_Top_and_Solve_along_Two_Chains
( file : in file_type; nd : in Node;
n,l_bot,k_bot,ind : in natural; cod : in Bracket;
poset : in out Array_of_Array_of_VecMats;
bot_pivots,bot_columns : in Bracket;
bot_bs : in Bracket_System;
bot_special,bot_start,bot_target
: in Standard_Complex_Matrices.Matrix;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Assumes that l_top = k_top, l_bot < k_bot, and ind > cod'first so
-- that first a new top chain can be started, which is then solved
-- along with the existing chain for decrementing bottom pivots.
m : constant natural := n - nd.p;
new_k_top : constant natural := cod(ind);
kd : constant natural := n+1-new_k_top;
new_top_pivots : Bracket(1..m) := Complement(n,nd.top);
new_top_special : Standard_Complex_Matrices.Matrix(1..n,1..m)
:= Special_Top_Plane(m,nd.top);
new_top_target : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_top)
:= planes(ind).all;
new_top_columns : Bracket(1..m);
new_top_start : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_top);
new_top_bm : Bracket_Monomial := Maximal_Minors(n,kd);
new_top_bs : Bracket_System(0..Number_of_Brackets(new_top_bm))
:= Minor_Equations(kd,kd-nd.p,new_top_bm);
begin
for i in new_top_columns'range loop
new_top_columns(i) := i;
end loop;
new_top_start
:= Special_Plane(n,m,new_k_top,new_top_columns,new_top_special);
Solve_along_Two_Chains_Deforming_Top_and_Bottom
(file,nd,n,0,new_k_top,l_bot,k_bot,ind,cod,poset,
new_top_pivots,new_top_columns,bot_pivots,bot_columns,
new_top_bs,bot_bs,new_top_special,new_top_start,new_top_target,
bot_special,bot_start,bot_target,planes,
report,outlog,npaths,timings);
Clear(new_top_bm); Clear(new_top_bs);
end Switch_Top_and_Solve_along_Two_Chains;
procedure Switch_Top_and_Solve_along_One_Chain
( file : in file_type; nd : in Node; n,ind : in natural;
cod : in Bracket; poset : in out Array_of_Array_of_VecMats;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Assumes that l_top = k_top, l_bot = k_bot, with nd.tp = top
-- and ind > cod'first so that a new top chain can be started.
m : constant natural := n - nd.p;
new_k_top : constant natural := cod(ind);
kd : constant natural := n+1-new_k_top;
new_top_pivots : Bracket(1..m) := Complement(n,nd.top);
new_top_special : Standard_Complex_Matrices.Matrix(1..n,1..m)
:= Special_Top_Plane(m,nd.top);
new_top_target : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_top)
:= planes(ind).all;
new_top_columns : Bracket(1..m);
new_top_start : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_top);
new_top_bm : Bracket_Monomial := Maximal_Minors(n,kd);
new_top_bs : Bracket_System(0..Number_of_Brackets(new_top_bm))
:= Minor_Equations(kd,kd-nd.p,new_top_bm);
begin
for i in new_top_columns'range loop
new_top_columns(i) := i;
end loop;
new_top_start
:= Special_Plane(n,m,new_k_top,new_top_columns,new_top_special);
Solve_along_One_Chain
(file,nd,n,0,new_k_top,ind,cod(cod'first..ind),poset,
new_top_pivots,new_top_columns,new_top_bs,
new_top_special,new_top_start,new_top_target,
planes,report,outlog,npaths,timings);
Clear(new_top_bm); Clear(new_top_bs);
end Switch_Top_and_Solve_along_One_Chain;
procedure Switch_Bottom_and_Solve_along_Two_Chains
( file : in file_type; nd : in Node;
n,l_top,k_top,ind : in natural; cod : in Bracket;
poset : in out Array_of_Array_of_VecMats;
top_pivots,top_columns : in Bracket;
top_bs : in Bracket_System;
top_special,top_start,top_target
: in Standard_Complex_Matrices.Matrix;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Assumes that l_top < k_top, l_bot = k_bot, and ind > cod'first so
-- that first a new bottom chain can be launched, which is then solved
-- along with the existing chain for incrementing top pivots.
m : constant natural := n - nd.p;
new_k_bot : constant natural := cod(ind);
kd : constant natural := n+1-new_k_bot;
new_bot_pivots : Bracket(1..m) := Complement(n,nd.bottom);
new_bot_special : Standard_Complex_Matrices.Matrix(1..n,1..m)
:= Special_Bottom_Plane(m,nd.bottom);
new_bot_target : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_bot)
:= planes(ind).all;
new_bot_columns : Bracket(1..m);
new_bot_start : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_bot);
new_bot_bm : Bracket_Monomial := Maximal_Minors(n,kd);
new_bot_bs : Bracket_System(0..Number_of_Brackets(new_bot_bm))
:= Minor_Equations(kd,kd-nd.p,new_bot_bm);
begin
for i in new_bot_columns'range loop
new_bot_columns(i) := i;
end loop;
new_bot_start
:= Special_Plane(n,m,new_k_bot,new_bot_columns,new_bot_special);
Solve_along_Two_Chains_Deforming_Top_and_Bottom
(file,nd,n,l_top,k_top,0,new_k_bot,ind,cod,poset,
top_pivots,top_columns,new_bot_pivots,new_bot_columns,
top_bs,new_bot_bs,top_special,top_start,top_target,
new_bot_special,new_bot_start,new_bot_target,
planes,report,outlog,npaths,timings);
Clear(new_bot_bm); Clear(new_bot_bs);
end Switch_Bottom_and_Solve_along_Two_Chains;
procedure Switch_Bottom_and_Solve_along_One_Chain
( file : in file_type; nd : in Node; n,ind : in natural;
cod : in Bracket; poset : in out Array_of_Array_of_VecMats;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Assumes that l_top < k_top, l_bot = k_bot, and ind > cod'first so
-- that first a new bottom chain can be launched, which is then solved
-- along with the existing chain for incrementing top pivots.
m : constant natural := n - nd.p;
new_k_bot : constant natural := cod(ind);
kd : constant natural := n+1-new_k_bot;
new_bot_pivots : Bracket(1..m) := Complement(n,nd.bottom);
new_bot_special : Standard_Complex_Matrices.Matrix(1..n,1..m)
:= Special_Bottom_Plane(m,nd.bottom);
new_bot_target : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_bot)
:= planes(ind).all;
new_bot_columns : Bracket(1..m);
new_bot_start : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_bot);
new_bot_bm : Bracket_Monomial := Maximal_Minors(n,kd);
new_bot_bs : Bracket_System(0..Number_of_Brackets(new_bot_bm))
:= Minor_Equations(kd,kd-nd.p,new_bot_bm);
begin
for i in new_bot_columns'range loop
new_bot_columns(i) := i;
end loop;
new_bot_start
:= Special_Plane(n,m,new_k_bot,new_bot_columns,new_bot_special);
Solve_along_One_Chain
(file,nd,n,0,new_k_bot,ind,cod(cod'first..ind),poset,
new_bot_pivots,new_bot_columns,new_bot_bs,
new_bot_special,new_bot_start,new_bot_target,
planes,report,outlog,npaths,timings);
Clear(new_bot_bm); Clear(new_bot_bs);
end Switch_Bottom_and_Solve_along_One_Chain;
procedure Switch_Top_Bottom_and_Solve_along_Two_Chains
( file : in file_type; nd : in Node; n,ind : in natural;
cod : in Bracket; poset : in out Array_of_Array_of_VecMats;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Assumes that l_top = k_top, l_bot = k_bot, and ind > cod'first+1
-- so that first new top and bottom chains can be started which are
-- then solved along.
m : constant natural := n - nd.p;
new_k_top : constant natural := cod(ind);
new_k_bot : constant natural := cod(ind+1);
kd_top : constant natural := n+1-new_k_top;
kd_bot : constant natural := n+1-new_k_bot;
new_top_pivots : Bracket(1..m) := Complement(n,nd.top);
new_bot_pivots : Bracket(1..m) := Complement(n,nd.bottom);
new_top_special : Standard_Complex_Matrices.Matrix(1..n,1..m)
:= Special_Top_Plane(m,nd.top);
new_bot_special : Standard_Complex_Matrices.Matrix(1..n,1..m)
:= Special_Bottom_Plane(m,nd.bottom);
new_top_target : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_top)
:= planes(ind).all;
new_bot_target : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_bot)
:= planes(ind+1).all;
new_top_columns,new_bot_columns : Bracket(1..m);
new_top_start : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_top);
new_bot_start : Standard_Complex_Matrices.Matrix(1..n,1..m+1-new_k_bot);
new_top_bm : Bracket_Monomial := Maximal_Minors(n,kd_top);
new_bot_bm : Bracket_Monomial := Maximal_Minors(n,kd_bot);
new_top_bs : Bracket_System(0..Number_of_Brackets(new_top_bm))
:= Minor_Equations(kd_top,kd_top-nd.p,new_top_bm);
new_bot_bs : Bracket_System(0..Number_of_Brackets(new_bot_bm))
:= Minor_Equations(kd_bot,kd_bot-nd.p,new_bot_bm);
begin
for i in new_bot_columns'range loop
new_top_columns(i) := i;
new_bot_columns(i) := i;
end loop;
new_top_start
:= Special_Plane(n,m,new_k_top,new_top_columns,new_top_special);
new_bot_start
:= Special_Plane(n,m,new_k_bot,new_bot_columns,new_bot_special);
Solve_along_Two_Chains_Deforming_Top_and_Bottom
(file,nd,n,0,new_k_top,0,new_k_bot,ind,cod,poset,
new_top_pivots,new_top_columns,new_bot_pivots,new_bot_columns,
new_top_bs,new_bot_bs,new_top_special,new_top_start,new_top_target,
new_bot_special,new_bot_start,new_bot_target,
planes,report,outlog,npaths,timings);
Clear(new_top_bm); Clear(new_top_bs);
Clear(new_bot_bm); Clear(new_bot_bs);
end Switch_Top_Bottom_and_Solve_along_Two_Chains;
procedure Solve_along_Two_Chains
( file : in file_type; nd : in Node;
n,l_top,k_top,l_bot,k_bot,ind : in natural;
cod : in Bracket; poset : in out Array_of_Array_of_VecMats;
top_pivots,top_columns,bot_pivots,bot_columns : in Bracket;
top_bs,bot_bs : in Bracket_System;
top_special,top_start,top_target,bot_special,bot_start,
bot_target : in Standard_Complex_Matrices.Matrix;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- Does the dispatching according to the case analysis.
-- IMPORTANT :
-- The control structure in the case analysis matches the structure
-- in Localization_Posets.Recursive_Top_Bottom_Create.
begin
if empty(poset,nd.level,nd.label)
then
if nd.level = 0
then poset(nd.level)(nd.label)(1)
:= new Standard_Complex_Matrices.Matrix'(Leaf_Plane(n,nd));
elsif nd.roco > 0
then
if ((l_top < k_top) and (l_bot < k_bot))
then Solve_along_Two_Chains_Deforming_Top_and_Bottom
(file,nd,n,l_top,k_top,l_bot,k_bot,ind,cod,poset,
top_pivots,top_columns,bot_pivots,bot_columns,
top_bs,bot_bs,top_special,top_start,top_target,
bot_special,bot_start,bot_target,planes,
report,outlog,npaths,timings);
elsif ((l_top = k_top) and (l_bot < k_bot))
then
if ind = cod'first
then Solve_along_One_Chain
(file,nd,n,l_bot,k_bot,ind,cod,poset,
bot_pivots,bot_columns,bot_bs,bot_special,
bot_start,bot_target,planes,report,outlog,
npaths,timings);
else Switch_Top_and_Solve_along_Two_Chains
(file,nd,n,l_bot,k_bot,ind-1,cod,poset,
bot_pivots,bot_columns,bot_bs,bot_special,
bot_start,bot_target,planes,report,outlog,
npaths,timings);
end if;
elsif ((l_top < k_top) and (l_bot = k_bot))
then
if ind = cod'first
then Solve_along_One_Chain
(file,nd,n,l_top,k_top,ind,cod,poset,
top_pivots,top_columns,top_bs,top_special,
top_start,top_target,planes,report,outlog,
npaths,timings);
else Switch_Bottom_and_Solve_along_Two_Chains
(file,nd,n,l_top,k_top,ind-1,cod,poset,
top_pivots,top_columns,top_bs,top_special,
top_start,top_target,planes,report,outlog,
npaths,timings);
end if;
else -- ((l_top = k_top) and (l_bot = k_bot))
if ind > cod'first+1
then Switch_Top_Bottom_and_Solve_along_Two_Chains
(file,nd,n,ind-2,cod,poset,planes,
report,outlog,npaths,timings);
elsif ind > cod'first
then
if nd.tp = bottom
then
Switch_Bottom_and_Solve_along_One_Chain
(file,nd,n,ind-1,cod,poset,planes,
report,outlog,npaths,timings);
else
Switch_Top_and_Solve_along_One_Chain
(file,nd,n,ind-1,cod,poset,planes,
report,outlog,npaths,timings);
end if;
end if;
end if;
end if;
end if;
end Solve_along_Two_Chains;
procedure One_Solve
( file : in file_type; n : in natural; cod : in Bracket;
poset : in out Array_of_Array_of_VecMats;
nd : in Node; planes : in VecMat;
report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- NOTE :
-- We assume that we are always folding in the last condition, with
-- codimension cod(cod'last). In passing recursively to higher levels
-- in the deformation poset, we take a slice of k, omitting the last one.
-- Applies the solver first to all grandchildren of the current node,
-- which is the additional layer compared to the hypersurface case.
-- This implementation will only work in the non-mixed case.
m : constant natural := n - nd.p;
kk : constant natural := cod(cod'last);
kd : constant natural := n+1-kk;
procedure Solve_Grand_Child
( lnd : in Link_to_Node; continue : out boolean ) is
-- DESCRIPTION :
-- This node lnd is a grandchild of the current node.
begin
if Empty(poset,lnd.level,lnd.label)
then if lnd.level = 0
then poset(lnd.level)(lnd.label)(1)
:= new Standard_Complex_Matrices.Matrix'
(Leaf_Plane(n,lnd.all));
elsif lnd.roco > 0
then One_Solve(file,n,cod(cod'first..cod'last-1),poset,
lnd.all,planes,report,outlog,npaths,timings);
end if;
end if;
continue := true;
end Solve_Grand_Child;
procedure Solve_Grand_Children is
new Enumerate_Grand_Children(Solve_Grand_Child);
begin
if (Empty(poset,nd.level,nd.label) and (nd.roco > 0))
then
if cod'last >= cod'first
then Solve_Grand_Children(nd,kk);
end if;
declare
pivots,columns : Bracket(1..m);
special : Standard_Complex_Matrices.Matrix(1..n,1..m);
target : constant Standard_Complex_Matrices.Matrix
:= planes(cod'last).all;
start : Standard_Complex_Matrices.Matrix(1..n,1..m+1-kk);
bm : Bracket_Monomial := Maximal_Minors(n,kd);
bs : Bracket_System(0..Number_of_Brackets(bm))
:= Minor_Equations(kd,kd-nd.p,bm);
begin
for i in columns'range loop
columns(i) := i;
end loop;
if nd.tp = top
then pivots := Complement(n,nd.top);
special := Special_Top_Plane(m,nd.top);
else pivots := Complement(n,nd.bottom);
special := Special_Bottom_Plane(m,nd.bottom);
end if;
start := Special_Plane(n,m,kk,columns,special);
One_Solve_along_Chains
(file,nd,n,0,kk,cod'last,poset,pivots,columns,bs,
special,start,target,planes,report,outlog,npaths,timings);
Clear(bm); Clear(bs);
end;
end if;
end One_Solve;
procedure Chain_Solve
( file : in file_type; n : in natural; cod : in Bracket;
poset : in out Array_of_Array_of_VecMats; nd : in Node;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- NOTE :
-- The convention is that the last co-dimension condition is treated
-- when the type of the node is not mixed, otherwise the last two entries
-- of the vector of co-dimension conditions are sliced off when moving
-- to the upper levels.
-- This is another organization of One_Solve and only works when the
-- type of the nodes are not mixed.
begin
if Empty(poset,nd.level,nd.label)
then
if nd.level = 0
then poset(nd.level)(nd.label)(1)
:= new Standard_Complex_Matrices.Matrix'(Leaf_Plane(n,nd));
elsif nd.roco > 0
then declare
m : constant natural := n - nd.p;
pivots,columns : Bracket(1..m);
special : Standard_Complex_Matrices.Matrix(1..n,1..m);
kk : constant natural := cod(cod'last);
kd : constant natural := n+1-kk;
start : Standard_Complex_Matrices.Matrix(1..n,1..m+1-kk);
target : constant Standard_Complex_Matrices.Matrix
:= planes(planes'last).all;
bm : Bracket_Monomial := Maximal_Minors(n,kd);
bs : Bracket_System(0..Number_of_Brackets(bm))
:= Minor_Equations(kd,kd-nd.p,bm);
begin
for i in columns'range loop
columns(i) := i;
end loop;
if nd.tp = top
then pivots := Complement(n,nd.top);
special := Special_Top_Plane(m,nd.top);
else pivots := Complement(n,nd.bottom);
special := Special_Bottom_Plane(m,nd.bottom);
end if;
start := Special_Plane(n,m,kk,columns,special);
Solve_along_One_Chain
(file,nd,n,0,cod(cod'last),cod'last,cod,poset,
pivots,columns,bs,special,start,target,planes,
report,outlog,npaths,timings);
Clear(bm); Clear(bs);
end;
end if;
end if;
end Chain_Solve;
procedure Solve ( file : in file_type; n : in natural; cod : in Bracket;
poset : in out Array_of_Array_of_VecMats; nd : in Node;
planes : in VecMat; report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- REQUIREMENT ONE :
-- The convention is that the last co-dimension condition is treated
-- when the type of the node is not mixed, otherwise the last two entries
-- of the vector of co-dimension conditions are sliced off when moving
-- to the upper levels.
-- REQUIREMENT TWO :
-- The nodes that are not mixed appear at the top of the poset.
begin
if Empty(poset,nd.level,nd.label)
then
if nd.level = 0
then poset(nd.level)(nd.label)(1)
:= new Standard_Complex_Matrices.Matrix'(Leaf_Plane(n,nd));
elsif nd.roco > 0
then
if nd.tp /= mixed
then
declare
m : constant natural := n - nd.p;
pivots,columns : Bracket(1..m);
special : Standard_Complex_Matrices.Matrix(1..n,1..m);
kk : constant natural := cod(cod'last);
kd : constant natural := n+1-kk;
start : Standard_Complex_Matrices.Matrix(1..n,1..m+1-kk);
target : constant Standard_Complex_Matrices.Matrix
:= planes(planes'last).all;
bm : Bracket_Monomial := Maximal_Minors(n,kd);
bs : Bracket_System(0..Number_of_Brackets(bm))
:= Minor_Equations(kd,kd-nd.p,bm);
begin
for i in columns'range loop
columns(i) := i;
end loop;
if nd.tp = top
then pivots := Complement(n,nd.top);
special := Special_Top_Plane(m,nd.top);
else pivots := Complement(n,nd.bottom);
special := Special_Bottom_Plane(m,nd.bottom);
end if;
start := Special_Plane(n,m,kk,columns,special);
Solve_along_One_Chain
(file,nd,n,0,cod(cod'last),cod'last,cod,poset,
pivots,columns,bs,special,start,target,planes,
report,outlog,npaths,timings);
Clear(bm); Clear(bs);
end;
else
declare
m : constant natural := n - nd.p;
top_col,bot_col : Bracket(1..m);
kk_top : constant natural := cod(cod'last-1);
kk_bot : constant natural := cod(cod'last);
kd_top : constant natural := n+1-kk_top;
kd_bot : constant natural := n+1-kk_bot;
top_bm : Bracket_Monomial := Maximal_Minors(n,kd_top);
bot_bm : Bracket_Monomial := Maximal_Minors(n,kd_bot);
top_bs : Bracket_System(0..Number_of_Brackets(top_bm))
:= Minor_Equations(kd_top,kd_top-nd.p,top_bm);
bot_bs : Bracket_System(0..Number_of_Brackets(bot_bm))
:= Minor_Equations(kd_bot,kd_bot-nd.p,bot_bm);
top_special : Standard_Complex_Matrices.Matrix(1..n,1..m)
:= Special_Top_Plane(m,nd.top);
bot_special : Standard_Complex_Matrices.Matrix(1..n,1..m)
:= Special_Bottom_Plane(m,nd.bottom);
top_piv : Bracket(1..m) := Complement(n,nd.top);
bot_piv : Bracket(1..m) := Complement(n,nd.bottom);
top_start,top_target
: Standard_Complex_Matrices.Matrix(1..n,1..m+1-kk_top);
bot_start,bot_target
: Standard_Complex_Matrices.Matrix(1..n,1..m+1-kk_bot);
begin
for i in top_col'range loop
top_col(i) := i;
end loop;
top_start := Special_Plane(n,m,kk_top,top_col,top_special);
top_target := planes(planes'last-1).all;
for i in bot_col'range loop
bot_col(i) := i;
end loop;
bot_start := Special_Plane(n,m,kk_bot,bot_col,bot_special);
bot_target := planes(planes'last).all;
Solve_along_Two_Chains
(file,nd,n,0,kk_top,0,kk_bot,cod'last-1,cod,poset,
top_piv,top_col,bot_piv,bot_col,top_bs,bot_bs,
top_special,top_start,top_target,bot_special,bot_start,
bot_target,planes,report,outlog,npaths,timings);
Clear(top_bm); Clear(top_bs); Clear(bot_bm); Clear(bot_bs);
end;
end if;
end if;
end if;
end Solve;
procedure Recursive_Quantum_Solve
( file : in file_type; n,q : in natural;
nd : in Node; expbp : in Bracket_Polynomial;
poset : in out Array_of_Array_of_VecMats;
planes : in VecMat; s : in Standard_Complex_Vectors.Vector;
report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
-- DESCRIPTION :
-- This additional layer is added to avoid the repeated construction
-- of the structure of the equations, that is now in expbp.
-- This is the q-analogue to the Recursive Hypersurface Solver.
begin
if Empty(poset,nd.level,nd.label)
then
if nd.level = 0
then
poset(nd.level)(nd.label)(1)
:= new Standard_Complex_Matrices.Matrix'(Leaf_Plane(n*(q+1),nd));
else
for i in nd.children'range(1) loop
for j in nd.children'range(2) loop
if nd.children(i,j) /= null
then Recursive_Quantum_Solve
(file,n,q,nd.children(i,j).all,expbp,poset,
planes,s,report,outlog,npaths,timings);
end if;
end loop;
end loop;
Quantum_Deform(file,n,q,poset,nd,expbp,planes,s,report,outlog,
npaths,timings);
end if;
end if;
end Recursive_Quantum_Solve;
procedure Solve ( file : in file_type; n,q : in natural;
poset : in out Array_of_Array_of_VecMats; nd : in Node;
planes : in VecMat; s : in Standard_Complex_Vectors.Vector;
report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
bm : Bracket_Monomial := Maximal_Minors(n,n);
bs : Bracket_System(0..Number_of_Brackets(bm))
:= Minor_Equations(n,n-nd.p,bm);
begin
Recursive_Quantum_Solve(file,n,q,nd,bs(1),poset,planes,s,report,outlog,
npaths,timings);
Clear(bm); Clear(bs);
end Solve;
procedure One_Solve
( file : in file_type; n,q : in natural; cod : in Bracket;
poset : in out Array_of_Array_of_VecMats; nd : in Node;
planes : in VecMat; s : in Standard_Complex_Vectors.Vector;
report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
m : constant natural := n - nd.p;
kk : constant natural := cod(cod'last);
kd : constant natural := n+1-kk;
procedure Solve_Grand_Child
( lnd : in Link_to_Node; continue : out boolean ) is
-- DESCRIPTION :
-- This node lnd is a grandchild of the current node.
begin
if Empty(poset,lnd.level,lnd.label)
then if lnd.level = 0
then poset(lnd.level)(lnd.label)(1)
:= new Standard_Complex_Matrices.Matrix'
(Leaf_Plane(n*(q+1),lnd.all));
elsif lnd.roco > 0
then One_Solve(file,n,q,cod(cod'first..cod'last-1),poset,
lnd.all,planes,s,report,outlog,npaths,timings);
end if;
end if;
continue := true;
end Solve_Grand_Child;
procedure Solve_Grand_Children is
new Enumerate_Grand_Children(Solve_Grand_Child);
begin
if (Empty(poset,nd.level,nd.label) and (nd.roco > 0))
then
if cod'last >= cod'first
then Solve_Grand_Children(nd,kk);
end if;
declare
pivots,columns : Bracket(1..m);
mod_piv : Bracket(1..nd.p);
special : Standard_Complex_Matrices.Matrix(1..n,1..m);
target : constant Standard_Complex_Matrices.Matrix
:= planes(cod'last).all;
start : Standard_Complex_Matrices.Matrix(1..n,1..m+1-kk);
bm : Bracket_Monomial := Maximal_Minors(n,kd);
bs : Bracket_System(0..Number_of_Brackets(bm))
:= Minor_Equations(kd,kd-nd.p,bm);
begin
for i in columns'range loop
columns(i) := i;
end loop;
if nd.tp = top
then mod_piv := Modulo(nd.top,n);
pivots := Complement(n,mod_piv);
special := Special_Top_Plane(m,mod_piv);
else mod_piv := Modulo(nd.bottom,n);
pivots := Complement(n,mod_piv);
special := Special_Bottom_Plane(m,mod_piv);
end if;
start := Special_Plane(n,m,kk,columns,special);
One_Quantum_Solve_along_Chains
(file,nd,n,q,0,kk,cod'last,poset,pivots,columns,bs,
special,start,target,planes,s,report,outlog,npaths,timings);
Clear(bm); Clear(bs);
end;
end if;
end One_Solve;
procedure Solve ( file : in file_type; n,q : in natural; cod : in Bracket;
poset : in out Array_of_Array_of_VecMats; nd : in Node;
planes : in VecMat; s : in Standard_Complex_Vectors.Vector;
report,outlog : in boolean;
npaths : in out Standard_Natural_Vectors.Vector;
timings : in out Duration_Array ) is
begin
null;
end Solve;
-- DESTRUCTORS :
procedure Clear ( avm : in out Array_of_VecMats ) is
begin
for i in avm'range loop
Deep_Clear(avm(i));
end loop;
end Clear;
procedure Clear ( avm : in out Link_to_Array_of_VecMats ) is
procedure free is
new unchecked_deallocation(Array_of_VecMats,Link_to_Array_of_VecMats);
begin
if avm /= null
then Clear(avm.all);
free(avm);
end if;
end Clear;
procedure Clear ( avm : in out Array_of_Array_of_VecMats ) is
begin
for i in avm'range loop
if avm(i) /= null
then Clear(avm(i).all);
end if;
end loop;
end Clear;
end Deformation_Posets;