Return to ts_givrot.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports

Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Supports/ts_givrot.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_Integer_Vectors;
                      4: with Standard_Integer_Vectors_io;        use Standard_Integer_Vectors_io;
                      5: with Standard_Floating_Vectors;          use Standard_Floating_Vectors;
                      6: with Standard_Floating_Vectors_io;       use Standard_Floating_Vectors_io;
                      7: with Standard_Floating_Matrices;         use Standard_Floating_Matrices;
                      8: with Standard_Floating_Matrices_io;      use Standard_Floating_Matrices_io;
                      9: with Givens_Rotations;                   use Givens_Rotations;
                     10:
                     11: procedure ts_givrot is
                     12:
                     13: -- DESCRIPTION :
                     14: --   Interactive testing of Givens rotations.
                     15:
                     16:   n,m : natural;
                     17:   ans : character;
                     18:   tol : constant double_float := 10.0**(-12);
                     19:
                     20:   function Compute_Residual ( mat : in matrix; rhs,x : in vector;
                     21:                               ipvt : in Standard_Integer_Vectors.Vector )
                     22:                             return Vector is
                     23:
                     24:     res : Vector(rhs'range) := rhs;
                     25:
                     26:   begin
                     27:     for i in mat'range(1) loop
                     28:       for j in mat'range(2) loop
                     29:         res(i) := res(i) - mat(i,ipvt(j))*x(j);
                     30:         exit when j > mat'last(2);
                     31:       end loop;
                     32:     end loop;
                     33:     return res;
                     34:   end Compute_Residual;
                     35:
                     36:   procedure Row_Wise_Triangulation ( mat : in out Matrix ) is
                     37:
                     38:     pivots : Standard_Integer_Vectors.Vector(mat'range(2));
                     39:     ipiv : integer;
                     40:
                     41:   begin
                     42:     for i in pivots'range loop
                     43:       pivots(i) := i;
                     44:     end loop;
                     45:     for i in mat'range(1) loop
                     46:       Upper_Triangulate(i,mat,tol,pivots,ipiv);
                     47:       put("At row "); put(i,1); put(" we have pivot : "); put(ipiv,1);
                     48:       new_line;
                     49:       put_line("with matrix : "); put(mat);
                     50:       put("and pivot vector "); put(pivots); new_line;
                     51:       exit when (ipiv = 0);
                     52:     end loop;
                     53:   end Row_Wise_Triangulation;
                     54:
                     55: begin
                     56:   loop
                     57:     put("Give the number of rows : "); get(n);
                     58:     put("Give the number of columns : "); get(m);
                     59:     declare
                     60:       mat,wrkmat : matrix(1..n,1..m);
                     61:       rhs,wrkrhs : vector(1..n);
                     62:       ipvt : Standard_Integer_Vectors.Vector(1..m);
                     63:     begin
                     64:       put("Give the elements of the ");
                     65:       put(n,1); put("x"); put(m,1); put_line("-matrix :"); get(mat);
                     66:       put_line("The matrix :"); put(mat);
                     67:       wrkmat := mat;
                     68:       put("Do you have a right hand side ? (y/n) "); get(ans);
                     69:       if ans = 'y'
                     70:        then put("Give "); put(n,1); put(" floating point numbers : ");
                     71:             rhs := (1..n => 0.0); get(rhs);
                     72:             wrkrhs := rhs;
                     73:             Upper_Triangulate(wrkmat,wrkrhs,tol,ipvt);
                     74:        else Row_Wise_Triangulation(wrkmat);
                     75:             wrkmat := mat;
                     76:             Upper_Triangulate(wrkmat,tol,ipvt);
                     77:       end if;
                     78:       put_line("The matrix in upper triangular form : "); put(wrkmat);
                     79:       put(" with pivoting information : "); put(ipvt); new_line;
                     80:       if ans = 'y'
                     81:        then put(" and transformed right hand side : ");
                     82:             put(wrkrhs); new_line;
                     83:             declare
                     84:               sol,res : Vector(1..n);
                     85:             begin
                     86:               Solve(wrkmat,wrkrhs,tol,sol);
                     87:               put("The solution : "); put(sol); new_line;
                     88:               res := Compute_Residual(mat,rhs,sol,ipvt);
                     89:               put(" with residual : "); put(res); new_line;
                     90:             end;
                     91:       end if;
                     92:     end;
                     93:     put("Do you want more tests ? (y/n) "); get(ans);
                     94:     exit when ans /= 'y';
                     95:   end loop;
                     96: end ts_givrot;