[BACK]Return to ts_fltdls.adb CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices

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>