# source:git/kernel/linearAlgebra.h@fec53d

fieker-DuValspielwiese
Last change on this file since fec53d was f86c48, checked in by Hans Schoenemann <hannes@…>, 14 years ago
• Property mode set to `100644`
File size: 11.6 KB
Line
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file lineareAlgebra.h
5 *
6 * This file provides basic linear algebra functionality.
7 *
8 * ABSTRACT: This file provides basic algorithms from linear algebra over
9 * any SINGULAR-supported field.
10 * For the time being, the procedures defined in this file expect matrices
11 * containing objects of the SINGULAR type 'number'. This means that, when
12 * 'currentRing' represents some polynomial ring K[x_1, x_2, ..., x_n], then
13 * the entries of the matrices are 'numbers' representing elements of K (and
14 * NOT 'polys' in K[x_1, x_2, ..., x_n]).
15 * This restriction may become obselete in the future.
16 *
17 * @author Frank Seelisch
18 *
19 * @internal @version \\$Id\$
20 *
21 **/
22/*****************************************************************************/
23
24#ifndef LINEAR_ALGEBRA_H
25#define LINEAR_ALGEBRA_H
26
27// include basic SINGULAR structures
28#include <kernel/structs.h>
29
30/**
31 * LU-decomposition of a given (m x n)-matrix.
32 *
33 * Given an (m x n) matrix A, the method computes its LU-decomposition,
34 * that is, it computes matrices P, L, and U such that<br>
35 * - P * A = L * U,<br>
36 * - P is an (m x m) permutation matrix, i.e., its row/columns form the
37 *   standard basis of K^m,<br>
38 * - L is an (m x m) matrix in lower triangular form with all diagonal
39 *   entries equal to 1,<br>
40 * - U is an (m x n) matrix in upper row echelon form.<br>
41 * From these conditions, it easily follows that also A = P * L * U,
42 * since P is self-inverse.<br>
43 * The method will modify all argument matrices except aMat, so that
44 * they will finally contain the matrices P, L, and U as specified
45 * above.
46 **/
47void luDecomp(
48       const matrix aMat, /**< [in]  the initial matrix A,          */
49       matrix &pMat,      /**< [out] the row permutation matrix P   */
50       matrix &lMat,      /**< [out] the lower triangular matrix L  */
51       matrix &uMat       /**< [out] the upper row echelon matrix U */
52             );
53
54/**
55 * Returns a pivot score for any non-zero matrix entry.
56 *
57 * The smaller the score the better will n serve as a pivot element
58 * in subsequent Gauss elimination steps.
59 *
60 * @return the pivot score of n
61 **/
62int pivotScore(
63               number n  /**< [in] a non-zero matrix entry */
64              );
65
66/**
67 * Finds the best pivot element in some part of a given matrix.
68 *
69 * Given any matrix A with valid row indices r1..r2 and valid column
70 * indices c1..c2, this method finds the best pivot element for
71 * subsequent Gauss elimination steps in A[r1..r2, c1..c2]. "Best"
72 * here means best with respect to the implementation of the method
73 * 'pivotScore(number n)'.<br>
74 * In the case that all elements in A[r1..r2, c1..c2] are zero, the
75 * method returns false, otherwise true.
76 *
77 * @return false if all relevant matrix entries are zero, true otherwise
78 * @sa pivotScore(number n)
79 **/
80bool pivot(
81           const matrix aMat, /**< [in]  any matrix with number entries */
82           const int r1,      /**< [in]  lower row index                */
83           const int r2,      /**< [in]  upper row index                */
84           const int c1,      /**< [in]  lower column index             */
85           const int c2,      /**< [in]  upper column index             */
86           int* bestR,        /**< [out] address of row index of best
87                                         pivot element                  */
88           int* bestC         /**< [out] address of column index of
89                                         best pivot element             */
90          );
91
92/**
93 * Computes the inverse of a given (n x n)-matrix.
94 *
95 * This method expects an (n x n)-matrix, that is, it must have as many
96 * rows as columns. Inversion of the first argument is attempted via the
97 * LU-decomposition. There are two cases:<br>
98 * 1) The matrix is not invertible. Then the method returns false, and
99 *    &iMat remains unchanged.<br>
100 * 2) The matrix is invertible. Then the method returns true, and the
101 *    content of iMat is the inverse of aMat.
102 *
103 * @return true iff aMat is invertible, false otherwise
104 * @sa luInverseFromLUDecomp(const matrix pMat, const matrix lMat,
105 *                           const matrix uMat, matrix &iMat)
106 **/
107bool luInverse(
108               const matrix aMat, /**< [in]  matrix to be inverted */
109               matrix &iMat       /**< [out] inverse of aMat if
110                                             invertible            */
111              );
112
113/**
114 * Computes the inverse of a given (n x n)-matrix in upper right
115 * triangular form.
116 *
117 * This method expects a quadratic matrix, that is, it must have as
118 * many rows as columns. Moreover, uMat[i, j] = 0, at least for all
119 * i > j, that is, u is in upper right triangular form.<br>
120 * If the argument diagonalIsOne is true, then we know additionally,
121 * that uMat[i, i] = 1, for all i. In this case uMat is invertible.
122 * Contrariwise, if diagonalIsOne is false, we do not know anything
123 * about the diagonal entries. (Note that they may still all be
124 * 1.)<br>
125 * In general, there are two cases:<br>
126 * 1) The matrix is not invertible. Then the method returns false,
127 *    and &iMat remains unchanged.<br>
128 * 2) The matrix is invertible. Then the method returns true, and
129 *    the content of iMat is the inverse of uMat.
130 *
131 * @return true iff uMat is invertible, false otherwise
132 **/
133bool upperRightTriangleInverse(
134       const matrix uMat, /**< [in]  (n x n)-matrix in upper right
135                                     triangular form               */
136       matrix &iMat,      /**< [out] inverse of uMat if invertible */
137       bool diagonalIsOne /**< [in]  if true, then all diagonal
138                                     entries of uMat are 1         */
139                              );
140
141/**
142 * Computes the inverse of a given (n x n)-matrix in lower left
143 * triangular form.
144 *
145 * This method expects an (n x n)-matrix, that is, it must have as
146 * many rows as columns. Moreover, lMat[i,j] = 0, at least for all
147 * j > i, that ism lMat is in lower left triangular form.<br>
148 * If the argument diagonalIsOne is true, then we know additionally,
149 * that lMat[i, i] = 1, for all i. In this case lMat is invertible.
150 * Contrariwise, if diagonalIsOne is false, we do not know anything
151 * about the diagonal entries. (Note that they may still all be
152 * 1.)<br>
153 * In general, there are two cases:<br>
154 * 1) The matrix is not invertible. Then the method returns false,
155 *    and &iMat remains unchanged.<br>
156 * 2) The matrix is invertible. Then the method returns true, and
157 *    the content of iMat is the inverse of lMat.
158 *
159 * @return true iff lMat is invertible, false otherwise
160 **/
161bool lowerLeftTriangleInverse(
162       const matrix lMat, /**< [in]  (n x n)-matrix in lower left
163                                     triangular form               */
164       matrix &iMat,      /**< [out] inverse of lMat if invertible */
165       bool diagonalIsOne /**< [in]  if true, then all diagonal
166                                     entries of lMat are 1         */
167                              );
168
169/**
170 * Computes the inverse of an (n x n)-matrix which is given by its LU-
171 * decomposition.
172 *
173 * With A denoting the matrix to be inverted, the method expects the
174 * LU-decomposition of A, that is, pMat * A = lMat * uMat, where
175 * the argument matrices have the appropriate proteries; see method
176 * 'luDecomp(const matrix aMat, matrix &pMat, matrix &lMat,
177 * matrix &uMat)'.<br>
178 * Furthermore, uMat is expected to be an (n x n)-matrix. Then A^(-1)
179 * is computed according to A^(-1) = uMat^(-1) * lMat^(-1) * pMat,
180 * since pMat is self-inverse. This will work if and only if uMat is
181 * invertible, because lMat and pMat are in any case invertible.<br>
182 * In general, there are two cases:<br>
183 * 1) uMat and hence A is not invertible. Then the method returns
184 *    false, and &iMat remains unchanged.<br>
185 * 2) uMat and hence A is invertible. Then the method returns true,
186 *    and the content of iMat is the inverse of A.
187 *
188 * @return true if A is invertible, false otherwise
189 * @sa luInverse(const matrix aMat, matrix &iMat)
190 **/
191bool luInverseFromLUDecomp(
192       const matrix pMat, /**< [in]  permutation matrix of an LU-
193                                     decomposition                */
194       const matrix lMat, /**< [in]  lower left matrix of an LU-
195                                     decomposition                */
196       const matrix uMat, /**< [in]  upper right matrix of an LU-
197                                     decomposition                */
198       matrix &iMat       /**< [out] inverse of A if invertible   */
199                          );
200
201/**
202 * Computes the rank of a given (m x n)-matrix.
203 *
204 * The matrix may already be given in row echelon form, e.g., when
205 * the user has before called luDecomp and passes the upper triangular
206 * matrix U to luRank. In this case, the second argument can be set to
207 * true in order to pass this piece of information to the method.
208 * Otherwise, luDecomp will be called first to compute the matrix U.
209 * The rank is then read off the matrix U.
210 *
211 * @return the rank as an int
212 * @sa luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)
213 **/
214int luRank(
215       const matrix aMat,      /**< [in]  input matrix */
216       const bool isRowEchelon /**< [in]  if true then aMat is assumed to be
217                                          already in row echelon form */
218          );
219
220/**
221 * Solves the linear system A*x = b, where A is an (n x n)-matrix
222 * which is given by its LU-decomposition.
223 *
224 * The method expects the LU-decomposition of A, that is,
225 * pMat * A = lMat * uMat, where the argument matrices have the
226 * appropriate proteries; see method
227 * 'luDecomp(const matrix aMat, matrix &pMat, matrix &lMat,
228 * matrix &uMat)'.<br>
229 * Instead of trying to invert A and return x = A^(-1)*b, this
230 * method
231 * 1) computes b' = pMat * b,
232 * 2) solves the simple system L * y = b', and then
233 * 3) solves the simple system U * x = y.
234 * Note that steps 1) and 2) will always work, as L is in any case
235 * invertible. Moreover, y is uniquely determined. Step 3) will only
236 * work if and only if y is in the column span of U. In that case,
237 * there may be infinitely many solutions.
238 * Thus, there are three cases:<br>
239 * 1) There is no solution. Then the method returns false, and &xVec
240 *    as well as &H remain unchanged.<br>
241 * 2) There is a unique solution. Then the method returns true and
242 *    H is the 1x1 matrix with zero entry.<br>
243 * 3) There are infinitely many solutions. Then the method returns
244 *    true and some solution of the given original linear system.
245 *    Furthermore, the columns of H span the solution space of the
246 *    homogeneous linear system. The dimension of the solution
247 *    space is then the number of columns of H.
248 *
249 * @return true if there is at least one solution, false otherwise
250 **/
251bool luSolveViaLUDecomp(
252       const matrix pMat, /**< [in]  permutation matrix of an LU-
253                                     decomposition                */
254       const matrix lMat, /**< [in]  lower left matrix of an LU-
255                                     decomposition                */
256       const matrix uMat, /**< [in]  upper right matrix of an LU-
257                                     decomposition                */
258       const matrix bVec, /**< [in]  right-hand side of the linear
259                                     system to be solved          */
260       matrix &xVec,      /**< [out] solution of A*x = b          */
261       matrix &H          /**< [out] matrix with columns spanning
262                                     homogeneous solution space   */
263                          );
264
265#endif
266/* LINEAR_ALGEBRA_H */
267
Note: See TracBrowser for help on using the repository browser.