Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Supports/ts_givrot.adb, Revision 1.1.1.1
1.1 maekawa 1: with text_io,integer_io; use text_io,integer_io;
2: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
3: with Standard_Integer_Vectors;
4: with Standard_Integer_Vectors_io; use Standard_Integer_Vectors_io;
5: with Standard_Floating_Vectors; use Standard_Floating_Vectors;
6: with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
7: with Standard_Floating_Matrices; use Standard_Floating_Matrices;
8: with Standard_Floating_Matrices_io; use Standard_Floating_Matrices_io;
9: with Givens_Rotations; use Givens_Rotations;
10:
11: procedure ts_givrot is
12:
13: -- DESCRIPTION :
14: -- Interactive testing of Givens rotations.
15:
16: n,m : natural;
17: ans : character;
18: tol : constant double_float := 10.0**(-12);
19:
20: function Compute_Residual ( mat : in matrix; rhs,x : in vector;
21: ipvt : in Standard_Integer_Vectors.Vector )
22: return Vector is
23:
24: res : Vector(rhs'range) := rhs;
25:
26: begin
27: for i in mat'range(1) loop
28: for j in mat'range(2) loop
29: res(i) := res(i) - mat(i,ipvt(j))*x(j);
30: exit when j > mat'last(2);
31: end loop;
32: end loop;
33: return res;
34: end Compute_Residual;
35:
36: procedure Row_Wise_Triangulation ( mat : in out Matrix ) is
37:
38: pivots : Standard_Integer_Vectors.Vector(mat'range(2));
39: ipiv : integer;
40:
41: begin
42: for i in pivots'range loop
43: pivots(i) := i;
44: end loop;
45: for i in mat'range(1) loop
46: Upper_Triangulate(i,mat,tol,pivots,ipiv);
47: put("At row "); put(i,1); put(" we have pivot : "); put(ipiv,1);
48: new_line;
49: put_line("with matrix : "); put(mat);
50: put("and pivot vector "); put(pivots); new_line;
51: exit when (ipiv = 0);
52: end loop;
53: end Row_Wise_Triangulation;
54:
55: begin
56: loop
57: put("Give the number of rows : "); get(n);
58: put("Give the number of columns : "); get(m);
59: declare
60: mat,wrkmat : matrix(1..n,1..m);
61: rhs,wrkrhs : vector(1..n);
62: ipvt : Standard_Integer_Vectors.Vector(1..m);
63: begin
64: put("Give the elements of the ");
65: put(n,1); put("x"); put(m,1); put_line("-matrix :"); get(mat);
66: put_line("The matrix :"); put(mat);
67: wrkmat := mat;
68: put("Do you have a right hand side ? (y/n) "); get(ans);
69: if ans = 'y'
70: then put("Give "); put(n,1); put(" floating point numbers : ");
71: rhs := (1..n => 0.0); get(rhs);
72: wrkrhs := rhs;
73: Upper_Triangulate(wrkmat,wrkrhs,tol,ipvt);
74: else Row_Wise_Triangulation(wrkmat);
75: wrkmat := mat;
76: Upper_Triangulate(wrkmat,tol,ipvt);
77: end if;
78: put_line("The matrix in upper triangular form : "); put(wrkmat);
79: put(" with pivoting information : "); put(ipvt); new_line;
80: if ans = 'y'
81: then put(" and transformed right hand side : ");
82: put(wrkrhs); new_line;
83: declare
84: sol,res : Vector(1..n);
85: begin
86: Solve(wrkmat,wrkrhs,tol,sol);
87: put("The solution : "); put(sol); new_line;
88: res := Compute_Residual(mat,rhs,sol,ipvt);
89: put(" with residual : "); put(res); new_line;
90: end;
91: end if;
92: end;
93: put("Do you want more tests ? (y/n) "); get(ans);
94: exit when ans /= 'y';
95: end loop;
96: end ts_givrot;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>