File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports / ts_givrot.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:28 2000 UTC (23 years, 10 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
with text_io,integer_io; use text_io,integer_io;
with Standard_Floating_Numbers; use Standard_Floating_Numbers;
with Standard_Integer_Vectors;
with Standard_Integer_Vectors_io; use Standard_Integer_Vectors_io;
with Standard_Floating_Vectors; use Standard_Floating_Vectors;
with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
with Standard_Floating_Matrices; use Standard_Floating_Matrices;
with Standard_Floating_Matrices_io; use Standard_Floating_Matrices_io;
with Givens_Rotations; use Givens_Rotations;
procedure ts_givrot is
-- DESCRIPTION :
-- Interactive testing of Givens rotations.
n,m : natural;
ans : character;
tol : constant double_float := 10.0**(-12);
function Compute_Residual ( mat : in matrix; rhs,x : in vector;
ipvt : in Standard_Integer_Vectors.Vector )
return Vector is
res : Vector(rhs'range) := rhs;
begin
for i in mat'range(1) loop
for j in mat'range(2) loop
res(i) := res(i) - mat(i,ipvt(j))*x(j);
exit when j > mat'last(2);
end loop;
end loop;
return res;
end Compute_Residual;
procedure Row_Wise_Triangulation ( mat : in out Matrix ) is
pivots : Standard_Integer_Vectors.Vector(mat'range(2));
ipiv : integer;
begin
for i in pivots'range loop
pivots(i) := i;
end loop;
for i in mat'range(1) loop
Upper_Triangulate(i,mat,tol,pivots,ipiv);
put("At row "); put(i,1); put(" we have pivot : "); put(ipiv,1);
new_line;
put_line("with matrix : "); put(mat);
put("and pivot vector "); put(pivots); new_line;
exit when (ipiv = 0);
end loop;
end Row_Wise_Triangulation;
begin
loop
put("Give the number of rows : "); get(n);
put("Give the number of columns : "); get(m);
declare
mat,wrkmat : matrix(1..n,1..m);
rhs,wrkrhs : vector(1..n);
ipvt : Standard_Integer_Vectors.Vector(1..m);
begin
put("Give the elements of the ");
put(n,1); put("x"); put(m,1); put_line("-matrix :"); get(mat);
put_line("The matrix :"); put(mat);
wrkmat := mat;
put("Do you have a right hand side ? (y/n) "); get(ans);
if ans = 'y'
then put("Give "); put(n,1); put(" floating point numbers : ");
rhs := (1..n => 0.0); get(rhs);
wrkrhs := rhs;
Upper_Triangulate(wrkmat,wrkrhs,tol,ipvt);
else Row_Wise_Triangulation(wrkmat);
wrkmat := mat;
Upper_Triangulate(wrkmat,tol,ipvt);
end if;
put_line("The matrix in upper triangular form : "); put(wrkmat);
put(" with pivoting information : "); put(ipvt); new_line;
if ans = 'y'
then put(" and transformed right hand side : ");
put(wrkrhs); new_line;
declare
sol,res : Vector(1..n);
begin
Solve(wrkmat,wrkrhs,tol,sol);
put("The solution : "); put(sol); new_line;
res := Compute_Residual(mat,rhs,sol,ipvt);
put(" with residual : "); put(res); new_line;
end;
end if;
end;
put("Do you want more tests ? (y/n) "); get(ans);
exit when ans /= 'y';
end loop;
end ts_givrot;