![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to minor_computations.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Schubert|
Return to minor_computations.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Schubert |
1.1 maekawa 1: with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
2: with Standard_Natural_Vectors;
3: with Standard_Natural_Vectors_io; use Standard_Natural_Vectors_io;
4: with Standard_Floating_Linear_Solvers; use Standard_Floating_Linear_Solvers;
5:
6: package body Minor_Computations is
7:
8: function Determinant
9: ( mat : Matrix; rows,cols : Standard_Natural_Vectors.Vector )
10: return double_float is
11:
12: -- DESCRIPTION :
13: -- Computes the determinant of the matrix selecting d rows and columns,
14: -- where d = rows'length = cols'length.
15:
16: res : double_float := 1.0;
17: sqm : Matrix(rows'range,rows'range);
18: piv : Standard_Natural_Vectors.Vector(rows'range);
19: inf : natural;
20:
21: begin
22: for i in rows'range loop
23: piv(i) := i;
24: for j in cols'range loop
25: sqm(i,j) := mat(rows(i),cols(j));
26: end loop;
27: end loop;
28: lufac(sqm,rows'last,piv,inf);
29: for i in rows'range loop
30: res := res*sqm(i,i);
31: end loop;
32: for i in piv'range loop
33: if piv(i) > i
34: then res := -res;
35: end if;
36: end loop;
37: return res;
38: end Determinant;
39:
40: procedure Minors ( file : in file_type;
41: n,m : in natural; mat : in Matrix ) is
42:
43: procedure Minors ( d : in natural ) is
44:
45: rows,cols : Standard_Natural_Vectors.Vector(1..d);
46:
47: procedure Select_Rows ( k,start : in natural ) is
48:
49: det : double_float;
50:
51: begin
52: if k > d
53: then det := Determinant(mat,rows,cols);
54: put(file,"Minor ");
55: put(file,rows); put(file," X"); put(file,cols);
56: put(file," equals ");
57: put(file,det); new_line(file);
58: else for j in start..n loop
59: rows(k) := j;
60: Select_Rows(k+1,j+1);
61: end loop;
62: end if;
63: end Select_Rows;
64:
65: procedure Select_Columns ( k,start : in natural ) is
66: begin
67: if k > d
68: then Select_Rows(1,1);
69: else for j in start..m loop
70: cols(k) := j;
71: Select_Columns(k+1,j+1);
72: end loop;
73: end if;
74: end Select_Columns;
75:
76: begin
77: Select_Columns(1,1);
78: end Minors;
79:
80: begin
81: for d in 2..m loop
82: Minors(d);
83: end loop;
84: end Minors;
85:
86: function Number_of_Minors ( n,m : natural ) return natural is
87:
88: res : natural := 0;
89:
90: procedure Count_Minors ( d : in natural ) is
91:
92: procedure Select_Rows ( k,start : in natural ) is
93: begin
94: if k > d
95: then res := res + 1;
96: else for j in start..n loop
97: Select_Rows(k+1,j+1);
98: end loop;
99: end if;
100: end Select_Rows;
101:
102: procedure Select_Columns ( k,start : in natural ) is
103: begin
104: if k > d
105: then Select_Rows(1,1);
106: else for j in start..m loop
107: Select_Columns(k+1,j+1);
108: end loop;
109: end if;
110: end Select_Columns;
111:
112: begin
113: Select_Columns(1,1);
114: end Count_Minors;
115:
116: begin
117: for d in 2..m loop
118: Count_Minors(d);
119: end loop;
120: return res;
121: end Number_of_Minors;
122:
123: function Sign_of_Minors ( n,m : natural; mat : Matrix ) return Vector is
124:
125: res : Vector(1..Number_of_Minors(n,m));
126: tol : constant double_float := 10.0**(-10);
127: cnt : natural := 0;
128:
129: procedure Compute_Minors ( d : in natural ) is
130:
131: rows,cols : Standard_Natural_Vectors.Vector(1..d);
132:
133: procedure Select_Rows ( k,start : in natural ) is
134:
135: det : double_float;
136:
137: begin
138: if k > d
139: then det := Determinant(mat,rows,cols);
140: cnt := cnt + 1;
141: if abs(det) < tol
142: then res(cnt) := 0;
143: elsif det > 0.0
144: then res(cnt) := +1;
145: else res(cnt) := -1;
146: end if;
147: else for j in start..n loop
148: rows(k) := j;
149: Select_Rows(k+1,j+1);
150: end loop;
151: end if;
152: end Select_Rows;
153:
154: procedure Select_Columns ( k,start : in natural ) is
155: begin
156: if k > d
157: then Select_Rows(1,1);
158: else for j in start..m loop
159: cols(k) := j;
160: Select_Columns(k+1,j+1);
161: end loop;
162: end if;
163: end Select_Columns;
164:
165: begin
166: Select_Columns(1,1);
167: end Compute_Minors;
168:
169: begin
170: for d in 2..m loop
171: Compute_Minors(d);
172: end loop;
173: return res;
174: end Sign_of_Minors;
175:
176: end Minor_Computations;