File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices / ts_cmpmat.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:24 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 text_io,integer_io; use text_io,integer_io;
with Standard_Natural_Vectors;
with Standard_Floating_Numbers; use Standard_Floating_Numbers;
with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
with Standard_Complex_Numbers;
with Standard_Complex_Vectors;
with Standard_Complex_Vectors_io;
with Standard_Complex_Matrices;
with Standard_Complex_Matrices_io;
with Standard_Complex_VecMats;
with Standard_Complex_VecMats_io;
with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
with Standard_Random_Vectors; use Standard_Random_Vectors;
with Multprec_Random_Vectors; use Multprec_Random_Vectors;
with Multprec_Floating_Numbers; use Multprec_Floating_Numbers;
with Multprec_Floating_Numbers_io; use Multprec_Floating_Numbers_io;
with Multprec_Complex_Numbers;
with Multprec_Complex_Vectors;
with Multprec_Complex_Vectors_io;
with Multprec_Complex_Matrices;
with Multprec_Complex_Matrices_io;
with Multprec_Complex_Linear_Solvers; use Multprec_Complex_Linear_Solvers;
with Multprec_Complex_Norms_Equals; use Multprec_Complex_Norms_Equals;
procedure ts_cmpmat is
-- DESCRIPTION :
-- Tests the matrix packages of standard and multi-precision complex numbers.
procedure Test_Standard_io is
use Standard_Complex_Matrices,Standard_Complex_Matrices_io;
n,m : natural;
begin
put("Give the number of rows : "); get(n);
put("Give the number of columns : "); get(m);
declare
mat : Matrix(1..n,1..m);
begin
put("Give "); put(n,1); put("x"); put(m,1);
put_line(" complex matrix : "); get(mat);
put_line("Your matrix : "); put(mat); new_line;
end;
end Test_Standard_io;
procedure Test_Standard_VecMat_io is
use Standard_Complex_Matrices,Standard_Complex_Matrices_io;
use Standard_Complex_VecMats,Standard_Complex_VecMats_io;
n,n1,n2 : natural;
lv : Link_to_VecMat;
begin
put("Give the number of matrices : "); get(n);
put("Give #rows : "); get(n1);
put("Give #columns : "); get(n2);
put("Give "); put(n,1); put(" "); put(n1,1); put("-by-"); put(n2,1);
put_line(" integer matrices : ");
get(n,n1,n2,lv);
put_line("The vector of matrices :"); put(lv);
end Test_Standard_VecMat_io;
procedure Test_Multprec_io is
use Multprec_Complex_Matrices,Multprec_Complex_Matrices_io;
n,m : natural;
begin
put("Give the number of rows : "); get(n);
put("Give the number of columns : "); get(m);
declare
mat : Matrix(1..n,1..m);
begin
put("Give "); put(n,1); put("x"); put(m,1);
put_line(" complex matrix : "); get(mat);
put_line("Your matrix : "); put(mat); new_line;
end;
end Test_Multprec_io;
function Vdm_Matrix ( v : Standard_Complex_Vectors.Vector )
return Standard_Complex_Matrices.Matrix is
use Standard_Complex_Numbers;
use Standard_Complex_Matrices;
n : constant natural := v'length;
res : Matrix(1..n,1..n);
begin
for i in res'range(1) loop
for j in res'range(2) loop
res(i,j) := v(i)**(j-1);
end loop;
end loop;
return res;
end Vdm_Matrix;
procedure lufac_Solve ( n : in natural;
mat : in Standard_Complex_Matrices.Matrix;
rhs : in Standard_Complex_Vectors.Vector ) is
use Standard_Complex_Vectors;
use Standard_Complex_Vectors_io;
use Standard_Complex_Matrices;
use Standard_Complex_Matrices_io;
wrk : Matrix(mat'range(1),mat'range(2));
piv : Standard_Natural_Vectors.Vector(mat'range(2));
info : natural;
res,sol,acc : Vector(rhs'range);
nrm : double_float;
begin
put_line("Solving the linear system with lufac.");
wrk := mat;
lufac(wrk,n,piv,info);
put("info : "); put(info,1); new_line;
sol := rhs;
lusolve(wrk,n,piv,sol);
put_line("The solution vector :"); put(sol); new_line;
acc := mat*sol;
res := rhs - acc;
put_line("The residual : "); put(res); new_line;
nrm := Max_Norm(res);
put("Max norm of residual : "); put(nrm); new_line;
nrm := Sum_Norm(res);
put("Sum norm of residual : "); put(nrm); new_line;
end lufac_Solve;
procedure lufco_Solve ( n : in natural;
mat : in Standard_Complex_Matrices.Matrix;
rhs : in Standard_Complex_Vectors.Vector ) is
use Standard_Complex_Vectors;
use Standard_Complex_Vectors_io;
use Standard_Complex_Matrices;
use Standard_Complex_Matrices_io;
wrk : Matrix(mat'range(1),mat'range(2)) := mat;
piv : Standard_Natural_Vectors.Vector(mat'range(2));
rcond,nrm : double_float;
res,sol : Vector(rhs'range);
begin
put_line("Solving the linear system with lufco.");
lufco(wrk,n,piv,rcond);
put("inverse condition : "); put(rcond); new_line;
sol := rhs;
lusolve(wrk,n,piv,sol);
put_line("The solution vector :"); put(sol); new_line;
res := rhs - mat*sol;
put_line("The residual : "); put(res); new_line;
nrm := Max_Norm(res);
put("Max norm of residual : "); put(nrm); new_line;
nrm := Sum_Norm(res);
put("Sum norm of residual : "); put(nrm); new_line;
end lufco_Solve;
procedure Random_Test_Standard_Linear_Solvers is
use Standard_Complex_Vectors;
use Standard_Complex_Vectors_io;
use Standard_Complex_Matrices;
use Standard_Complex_Matrices_io;
n,nb : natural;
begin
new_line;
put_line("Testing of solving random standard linear systems.");
new_line;
put("Give the dimension : "); get(n);
put("Give the number of tests : "); get(nb);
for i in 1..nb loop
declare
mat : Matrix(1..n,1..n);
rhs : Vector(1..n);
begin
mat := Vdm_Matrix(Random_Vector(1,n));
rhs := Random_Vector(1,n);
lufac_Solve(n,mat,rhs);
lufco_Solve(n,mat,rhs);
end;
end loop;
end Random_Test_Standard_Linear_Solvers;
function Vdm_Matrix ( v : Multprec_Complex_Vectors.Vector )
return Multprec_Complex_Matrices.Matrix is
use Multprec_Complex_Numbers;
use Multprec_Complex_Matrices;
n : constant natural := v'length;
res : Matrix(1..n,1..n);
begin
for i in res'range(1) loop
for j in res'range(2) loop
res(i,j) := v(i)**(j-1);
end loop;
end loop;
return res;
end Vdm_Matrix;
procedure lufac_Solve ( n : in natural;
mat : in Multprec_Complex_Matrices.Matrix;
rhs : in Multprec_Complex_Vectors.Vector ) is
use Multprec_Complex_Vectors;
use Multprec_Complex_Vectors_io;
use Multprec_Complex_Matrices;
use Multprec_Complex_Matrices_io;
wrk : Matrix(mat'range(1),mat'range(2)) := +mat;
piv : Standard_Natural_Vectors.Vector(mat'range(2));
info : natural;
res,sol,acc : Vector(rhs'range);
nrm : Floating_Number;
begin
put_line("Solving the linear system with lufac.");
lufac(wrk,n,piv,info);
put("info : "); put(info,1); new_line;
sol := +rhs;
lusolve(wrk,n,piv,sol);
put_line("The solution vector :"); put(sol); new_line;
acc := mat*sol;
res := rhs - acc;
put_line("The residual : "); put(res); new_line;
nrm := Max_Norm(res);
put("Max norm of residual : "); put(nrm); new_line;
Clear(nrm);
nrm := Sum_Norm(res);
put("Sum norm of residual : "); put(nrm); new_line;
Clear(nrm);
Clear(wrk); Clear(acc); Clear(res); Clear(sol);
end lufac_Solve;
procedure lufco_Solve ( n : in natural;
mat : in Multprec_Complex_Matrices.Matrix;
rhs : in Multprec_Complex_Vectors.Vector ) is
use Multprec_Complex_Vectors;
use Multprec_Complex_Vectors_io;
use Multprec_Complex_Matrices;
use Multprec_Complex_Matrices_io;
wrk : Matrix(mat'range(1),mat'range(2)) := +mat;
piv : Standard_Natural_Vectors.Vector(mat'range(2));
rcond,nrm : Floating_Number;
res,sol,acc : Vector(rhs'range);
begin
put_line("Solving the linear system with lufco.");
lufco(wrk,n,piv,rcond);
put("inverse condition : "); put(rcond); new_line;
sol := +rhs;
lusolve(wrk,n,piv,sol);
put_line("The solution vector :"); put(sol); new_line;
acc := mat*sol;
res := rhs - acc;
put_line("The residual : "); put(res); new_line;
nrm := Max_Norm(res);
put("Max norm of residual : "); put(nrm); new_line;
Clear(nrm);
nrm := Sum_Norm(res);
put("Sum norm of residual : "); put(nrm); new_line;
Clear(wrk); Clear(acc); Clear(res); Clear(sol);
Clear(nrm); Clear(rcond);
end lufco_Solve;
procedure Random_Test_Multprec_Linear_Solvers is
use Multprec_Complex_Vectors;
use Multprec_Complex_Vectors_io;
use Multprec_Complex_Matrices;
use Multprec_Complex_Matrices_io;
n,sz,nb : natural;
begin
new_line;
put_line("Testing of solving random multi-precision linear systems.");
new_line;
put("Give the dimension : "); get(n);
put("Give the size of the numbers : "); get(sz);
put("Give the number of tests : "); get(nb);
for i in 1..nb loop
declare
rnd : Vector(1..n) := Random_Vector(1,n,sz);
mat : Matrix(1..n,1..n) := Vdm_Matrix(rnd);
rhs : Vector(1..n) := Random_Vector(1,n,sz);
begin
lufac_Solve(n,mat,rhs);
-- lufco_Solve(n,mat,rhs);
Clear(rnd); Clear(mat); Clear(rhs);
end;
end loop;
end Random_Test_Multprec_Linear_Solvers;
procedure Main is
ans : character;
begin
new_line;
put_line("Interactive testing of matrices of complex numbers");
loop
new_line;
put_line("Choose one of the following : ");
put_line(" 0. exit this program. ");
put_line(" 1. io of matrices of standard complex numbers. ");
put_line(" 2. io of vectors of matrices of standard complex numbers. ");
put_line(" 3. test on solving random standard linear systems. ");
put_line(" 4. io of matrices of multi-precision complex numbers. ");
put_line(" 5. test on solving random multi-precision complex numbers.");
put("Make your choice (0,1,2,3,4 or 5) : "); get(ans);
exit when (ans = '0');
case ans is
when '1' => Test_Standard_io;
when '2' => Test_Standard_VecMat_io;
when '3' => Random_Test_Standard_Linear_Solvers;
when '4' => Test_Multprec_io;
when '5' => Random_Test_Multprec_Linear_Solvers;
when others => null;
end case;
end loop;
end Main;
begin
Main;
end ts_cmpmat;