File: [local] / OpenXM_contrib / PHC / Ada / Schubert / drivers_for_input_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 integer_io,Numbers_io; use integer_io,Numbers_io;
with Communications_with_User; use Communications_with_User;
with Standard_Floating_Numbers; use Standard_Floating_Numbers;
with Standard_Complex_Numbers; use Standard_Complex_Numbers;
with Standard_Random_Numbers; use Standard_Random_Numbers;
with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
with Standard_Floating_Matrices;
with Standard_Floating_Matrices_io; use Standard_Floating_Matrices_io;
with Standard_Complex_Matrices;
with Standard_Random_Matrices; use Standard_Random_Matrices;
with Standard_Complex_VecMats_io; use Standard_Complex_VecMats_io;
with Osculating_Planes; use Osculating_Planes;
package body Drivers_for_Input_Planes is
function Finite ( dim : Bracket; fin_sum : natural ) return boolean is
-- DESCRIPTION :
-- Returns true if the co-dimensions in dim will lead to a finite
-- number of solutions.
sum : natural := 0;
begin
for i in dim'range loop
sum := sum + dim(i);
end loop;
if sum = fin_sum
then return true;
else return false;
end if;
end Finite;
function Read_Codimensions ( m,p,q : natural ) return Bracket is
mpq : constant natural := m*p + q*(m+p);
codim : Bracket(1..mpq);
n : natural;
begin
loop
put("Give number of intersection conditions : "); get(n);
put("Give "); put(n,1); put(" co-dimensions (sum = ");
put(mpq,1); put(") : ");
for i in 1..n loop
get(codim(i));
end loop;
for i in 1..n-1 loop
put(codim(i),1); put(" + ");
end loop;
put(codim(n),1);
if Finite(codim(1..n),mpq)
then put(" = "); put(mpq,1); put_line(" Finite #sols.");
exit;
else put(" /= "); put(mpq,1);
put_line(" Please try again.");
end if;
end loop;
return codim(1..n);
end Read_Codimensions;
function Convert ( mat : Standard_Floating_Matrices.Matrix )
return Standard_Complex_Matrices.Matrix is
-- DESCRIPTION :
-- Converts a real matrix into a complex one.
res : Standard_Complex_Matrices.Matrix(mat'range(1),mat'range(2));
begin
for i in mat'range(1) loop
for j in mat'range(2) loop
res(i,j) := Create(mat(i,j));
end loop;
end loop;
return res;
end Convert;
function Random_Complex_Planes ( m,p : natural ) return VecMat is
res : VecMat(1..m*p);
n : constant natural := m+p;
use Standard_Complex_Matrices;
begin
for i in res'range loop
res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m));
end loop;
return res;
end Random_Complex_Planes;
function Random_Complex_Planes ( m,p : natural; k : Bracket ) return VecMat is
res : VecMat(k'range);
n : constant natural := m+p;
use Standard_Complex_Matrices;
begin
for i in res'range loop
res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m+1-k(i)));
end loop;
return res;
end Random_Complex_Planes;
function Random_Complex_Planes ( m,p,q : natural ) return VecMat is
res : VecMat(1..(m*p+q*(m+p)));
n : constant natural := m+p;
use Standard_Complex_Matrices;
begin
for i in res'range loop
res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m));
end loop;
return res;
end Random_Complex_Planes;
function Random_Real_Planes ( m,p : natural ) return VecMat is
res : VecMat(1..m*p);
n : constant natural := m+p;
fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m);
begin
for i in res'range loop
fltmat := Random_Orthogonal_Matrix(n,m);
cmpmat := Convert(fltmat);
res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
end loop;
return res;
end Random_Real_Planes;
function Random_Real_Planes ( m,p : natural; k : Bracket ) return VecMat is
res : VecMat(k'range);
n : constant natural := m+p;
begin
for i in res'range loop
declare
fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m+1-k(i));
begin
fltmat := Random_Orthogonal_Matrix(n,m+1-k(i));
cmpmat := Convert(fltmat);
res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
end;
end loop;
return res;
end Random_Real_Planes;
function Random_Real_Planes ( m,p,q : natural ) return VecMat is
res : VecMat(1..(m*p+q*(m+p)));
n : constant natural := m+p;
fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m);
begin
for i in res'range loop
fltmat := Random_Orthogonal_Matrix(n,m);
cmpmat := Convert(fltmat);
res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
end loop;
return res;
end Random_Real_Planes;
function Equidistant_Interpolation_Points ( n : natural ) return Vector is
res : Vector(1..n);
s : double_float := Random;
inc : constant double_float := 2.0/double_float(n);
begin
for i in res'range loop
res(i) := s;
s := s+inc;
if s >= 1.0
then s := s - 2.0;
end if;
end loop;
return res;
end Equidistant_Interpolation_Points;
function Read_Interpolation_Points ( n : natural ) return Vector is
res : Vector(1..n);
begin
new_line;
put("Give "); put(n,1); put_line(" distinct real interpolation points : ");
for i in res'range loop
Read_Double_Float(res(i));
end loop;
return res;
end Read_Interpolation_Points;
function Osculating_Input_Planes ( m,p : natural ) return VecMat is
s : constant Vector := Equidistant_Interpolation_Points(m*p);
begin
return Osculating_Input_Planes(m,p,s);
end Osculating_Input_Planes;
function Osculating_Input_Planes
( m,p : natural; s : Vector ) return VecMat is
res : VecMat(1..m*p);
n : constant natural := m+p;
realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
begin
for i in res'range loop
realmat := Orthogonal_Basis(n,m,s(i));
res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
end loop;
return res;
end Osculating_Input_Planes;
function Osculating_Input_Planes
( m,p : natural; k : bracket ) return VecMat is
s : constant Vector := Equidistant_Interpolation_Points(k'length);
begin
return Osculating_Input_Planes(m,p,k,s);
end Osculating_Input_Planes;
function Osculating_Input_Planes
( m,p : natural; k : bracket; s : Vector ) return VecMat is
res : VecMat(k'range);
n : constant natural := m+p;
begin
for i in res'range loop
declare
realmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
begin
realmat := Orthogonal_Basis(n,m+1-k(i),s(i));
res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
end;
end loop;
return res;
end Osculating_Input_Planes;
function Osculating_Input_Planes ( m,p,q : natural ) return VecMat is
dim : constant natural := m*p + q*(m+p);
s : constant Vector := Equidistant_Interpolation_Points(dim);
begin
return Osculating_Input_Planes(m,p,q,s);
end Osculating_Input_Planes;
function Osculating_Input_Planes
( m,p,q : natural; s : Vector ) return VecMat is
res : VecMat(1..m*p+q*(m+p));
n : constant natural := m+p;
realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
begin
for i in res'range loop
realmat := Orthogonal_Basis(n,m,s(i));
res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
end loop;
return res;
end Osculating_Input_Planes;
function Read_Input_Planes ( m,p : natural ) return VecMat is
res : VecMat(1..m*p);
planesfile : file_type;
n : constant natural := m+p;
realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
begin
new_line;
put_line("Reading the name of the file with the input planes.");
Read_Name_and_Open_File(planesfile);
for i in res'range loop
get(planesfile,realmat);
realmat := Orthogonalize(realmat);
res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
end loop;
Close(planesfile);
return res;
end Read_Input_Planes;
function Read_Input_Planes ( m,p : natural; k : in Bracket ) return VecMat is
res : VecMat(k'range);
planesfile : file_type;
n : constant natural := m+p;
begin
new_line;
put_line("Reading the name of the file with the input planes.");
Read_Name_and_Open_File(planesfile);
for i in res'range loop
declare
realmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
begin
get(planesfile,realmat);
realmat := Orthogonalize(realmat);
res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
end;
end loop;
Close(planesfile);
return res;
end Read_Input_Planes;
function Read_Input_Planes ( m,p,q : natural ) return VecMat is
res : VecMat(1..m*p+q*(m+p));
planesfile : file_type;
n : constant natural := m+p;
realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
begin
new_line;
put_line("Reading the name of the file with the input planes.");
Read_Name_and_Open_File(planesfile);
for i in res'range loop
get(planesfile,realmat);
realmat := Orthogonalize(realmat);
res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
end loop;
Close(planesfile);
return res;
end Read_Input_Planes;
-- MAIN INTERACTIVE DRIVERS :
function Select_Input_Choice return character is
-- DESCRIPTION :
-- Displays the menu to obtain a choice for the kind of input.
res : character;
begin
new_line;
put_line("MENU for generating real input planes : ");
put_line(" 0. Random real planes at equidistant interpolation points.");
put_line(" 1. Generate input planes osculating a rational normal curve.");
put_line(" 2. Interactively give s-values for the "
& "osculating input planes.");
put_line(" 3. Give the name of the file with input planes.");
put("Type 0, 1, 2, or 3 to select : "); Ask_Alternative(res,"0123");
return res;
end Select_Input_Choice;
procedure Driver_for_Input_Planes
( file : in file_type; m,p : in natural; planes : out VecMat ) is
input_choice : character := Select_Input_Choice;
dim : constant natural := m*p;
svals : Vector(1..dim);
begin
case input_choice is
when '0' => planes := Random_Real_Planes(m,p);
when '1' => planes := Osculating_Input_Planes(m,p);
when '2' => svals := Read_Interpolation_Points(dim);
planes := Osculating_Input_Planes(m,p,svals);
when '3' => planes := Read_Input_Planes(m,p);
when others => null;
end case;
end Driver_for_Input_Planes;
procedure Driver_for_Input_Planes
( file : in file_type; m,p : in natural; k : in Bracket;
planes : out VecMat ) is
input_choice : character := Select_Input_Choice;
svals : Vector(k'range);
begin
case input_choice is
when '0' => planes := Random_Real_Planes(m,p,k);
when '1' => planes := Osculating_Input_Planes(m,p,k);
when '2' => svals := Read_Interpolation_Points(k'length);
planes := Osculating_Input_Planes(m,p,k,svals);
when '3' => planes := Read_Input_Planes(m,p,k);
when others => null;
end case;
end Driver_for_Input_Planes;
procedure Driver_for_Input_Planes
( file : in file_type; m,p,q : in natural;
s : out Vector; planes : out VecMat ) is
input_choice : character := Select_Input_Choice;
dim : constant natural := m*p + q*(m+p);
svals : Vector(1..dim);
begin
case input_choice is
when '0' => svals := Equidistant_Interpolation_Points(dim);
s := svals;
planes := Random_Real_Planes(m,p,q);
when '1' => svals := Equidistant_Interpolation_Points(dim);
s := svals;
planes := Osculating_Input_Planes(m,p,q,svals);
when '2' => svals := Read_Interpolation_Points(dim);
s := svals;
planes := Osculating_Input_Planes(m,p,q,svals);
when '3' => svals := Read_Interpolation_Points(dim);
s := svals;
planes := Read_Input_Planes(m,p,q);
when others => null;
end case;
put_line(file,"The interpolation points : ");
put_line(file,svals);
put_line(file,"The target planes : ");
put(file,planes,3);
end Driver_for_Input_Planes;
end Drivers_for_Input_Planes;