Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/ts_qrd.adb, Revision 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