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;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>