Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/ts_qrd.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_Complex_Numbers; use Standard_Complex_Numbers;
5: with Standard_Natural_Vectors;
6: with Standard_Natural_Vectors_io; use Standard_Natural_Vectors_io;
7: with Standard_Floating_Vectors;
8: with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
9: with Standard_Floating_Matrices;
10: with Standard_Floating_Matrices_io; use Standard_Floating_Matrices_io;
11: with Standard_Random_Vectors; use Standard_Random_Vectors;
12: with Standard_Random_Matrices; use Standard_Random_Matrices;
13: with Standard_Floating_QR_Decomposition; use Standard_Floating_QR_Decomposition;
14: with Standard_Floating_Least_Squares; use Standard_Floating_Least_Squares;
15: with Standard_Complex_Vectors;
16: with Standard_Complex_Vectors_io; use Standard_Complex_Vectors_io;
17: with Standard_Complex_Matrices;
18: with Standard_Complex_Matrices_io; use Standard_Complex_Matrices_io;
19: with Standard_Complex_QR_Decomposition; use Standard_Complex_QR_Decomposition;
20: with Standard_Complex_Least_Squares; use Standard_Complex_Least_Squares;
21:
22: procedure ts_qrd is
23:
24: -- DESCRIPTION :
25: -- This program tests the implementation of the QR decomposition
26: -- and least squares approximation.
27:
28: -- GENERAL TESTS ON QR DECOMPOSITION :
29:
30: function Extract_Upper_Triangular
31: ( a : Standard_Floating_Matrices.Matrix )
32: return Standard_Floating_Matrices.Matrix is
33:
34: -- DESCRIPTION :
35: -- Returns the upper triangular part of the matrix a.
36:
37: res : Standard_Floating_Matrices.Matrix(a'range(1),a'range(2));
38:
39: begin
40: for i in a'range(1) loop
41: for j in a'first(2)..(i-1) loop
42: res(i,j) := 0.0;
43: end loop;
44: for j in i..a'last(2) loop
45: res(i,j) := a(i,j);
46: end loop;
47: end loop;
48: return res;
49: end Extract_Upper_Triangular;
50:
51: function Extract_Upper_Triangular
52: ( a : Standard_Complex_Matrices.Matrix )
53: return Standard_Complex_Matrices.Matrix is
54:
55: -- DESCRIPTION :
56: -- Returns the upper triangular part of the matrix a.
57:
58: res : Standard_Complex_Matrices.Matrix(a'range(1),a'range(2));
59:
60: begin
61: for i in a'range(1) loop
62: for j in a'first(2)..(i-1) loop
63: res(i,j) := Create(0.0);
64: end loop;
65: for j in i..a'last(2) loop
66: res(i,j) := a(i,j);
67: end loop;
68: end loop;
69: return res;
70: end Extract_Upper_Triangular;
71:
72: function Differences ( a,b : in Standard_Floating_Matrices.Matrix )
73: return double_float is
74:
75: -- DESCRIPTION :
76: -- Returns the sum of the differences of all elements |a(i,j)-b(i,j)|.
77:
78: sum : double_float := 0.0;
79:
80: begin
81: for i in a'range(1) loop
82: for j in a'range(2) loop
83: sum := sum + abs(a(i,j)-b(i,j));
84: end loop;
85: end loop;
86: return sum;
87: end Differences;
88:
89: function Differences ( a,b : in Standard_Complex_Matrices.Matrix )
90: return double_float is
91:
92: -- DESCRIPTION :
93: -- Returns the sum of the differences of all elements |a(i,j)-b(i,j)|.
94:
95: sum : double_float := 0.0;
96:
97: begin
98: for i in a'range(1) loop
99: for j in a'range(2) loop
100: sum := sum + AbsVal(a(i,j)-b(i,j));
101: end loop;
102: end loop;
103: return sum;
104: end Differences;
105:
106: procedure Orthogonality ( q : in Standard_Floating_Matrices.Matrix ) is
107:
108: -- DESCRIPTION :
109: -- Tests whether the columns are orthogonal w.r.t. each other.
110:
111: sum,ip : double_float;
112:
113: begin
114: sum := 0.0;
115: for j in q'range(2) loop
116: for k in j+1..q'last(2) loop
117: ip := 0.0;
118: for i in q'range(1) loop
119: ip := ip + q(i,j)*q(i,k);
120: end loop;
121: sum := sum + abs(ip);
122: end loop;
123: end loop;
124: put("Orthogonality check : "); put(sum,3,3,3); new_line;
125: end Orthogonality;
126:
127: procedure Orthogonality ( q : in Standard_Complex_Matrices.Matrix ) is
128:
129: -- DESCRIPTION :
130: -- Tests whether the columns are orthogonal w.r.t. each other.
131:
132: sum : double_float := 0.0;
133: ip : Complex_Number;
134:
135: begin
136: for j in q'range(2) loop
137: for k in j+1..q'last(2) loop
138: ip := Create(0.0);
139: for i in q'range(1) loop
140: ip := ip + Conjugate(q(i,j))*q(i,k);
141: end loop;
142: sum := sum + AbsVal(ip);
143: end loop;
144: end loop;
145: put("Orthogonality check : "); put(sum,3,3,3); new_line;
146: end Orthogonality;
147:
148: procedure Test_QRD ( a,q,r : in Standard_Floating_Matrices.Matrix ) is
149:
150: wrk : Standard_Floating_Matrices.Matrix(a'range(1),a'range(2));
151: use Standard_Floating_Matrices;
152:
153: begin
154: put_line("The upper triangular part R :"); put(r,3);
155: wrk := q*r;
156: put_line("q*r :"); put(wrk,3);
157: put("Difference in 1-norm between the matrix and q*r: ");
158: put(Differences(a,wrk),3,3,3); new_line;
159: Orthogonality(q);
160: end Test_QRD;
161:
162: procedure Test_QRD ( a,q,r : in Standard_Complex_Matrices.Matrix ) is
163:
164: wrk : Standard_Complex_Matrices.Matrix(a'range(1),a'range(2));
165: use Standard_Complex_Matrices;
166:
167: begin
168: put_line("The upper triangular part R :"); put(r,3);
169: wrk := q*r;
170: put_line("q*r :"); put(wrk,3);
171: put("Difference in 1-norm between the matrix and q*r : ");
172: put(Differences(a,wrk),3,3,3); new_line;
173: Orthogonality(q);
174: end Test_QRD;
175:
176: -- REAL TEST DRIVERS :
177:
178: procedure Real_LS_Test ( n,m : in natural; piv : in boolean;
179: a : in Standard_Floating_Matrices.Matrix;
180: b : in Standard_Floating_Vectors.Vector ) is
181:
182: wrk : Standard_Floating_Matrices.Matrix(1..n,1..m) := a;
183: qraux : Standard_Floating_Vectors.Vector(1..m) := (1..m => 0.0);
184: jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);
185: sol : Standard_Floating_Vectors.Vector(1..m);
186: rsd,dum : Standard_Floating_Vectors.Vector(1..n);
187: info : integer;
188: use Standard_Floating_Matrices;
189: use Standard_Floating_Vectors;
190:
191: begin
192: put_line("The matrix : "); put(a,3);
193: QRD(wrk,qraux,jpvt,piv);
194: put_line("The matrix after QR : "); put(wrk,3);
195: put_line("The vector qraux : "); put(qraux,3); new_line;
196: if piv
197: then put("The vector jpvt : "); put(jpvt); new_line;
198: Permute_Columns(wrk,jpvt);
199: end if;
200: QRLS(wrk,n,n,m,qraux,b,dum,dum,sol,rsd,dum,110,info);
201: if piv
202: then Permute(sol,jpvt);
203: end if;
204: put_line("The solution : "); put(sol,3); new_line;
205: put_line("The residual : "); put(rsd,3); new_line;
206: dum := b - a*sol;
207: put_line("right-hand size - matrix*solution : "); put(dum,3); new_line;
208: end Real_LS_Test;
209:
210: procedure Real_QR_Test ( n,m : in natural; piv : in boolean;
211: a : in Standard_Floating_Matrices.Matrix ) is
212:
213: wrk : Standard_Floating_Matrices.Matrix(1..n,1..m) := a;
214: bas : Standard_Floating_Matrices.Matrix(1..n,1..n);
215: qraux : Standard_Floating_Vectors.Vector(1..m) := (1..m => 0.0);
216: jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);
217:
218: begin
219: put_line("The matrix : "); put(a,3);
220: QRD(wrk,qraux,jpvt,piv);
221: put_line("The matrix after QR : "); put(wrk,3);
222: put_line("The vector qraux : "); put(qraux,3); new_line;
223: if piv
224: then put("The vector jpvt : "); put(jpvt); new_line;
225: Permute_Columns(wrk,jpvt);
226: end if;
227: for i in wrk'range(1) loop
228: for j in wrk'range(2) loop
229: bas(i,j) := wrk(i,j);
230: end loop;
231: for j in n+1..m loop
232: bas(i,j) := 0.0;
233: end loop;
234: end loop;
235: Basis(bas,a);
236: put_line("The orthogonal part Q of QR :"); put(bas,3);
237: Test_QRD(a,bas,Extract_Upper_Triangular(wrk));
238: end Real_QR_Test;
239:
240: procedure Interactive_Real_QR_Test ( n,m : in natural; piv : in boolean ) is
241:
242: a : Standard_Floating_Matrices.Matrix(1..n,1..m);
243:
244: begin
245: put("Give a "); put(n,1); put("x"); put(m,1);
246: put_line(" matrix : "); get(a);
247: Real_QR_Test(n,m,piv,a);
248: end Interactive_Real_QR_Test;
249:
250: procedure Interactive_Real_LS_Test ( n,m : in natural; piv : in boolean ) is
251:
252: a : Standard_Floating_Matrices.Matrix(1..n,1..m);
253: b : Standard_Floating_Vectors.Vector(1..n);
254:
255: begin
256: put("Give a "); put(n,1); put("x"); put(m,1);
257: put_line(" matrix : "); get(a);
258: put("Give right-hand size "); put(n,1);
259: put_line("-vector : "); get(b);
260: Real_LS_Test(n,m,piv,a,b);
261: end Interactive_Real_LS_Test;
262:
263: procedure Random_Real_QR_Test ( n,m : in natural; piv : in boolean ) is
264:
265: a : Standard_Floating_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);
266:
267: begin
268: Real_QR_Test(n,m,piv,a);
269: end Random_Real_QR_Test;
270:
271: procedure Random_Real_LS_Test ( n,m : in natural; piv : in boolean ) is
272:
273: a : Standard_Floating_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);
274: b : Standard_Floating_Vectors.Vector(1..n) := Random_Vector(n);
275:
276: begin
277: Real_LS_Test(n,m,piv,a,b);
278: end Random_Real_LS_Test;
279:
280: -- COMPLEX TEST DRIVERS :
281:
282: procedure Complex_QR_Test ( n,m : in natural; piv : in boolean;
283: a : Standard_Complex_Matrices.Matrix ) is
284:
285: wrk : Standard_Complex_Matrices.Matrix(1..n,1..m) := a;
286: bas : Standard_Complex_Matrices.Matrix(1..n,1..n);
287: qraux : Standard_Complex_Vectors.Vector(1..m) := (1..m => Create(0.0));
288: jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);
289:
290: begin
291: put_line("The matrix : "); put(a,3);
292: QRD(wrk,qraux,jpvt,piv);
293: put_line("The matrix after QR : "); put(wrk,3);
294: put_line("The vector qraux : "); put(qraux,3); new_line;
295: -- put("The vector jpvt : "); put(jpvt); new_line;
296: if not piv
297: then for i in wrk'range(1) loop
298: for j in wrk'range(2) loop
299: bas(i,j) := wrk(i,j);
300: end loop;
301: for j in n+1..m loop
302: bas(i,j) := Create(0.0);
303: end loop;
304: end loop;
305: Basis(bas,a);
306: put_line("The orthogonal part Q of QR :"); put(bas,3);
307: Test_QRD(a,bas,Extract_Upper_Triangular(wrk));
308: end if;
309: end Complex_QR_Test;
310:
311: procedure Complex_LS_Test ( n,m : in natural; piv : in boolean;
312: a : Standard_Complex_Matrices.Matrix;
313: b : Standard_Complex_Vectors.Vector ) is
314:
315: wrk : Standard_Complex_Matrices.Matrix(1..n,1..m) := a;
316: qraux : Standard_Complex_Vectors.Vector(1..m) := (1..m => Create(0.0));
317: jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);
318: sol : Standard_Complex_Vectors.Vector(1..m);
319: rsd,dum : Standard_Complex_Vectors.Vector(1..n);
320: info : integer;
321: use Standard_Complex_Matrices;
322: use Standard_Complex_Vectors;
323:
324: begin
325: put_line("The matrix : "); put(a,3);
326: QRD(wrk,qraux,jpvt,piv);
327: put_line("The matrix after QR : "); put(wrk,3);
328: put_line("The vector qraux : "); put(qraux,3); new_line;
329: -- put("The vector jpvt : "); put(jpvt); new_line;
330: QRLS(wrk,n,n,m,qraux,b,dum,dum,sol,rsd,dum,110,info);
331: put_line("The solution : "); put(sol,3); new_line;
332: put_line("The residual : "); put(rsd,3); new_line;
333: dum := b - a*sol;
334: put_line("right-hand size - matrix*solution : "); put(dum,3); new_line;
335: end Complex_LS_Test;
336:
337: procedure Interactive_Complex_QR_Test
338: ( n,m : in natural; piv : in boolean ) is
339:
340: a : Standard_Complex_Matrices.Matrix(1..n,1..m);
341:
342: begin
343: put("Give a "); put(n,1); put("x"); put(m,1);
344: put_line(" matrix : "); get(a);
345: Complex_QR_Test(n,m,piv,a);
346: end Interactive_Complex_QR_Test;
347:
348: procedure Interactive_Complex_LS_Test
349: ( n,m : in natural; piv : in boolean ) is
350:
351: a : Standard_Complex_Matrices.Matrix(1..n,1..m);
352: b : Standard_Complex_Vectors.Vector(1..n);
353:
354: begin
355: put("Give a "); put(n,1); put("x"); put(m,1);
356: put_line(" matrix : "); get(a);
357: put("Give right-hand size "); put(n,1);
358: put_line("-vector : "); get(b);
359: Complex_LS_Test(n,m,piv,a,b);
360: end Interactive_Complex_LS_Test;
361:
362: procedure Random_Complex_QR_Test ( n,m : in natural; piv : in boolean ) is
363:
364: a : Standard_Complex_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);
365:
366: begin
367: Complex_QR_Test(n,m,piv,a);
368: end Random_Complex_QR_Test;
369:
370: procedure Random_Complex_LS_Test ( n,m : in natural; piv : in boolean ) is
371:
372: a : Standard_Complex_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);
373: b : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(n);
374:
375: begin
376: Complex_LS_Test(n,m,piv,a,b);
377: end Random_Complex_LS_Test;
378:
379: -- MAIN PROGRAM :
380:
381: procedure Main is
382:
383: n,m : natural;
384: ans,choice : character;
385: piv : boolean := false;
386:
387: begin
388: new_line;
389: put_line("Test on the QR decomposition and Least Squares.");
390: loop
391: new_line;
392: put_line("Choose one of the following : ");
393: put_line(" 0. Exit this program.");
394: put_line(" 1. QR-decomposition on user-given floating-point matrix.");
395: put_line(" 2. complex matrix.");
396: put_line(" 3. on random floating-point matrix.");
397: put_line(" 4. complex matrix.");
398: put_line(" 5. Least Squares on user-given floating-point matrix.");
399: put_line(" 6. complex matrix.");
400: put_line(" 7. on random floating-point matrix.");
401: put_line(" 8. complex matrix.");
402: put("Make your choice (0,1,2,3,4,5,6,7 or 8) : "); get(choice);
403: exit when (choice = '0');
404: put("Give the number of rows : "); get(n);
405: put("Give the number of columns : "); get(m);
406: put("Do you want pivoting ? (y/n) "); get(ans);
407: piv := (ans = 'y');
408: case choice is
409: when '1' => Interactive_Real_QR_Test(n,m,piv);
410: when '2' => Interactive_Complex_QR_Test(n,m,piv);
411: when '3' => Random_Real_QR_Test(n,m,piv);
412: when '4' => Random_Complex_QR_Test(n,m,piv);
413: when '5' => Interactive_Real_LS_Test(n,m,piv);
414: when '6' => Interactive_Complex_LS_Test(n,m,piv);
415: when '7' => Random_Real_LS_Test(n,m,piv);
416: when '8' => Random_Complex_LS_Test(n,m,piv);
417: when others => null;
418: end case;
419: end loop;
420: end Main;
421:
422: begin
423: Main;
424: end ts_qrd;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ts_qrd.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 / Matrices