Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/generic_floating_linear_solvers.ads, Revision 1.1
1.1 ! maekawa 1: with Standard_Natural_Vectors; use Standard_Natural_Vectors;
! 2: with Abstract_Ring,Abstract_Ring.Field;
! 3: with Generic_Vectors,Generic_Matrices;
! 4:
! 5: generic
! 6:
! 7: with package Ring is new Abstract_Ring(<>);
! 8: with package Field is new Ring.Field(<>);
! 9: with package Vectors is new Generic_Vectors(Ring);
! 10: with package Matrices is new Generic_Matrices(Ring,Vectors);
! 11:
! 12: package Generic_Floating_Linear_Solvers is
! 13:
! 14: -- DESCRIPTION :
! 15: -- This package offers a few routines to solve linear systems of equations.
! 16: -- The code for lufac, lufco and lusolve is a literal translation from the
! 17: -- f77-linpack code. We distinguish between static triangulators, which
! 18: -- take the matrix as a whole at once, and dynamic triangulators, which
! 19: -- allow to triangulate row after row.
! 20:
! 21: use Ring,Matrices;
! 22:
! 23: -- STATIC TRIANGULATORS :
! 24:
! 25: procedure lufac ( a : in out Matrix; n : in natural;
! 26: ipvt : out Standard_Natural_Vectors.Vector;
! 27: info : out natural );
! 28:
! 29: -- DESCRIPTION :
! 30: -- lufac factors a float matrix by gaussian elimination
! 31:
! 32: -- lufac is usually called by lufco, but it can be called
! 33: -- directly with a saving of time if rcond is not needed.
! 34: -- (time for lufco) = (1 + 9/n)*(time for lufac).
! 35:
! 36: -- ON ENTRY :
! 37: -- a a floating point matrix(1..n,1..n) to be factored
! 38: -- n the dimension of the matrix a
! 39:
! 40: -- ON RETURN :
! 41: -- a an upper triangular matrix and the multipliers
! 42: -- which were used to obtain it.
! 43: -- The factorization can be written a = l*u where
! 44: -- l is a product of permutation and unit lower
! 45: -- triangular matrices and u is upper triangular.
! 46: -- ipvt a natural vector of pivot indices
! 47: -- info = 0 normal value
! 48: -- = k if u(k,k) = 0.0.
! 49: -- This is not an error for this routine,
! 50: -- but it does indicate that lusolve will
! 51: -- divide by zero if called. Use rcond in
! 52: -- lufco for a reliable indication of singularity.
! 53:
! 54: procedure lufco ( a : in out Matrix; n : in natural;
! 55: ipvt : out Standard_Natural_Vectors.Vector;
! 56: rcond : out number );
! 57:
! 58: -- DESCRIPTION :
! 59: -- lufco factors a floating point matrix by gaussian elimination
! 60: -- and estimates the condition of the matrix.
! 61:
! 62: -- If rcond is not needed, lufac is slightly faster.
! 63: -- To solve a*x = b, follow lufco by lusolve.
! 64:
! 65: -- ON ENTRY :
! 66: -- a a floating point matrix(1..n,1..n) to be factored;
! 67: -- n the dimension of the matrix a.
! 68:
! 69: -- ON RETURN :
! 70: -- a an upper triangular matrix and the multipliers
! 71: -- which are used to obtain it.
! 72: -- The factorization can be written a = l*u, where
! 73: -- l is a product of permutation and unit lower triangular
! 74: -- matrices and u is upper triangular.
! 75: -- ipvt a natural vector of pivot indices
! 76: -- rcond an estimate of the reciprocal condition of a.
! 77: -- For the system a*x = b, relative perturbations
! 78: -- in a and b of size epsilon may cause relative
! 79: -- perturbations in x of size epsilon/rcond.
! 80: -- If rcond is so small that the logical expression
! 81: -- 1.0 + rcond = 1.0
! 82: -- is true, than a may be singular to working precision.
! 83: -- In particular, rcond is zero if exact singularity is
! 84: -- detected or the estimate underflows.
! 85:
! 86: procedure lusolve ( a : in Matrix; n : in natural;
! 87: ipvt : in Standard_Natural_Vectors.Vector;
! 88: b : in out Vectors.Vector );
! 89:
! 90: -- DESCRIPTION :
! 91: -- lusolve solves the system a*x = b using the factors
! 92: -- computed by lufac or lufco
! 93:
! 94: -- ON ENTRY :
! 95: -- a a floating-point matrix(1..n,1..n), the output from
! 96: -- lufac or lufco;
! 97: -- n the dimension of the matrix a;
! 98: -- ipvt the pivot vector from lufac or lufco;
! 99: -- b the right hand side vector.
! 100:
! 101: -- ON RETURN :
! 102: -- b the solution vector x.
! 103:
! 104: procedure Triangulate ( a : in out Matrix; n,m : in integer );
! 105:
! 106: -- DESCRIPTION :
! 107: -- triangulate makes the n*m matrix a triangular using
! 108: -- Gaussian elimination.
! 109:
! 110: -- ON ENTRY :
! 111: -- a a floating-point matrix(1..n,1..m);
! 112: -- n the number of rows of a;
! 113: -- m the number of columns of a.
! 114:
! 115: -- ON RETURN :
! 116: -- a the triangulated matrix.
! 117:
! 118: procedure Diagonalize ( a : in out Matrix; n,m : in integer );
! 119:
! 120: -- DESCRIPTION :
! 121: -- Diagonalize makes the n*m floating-point matrix a diagonal.
! 122:
! 123: -- ON ENTRY :
! 124: -- a a floating-point matrix(1..n,1..m);
! 125: -- n the number of rows of a;
! 126: -- m the number of columns of a.
! 127:
! 128: -- ON RETURN :
! 129: -- a the diagonalized matrix.
! 130:
! 131: -- DYNAMIC TRIANGULATORS :
! 132:
! 133: procedure Upper_Triangulate
! 134: ( row : in natural; mat : in out Matrix; tol : in number;
! 135: ipvt : in out Standard_Natural_Vectors.Vector;
! 136: pivot : out integer );
! 137:
! 138: -- DESCRIPTION :
! 139: -- Makes the matrix upper triangular by updating the current row of
! 140: -- the matrix. If pivot = 0 on return, then the matrix is singular.
! 141:
! 142: -- REQUIRED :
! 143: -- The matrix is upper triangular up to current row, which means that
! 144: -- abs(max(i,i)) > tol, for i in mat'first(1)..row-1.
! 145: -- The matrix might have more columns than rows, but ipvt'range should
! 146: -- match the range of the first columns in the matrix.
! 147:
! 148: procedure Upper_Triangulate
! 149: ( roweli : in natural; elim : in Matrix; tol : in number;
! 150: rowmat : in natural; mat : in out Matrix );
! 151: procedure Upper_Triangulate
! 152: ( roweli : in natural; elim : in Matrix; tol : in number;
! 153: firstrow,lastrow : in natural; mat : in out Matrix );
! 154:
! 155: -- DESCRIPTION :
! 156: -- Using the matrix elim, the unknown roweli is eliminated in mat.
! 157:
! 158: -- ON ENTRY :
! 159: -- roweli current unknown to be eliminated;
! 160: -- elim elim(1..roweli,m'range(2)) is upper triangular;
! 161: -- firstrow indicates start in range of mat to be updated;
! 162: -- lastrow indicates end in range of mat to be updated;
! 163: -- roweli indicates the current unknown to be updated;
! 164: -- mat the unknows before roweli are already eliminated.
! 165:
! 166: -- ON RETURN :
! 167: -- mat the updated matrix after elimination of roweli.
! 168:
! 169: procedure Switch ( ipvt : in Standard_Natural_Vectors.Vector;
! 170: row : in integer; mat : in out Matrix );
! 171:
! 172: -- DESCRIPTION :
! 173: -- Applies the pivoting information ipvt to the row of the matrix.
! 174:
! 175: procedure Switch ( k,pivot,first,last : in integer; mat : in out Matrix );
! 176:
! 177: -- DESCRIPTION :
! 178: -- Swaps the column k with the pivot column for row first..last in mat.
! 179:
! 180: function Solve ( mat : Matrix; tol : number;
! 181: ipvt : Standard_Natural_Vectors.Vector )
! 182: return Vectors.Vector;
! 183:
! 184: -- DESCRIPTION :
! 185: -- Returns a solution to the homogenous linear system with inequalities
! 186: -- in the rows of the matrix.
! 187:
! 188: -- REQUIRED : the matrix is upper triangular.
! 189:
! 190: function Solve ( n,col : natural; mat : Matrix )
! 191: return Vectors.Vector;
! 192:
! 193: -- DESCRIPTION :
! 194: -- Solves the system with as right-hand side the column indicated
! 195: -- by the parameter col.
! 196:
! 197: -- REQUIRED :
! 198: -- The matrix must be upper triangular.
! 199:
! 200: end Generic_Floating_Linear_Solvers;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to generic_floating_linear_solvers.ads 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