Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/ts_fltdls.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_Floating_Numbers_io; use Standard_Floating_Numbers_io;
4: with Standard_Natural_Vectors; use Standard_Natural_Vectors;
5: with Standard_Natural_Vectors_io; use Standard_Natural_Vectors_io;
6: with Standard_Floating_Vectors; use Standard_Floating_Vectors;
7: with Standard_Floating_Matrices; use Standard_Floating_Matrices;
8: with Standard_Floating_Matrices_io; use Standard_Floating_Matrices_io;
9: with Standard_Floating_Linear_Solvers; use Standard_Floating_LInear_Solvers;
10: with Standard_Random_Matrices; use Standard_Random_Matrices;
11:
12: procedure ts_fltdls is
13:
14: -- DESCRIPTION :
15: -- Test the dynamic triangulators of floating-point matrices.
16:
17: tol : constant double_float := 10.0**(-8);
18:
19: procedure Write_Matrix ( mat : in Matrix ) is
20: begin
21: for i in mat'range(1) loop
22: for j in mat'range(2) loop
23: put(mat(i,j),3,3,3);
24: end loop;
25: new_line;
26: end loop;
27: end Write_Matrix;
28:
29: procedure Write_Vector ( v : in Standard_Floating_Vectors.Vector ) is
30: begin
31: for i in v'range loop
32: put(v(i),3,3,3);
33: end loop;
34: new_line;
35: end Write_Vector;
36:
37: function Product ( n,col : natural; mat : Matrix;
38: sol : Standard_Floating_Vectors.Vector )
39: return Standard_Floating_Vectors.Vector is
40:
41: -- DESCRIPTION :
42: -- Returns the product of the solution with the columns of the matrix.
43:
44: res : Standard_Floating_Vectors.Vector(1..n);
45:
46: begin
47: for i in 1..n loop
48: res(i) := mat(i,col)*sol(n+1);
49: for j in 1..n loop
50: res(i) := res(i) + mat(i,j)*sol(j);
51: end loop;
52: end loop;
53: return res;
54: end Product;
55:
56: function Residual ( n,col : natural; mat : Matrix;
57: ipvt : Standard_Natural_Vectors.Vector;
58: sol : Standard_Floating_Vectors.Vector )
59: return Standard_Floating_Vectors.Vector is
60:
61: -- DESCRIPTION :
62: -- Returns the residual of the solution.
63:
64: res : Standard_Floating_Vectors.Vector(1..n);
65:
66: begin
67: for i in 1..n loop
68: res(i) := mat(i,ipvt(col))*sol(n+1);
69: for j in 1..n loop
70: res(i) := res(i) + mat(i,ipvt(j))*sol(j);
71: end loop;
72: end loop;
73: return res;
74: end Residual;
75:
76: procedure Upper_Triangulate ( n,m : in natural; mat : in Matrix ) is
77:
78: -- DESCRIPTION :
79: -- Triangulates the matrix, computes a solution and checks residual.
80:
81: wrk : Matrix(1..n,1..m) := mat;
82: ipvt : Standard_Natural_Vectors.Vector(1..m);
83: pivot : integer;
84: sol : Standard_Floating_Vectors.Vector(1..n+1);
85: res,prod : Standard_Floating_Vectors.Vector(1..n);
86:
87: begin
88: put_line("The matrix on entry : "); Write_Matrix(mat);
89: for i in 1..m loop
90: ipvt(i) := i;
91: end loop;
92: for i in 1..n loop
93: Upper_Triangulate(i,wrk,tol,ipvt,pivot);
94: exit when (pivot = 0);
95: end loop;
96: put_line("The triangulated matrix : "); Write_Matrix(wrk);
97: put("Pivots : "); put(ipvt); new_line;
98: for i in n+1..m loop
99: sol := Solve(n,i,wrk);
100: put("Solution : "); Write_Vector(sol);
101: prod := Product(n,i,wrk,sol);
102: put("Product : "); Write_Vector(prod);
103: res := Residual(n,i,mat,ipvt,sol);
104: put("Residual : "); Write_Vector(res);
105: end loop;
106: end Upper_Triangulate;
107:
108: procedure Interactive_Testing is
109:
110: n,m : natural;
111:
112: begin
113: put("Give number of rows : "); get(n);
114: put("Give number of columns : "); get(m);
115: put("Give "); put(n,1); put("x"); put(m,1); put_line(" matrix : ");
116: declare
117: mat : Matrix(1..n,1..m);
118: begin
119: get(mat);
120: Upper_Triangulate(n,m,mat);
121: end;
122: end Interactive_Testing;
123:
124: procedure Random_Testing is
125:
126: n,m,nb : natural;
127:
128: begin
129: put("Give number of rows : "); get(n);
130: put("Give number of columns : "); get(m);
131: put("Give number of tests : "); get(nb);
132: for i in 1..nb loop
133: declare
134: mat : Matrix(1..n,1..m) := Random_Matrix(n,m);
135: begin
136: Upper_Triangulate(n,m,mat);
137: end;
138: end loop;
139: end Random_Testing;
140:
141: procedure Main is
142:
143: ans : character;
144:
145: begin
146: new_line;
147: put_line("Testing the dynamic triangulators of floating-point matrices.");
148: loop
149: new_line;
150: put_line("Choose one of the following : ");
151: put_line(" 0. Exit this program. ");
152: put_line(" 1. Interactive testing of dynamic triangulators. ");
153: put_line(" 2. Random testing of dynamic triangulators. ");
154: put("Type 0, 1, or 2 to select : "); get(ans);
155: exit when (ans = '0');
156: case ans is
157: when '1' => Interactive_Testing;
158: when '2' => Random_Testing;
159: when others => null;
160: end case;
161: end loop;
162: end Main;
163:
164: begin
165: Main;
166: end ts_fltdls;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>