Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Supports/givens_rotations.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
2:
3: package body Givens_Rotations is
4:
5: procedure Givens_Factors ( v: in Vector; i,j : in integer;
6: cosi,sine : out double_float ) is
7: norm : double_float;
8:
9: begin
10: norm := SQRT(v(i)*v(i) + v(j)*v(j));
11: cosi := v(i)/norm; sine := v(j)/norm;
12: end Givens_Factors;
13:
14: procedure Givens_Factors ( mat : in Matrix; i,j : in integer;
15: cosi,sine : out double_float ) is
16: norm : double_float;
17:
18: begin
19: norm := SQRT(mat(i,i)*mat(i,i) + mat(j,i)*mat(j,i));
20: cosi := mat(i,i)/norm; sine := mat(j,i)/norm;
21: end Givens_Factors;
22:
23: procedure Givens_Rotation ( v : in out Vector; i,j : in integer;
24: cosi,sine : in double_float ) is
25: temp : double_float;
26:
27: begin
28: temp := v(i);
29: v(i) := cosi*v(i) + sine*v(j);
30: v(j) := -sine*temp + cosi*v(j);
31: end Givens_Rotation;
32:
33: procedure Givens_Rotation ( mat : in out Matrix; i,j : in integer;
34: cosi,sine : in double_float ) is
35: temp : double_float;
36:
37: begin
38: -- for k in mat'range(2) loop -- certain angle preservation
39: for k in i..mat'last(2) loop -- only angle preservation if upper triangular
40: temp := mat(i,k);
41: mat(i,k) := cosi*mat(i,k) + sine*mat(j,k);
42: mat(j,k) := -sine*temp + cosi*mat(j,k);
43: end loop;
44: end Givens_Rotation;
45:
46: procedure Givens_Rotation ( mat : in out Matrix; lastcol,i,j : in integer;
47: cosi,sine : in double_float ) is
48: temp : double_float;
49:
50: begin
51: -- for k in mat'first(2)..lastcol loop -- certain angle preservation
52: for k in i..lastcol loop -- only angle preservation if upper triangular
53: temp := mat(i,k);
54: mat(i,k) := cosi*mat(i,k) + sine*mat(j,k);
55: mat(j,k) := -sine*temp + cosi*mat(j,k);
56: end loop;
57: end Givens_Rotation;
58:
59: procedure Givens_Rotation ( v : in out Vector; i,j : in integer ) is
60:
61: cosi,sine : double_float;
62:
63: begin
64: Givens_Factors(v,i,j,cosi,sine);
65: Givens_Rotation(v,i,j,cosi,sine);
66: end Givens_Rotation;
67:
68: procedure Givens_Rotation ( mat : in out Matrix; i,j : in integer ) is
69:
70: cosi,sine : double_float;
71:
72: begin
73: Givens_Factors(mat,i,j,cosi,sine);
74: Givens_Rotation(mat,i,j,cosi,sine);
75: end Givens_Rotation;
76:
77: procedure Givens_Rotation ( mat : in out Matrix;
78: lastcol,i,j : in integer ) is
79:
80: cosi,sine : double_float;
81:
82: begin
83: Givens_Factors(mat,i,j,cosi,sine);
84: Givens_Rotation(mat,lastcol,i,j,cosi,sine);
85: end Givens_Rotation;
86:
87: procedure Upper_Triangulate
88: ( row : in integer; mat : in out Matrix; tol : in double_float;
89: ipvt : in out Standard_Integer_Vectors.Vector;
90: pivot : out integer ) is
91:
92: piv,tpi : integer := 0;
93: tmp : double_float;
94:
95: begin
96: for j in row..mat'last(2) loop -- search pivot
97: if abs(mat(row,j)) > tol
98: then piv := j;
99: end if;
100: exit when (piv /= 0);
101: end loop;
102: if piv /= 0 -- zero row
103: then
104: if piv /= row -- interchange columns
105: then for k in mat'range(1) loop
106: tmp := mat(k,row); mat(k,row) := mat(k,piv); mat(k,piv) := tmp;
107: end loop;
108: tpi := ipvt(row); ipvt(row) := ipvt(piv); ipvt(piv) := tpi;
109: end if;
110: for k in row+1..mat'last(1) loop -- perform Givens rotations
111: if abs(mat(k,row)) > tol
112: then Givens_Rotation(mat,row,k);
113: end if;
114: end loop;
115: end if;
116: pivot := piv;
117: end Upper_Triangulate;
118:
119: procedure Upper_Triangulate ( mat : in out Matrix; tol : in double_float ) is
120:
121: pivot : integer;
122: temp : double_float;
123:
124: begin
125: for i in mat'range(1) loop -- mat(mat'first(1)..i-1,mat'first(2)..i-1)
126: pivot := 0; -- is already upper triangular
127: for j in i..mat'last(2) loop -- search pivot
128: if abs(mat(i,j)) > tol
129: then pivot := j;
130: end if;
131: exit when (pivot /= 0);
132: end loop;
133: exit when (pivot = 0); -- zero row
134: if pivot /= i -- interchange columns
135: then for k in mat'range(1) loop
136: temp := mat(k,i); mat(k,i) := mat(k,pivot); mat(k,pivot) := temp;
137: end loop;
138: end if;
139: for k in i+1..mat'last(1) loop -- perform Givens rotations
140: if abs(mat(k,i)) > tol
141: then Givens_Rotation(mat,i,k);
142: end if;
143: end loop; -- mat(mat'first(1)..i,mat'first(2)..i)
144: exit when i > mat'last(2); -- is upper triangular
145: end loop;
146: end Upper_Triangulate;
147:
148: procedure Upper_Triangulate ( mat : in out Matrix; tol : in double_float;
149: ipvt : out Standard_Integer_Vectors.Vector ) is
150:
151: pivot,tpi : integer;
152: pivots : Standard_Integer_Vectors.Vector(mat'range(2));
153: temp : double_float;
154:
155: begin
156: for i in pivots'range loop -- initialize vector of pivots
157: pivots(i) := i;
158: end loop;
159: for i in mat'range(1) loop -- mat(mat'first(1)..i-1,mat'first(2)..i-1)
160: pivot := 0; -- is already upper triangular
161: for j in i..mat'last(2) loop -- search pivot
162: if abs(mat(i,j)) > tol
163: then pivot := j;
164: end if;
165: exit when (pivot /= 0);
166: end loop;
167: exit when (pivot = 0); -- zero row
168: if pivot /= i -- interchange columns
169: then for k in mat'range(1) loop
170: temp := mat(k,i); mat(k,i) := mat(k,pivot); mat(k,pivot) := temp;
171: end loop;
172: tpi := pivots(i); pivots(i) := pivots(pivot); pivots(pivot) := tpi;
173: end if;
174: for k in i+1..mat'last(1) loop -- perform Givens rotations
175: if abs(mat(k,i)) > tol
176: then Givens_Rotation(mat,i,k);
177: end if;
178: end loop; -- mat(mat'first(1)..i,mat'first(2)..i)
179: exit when i > mat'last(2); -- is upper triangular
180: end loop;
181: ipvt := pivots;
182: end Upper_Triangulate;
183:
184: procedure Upper_Triangulate ( mat : in out Matrix; rhs : in out Vector;
185: tol : in double_float ) is
186:
187: pivot : integer;
188: temp,cosi,sine : double_float;
189:
190: begin
191: for i in mat'range(1) loop -- mat(mat'first(1)..i-1,mat'first(2)..i-1)
192: pivot := 0; -- is already upper triangular
193: for j in i..mat'last(2) loop -- search pivot
194: if abs(mat(i,j)) > tol
195: then pivot := j;
196: end if;
197: exit when (pivot /= 0);
198: end loop;
199: exit when (pivot = 0); -- zero row
200: if pivot /= i -- interchange columns
201: then for k in mat'range(1) loop
202: temp := mat(k,i); mat(k,i) := mat(k,pivot); mat(k,pivot) := temp;
203: end loop;
204: end if;
205: for k in i+1..mat'last(1) loop -- perform Givens rotations
206: if abs(mat(k,i)) > tol
207: then Givens_Factors(mat,i,k,cosi,sine);
208: Givens_Rotation(mat,i,k,cosi,sine);
209: Givens_Rotation(rhs,i,k,cosi,sine);
210: end if;
211: end loop; -- mat(mat'first(1)..i,mat'first(2)..i)
212: exit when i > mat'last(2); -- is upper triangular
213: end loop;
214: end Upper_Triangulate;
215:
216: procedure Upper_Triangulate ( mat : in out Matrix; rhs : in out Vector;
217: tol : in double_float;
218: ipvt : out Standard_Integer_Vectors.Vector ) is
219:
220: pivot,tpi : integer;
221: pivots : Standard_Integer_Vectors.Vector(mat'range(2));
222: temp,cosi,sine : double_float;
223:
224: begin
225: for i in pivots'range loop -- initialize vector of pivots
226: pivots(i) := i;
227: end loop;
228: for i in mat'range(1) loop -- mat(mat'first(1)..i-1,mat'first(2)..i-1)
229: pivot := 0; -- is already upper triangular
230: for j in i..mat'last(2) loop -- search pivot
231: if abs(mat(i,j)) > tol
232: then pivot := j;
233: end if;
234: exit when (pivot /= 0);
235: end loop;
236: exit when (pivot = 0); -- zero row
237: if pivot /= i -- interchange columns
238: then for k in mat'range(1) loop
239: temp := mat(k,i); mat(k,i) := mat(k,pivot); mat(k,pivot) := temp;
240: end loop;
241: tpi := pivots(i); pivots(i) := pivots(pivot); pivots(pivot) := tpi;
242: end if;
243: for k in i+1..mat'last(1) loop -- perform Givens rotations
244: if abs(mat(k,i)) > tol
245: then Givens_Factors(mat,i,k,cosi,sine);
246: Givens_Rotation(mat,i,k,cosi,sine);
247: Givens_Rotation(rhs,i,k,cosi,sine);
248: end if;
249: end loop; -- mat(mat'first(1)..i,mat'first(2)..i)
250: exit when i > mat'last(2); -- is upper triangular
251: end loop;
252: ipvt := pivots;
253: end Upper_Triangulate;
254:
255: procedure Solve ( mat : in Matrix; rhs : in Vector; tol : in double_float;
256: x : out Vector ) is
257:
258: rank : natural := 0;
259: sol : vector(mat'range(1)) := (mat'range(1) => 0.0);
260:
261: begin
262: for i in mat'range(1) loop -- determine the rank of the matrix
263: if abs(mat(i,i)) > tol
264: then rank := rank + 1;
265: end if;
266: exit when ((i >= mat'last(2)) or (abs(mat(i,i)) < tol));
267: end loop;
268: for i in reverse mat'first(1)..rank loop -- back substitution
269: for j in i+1..rank loop
270: sol(i) := sol(i) + mat(i,j)*sol(j);
271: end loop;
272: sol(i) := rhs(i) - sol(i);
273: if abs(mat(i,i)) > tol
274: then sol(i) := sol(i)/mat(i,i);
275: elsif abs(sol(i)) < tol
276: then sol(i) := 1.0;
277: else return;
278: end if;
279: end loop; -- invariant : sol(i..rank) computed
280: x := sol;
281: end Solve;
282:
283: procedure Solve ( mat : in Matrix; rhs : in Vector; tol : in double_float;
284: rank : in natural; x : out Vector ) is
285:
286: sol : vector(mat'range(1)) := (mat'range(1) => 0.0);
287:
288: begin
289: for i in reverse mat'first(1)..rank loop -- back substitution
290: for j in i+1..rank loop
291: sol(i) := sol(i) + mat(i,j)*sol(j);
292: end loop;
293: sol(i) := rhs(i) - sol(i);
294: if abs(mat(i,i)) > tol
295: then sol(i) := sol(i)/mat(i,i);
296: elsif abs(sol(i)) < tol
297: then sol(i) := 1.0;
298: else return;
299: end if;
300: end loop; -- invariant : sol(i..rank) computed
301: x := sol;
302: end Solve;
303:
304: end Givens_Rotations;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to givens_rotations.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports