Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/generic_floating_linear_solvers.ads, Revision 1.1.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