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     ! 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;