source: git/kernel/linearAlgebra.h @ 4f042c2

spielwiese
Last change on this file since 4f042c2 was a704475, checked in by Frank Seelisch <seelisch@…>, 13 years ago
slight changes in in-code docu of factmodd git-svn-id: file:///usr/local/Singular/svn/trunk@13854 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 26.5 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 follows immediately that also A = P * L * U,
42 * since P is self-inverse.<br>
43 *
44 * The method will modify all argument matrices except aMat, so that
45 * they will finally contain the matrices P, L, and U as specified
46 * above.
47 **/
48void luDecomp(
49       const matrix aMat, /**< [in]  the initial matrix A,          */
50       matrix &pMat,      /**< [out] the row permutation matrix P   */
51       matrix &lMat,      /**< [out] the lower triangular matrix L  */
52       matrix &uMat       /**< [out] the upper row echelon matrix U */
53             );
54
55/**
56 * Returns a pivot score for any non-zero matrix entry.
57 *
58 * The smaller the score the better will n serve as a pivot element
59 * in subsequent Gauss elimination steps.
60 *
61 * @return the pivot score of n
62 **/
63int pivotScore(
64               number n  /**< [in] a non-zero matrix entry */
65              );
66
67/**
68 * Finds the best pivot element in some part of a given matrix.
69 *
70 * Given any matrix A with valid row indices r1..r2 and valid column
71 * indices c1..c2, this method finds the best pivot element for
72 * subsequent Gauss elimination steps in A[r1..r2, c1..c2]. "Best"
73 * here means best with respect to the implementation of the method
74 * 'pivotScore(number n)'.<br>
75 * In the case that all elements in A[r1..r2, c1..c2] are zero, the
76 * method returns false, otherwise true.
77 *
78 * @return false if all relevant matrix entries are zero, true otherwise
79 * @sa pivotScore(number n)
80 **/
81bool pivot(
82           const matrix aMat, /**< [in]  any matrix with number entries */
83           const int r1,      /**< [in]  lower row index                */
84           const int r2,      /**< [in]  upper row index                */
85           const int c1,      /**< [in]  lower column index             */
86           const int c2,      /**< [in]  upper column index             */
87           int* bestR,        /**< [out] address of row index of best
88                                         pivot element                  */
89           int* bestC         /**< [out] address of column index of
90                                         best pivot element             */
91          );
92
93/**
94 * Computes the inverse of a given (n x n)-matrix.
95 *
96 * This method expects an (n x n)-matrix, that is, it must have as many
97 * rows as columns. Inversion of the first argument is attempted via the
98 * LU-decomposition. There are two cases:<br>
99 * 1) The matrix is not invertible. Then the method returns false, and
100 *    &iMat remains unchanged.<br>
101 * 2) The matrix is invertible. Then the method returns true, and the
102 *    content of iMat is the inverse of aMat.
103 *
104 * @return true iff aMat is invertible, false otherwise
105 * @sa luInverseFromLUDecomp(const matrix pMat, const matrix lMat,
106 *                           const matrix uMat, matrix &iMat)
107 **/
108bool luInverse(
109               const matrix aMat, /**< [in]  matrix to be inverted */
110               matrix &iMat       /**< [out] inverse of aMat if
111                                             invertible            */
112              );
113
114/**
115 * Computes the inverse of a given (n x n)-matrix in upper right
116 * triangular form.
117 *
118 * This method expects a quadratic matrix, that is, it must have as
119 * many rows as columns. Moreover, uMat[i, j] = 0, at least for all
120 * i > j, that is, u is in upper right triangular form.<br>
121 * If the argument diagonalIsOne is true, then we know additionally,
122 * that uMat[i, i] = 1, for all i. In this case uMat is invertible.
123 * Contrariwise, if diagonalIsOne is false, we do not know anything
124 * about the diagonal entries. (Note that they may still all be
125 * 1.)<br>
126 * In general, there are two cases:<br>
127 * 1) The matrix is not invertible. Then the method returns false,
128 *    and &iMat remains unchanged.<br>
129 * 2) The matrix is invertible. Then the method returns true, and
130 *    the content of iMat is the inverse of uMat.
131 *
132 * @return true iff uMat is invertible, false otherwise
133 **/
134bool upperRightTriangleInverse(
135       const matrix uMat, /**< [in]  (n x n)-matrix in upper right
136                                     triangular form               */
137       matrix &iMat,      /**< [out] inverse of uMat if invertible */
138       bool diagonalIsOne /**< [in]  if true, then all diagonal
139                                     entries of uMat are 1         */
140                              );
141
142/**
143 * Computes the inverse of a given (n x n)-matrix in lower left
144 * triangular form.
145 *
146 * This method expects an (n x n)-matrix, that is, it must have as
147 * many rows as columns. Moreover, lMat[i,j] = 0, at least for all
148 * j > i, that ism lMat is in lower left triangular form.<br>
149 * If the argument diagonalIsOne is true, then we know additionally,
150 * that lMat[i, i] = 1, for all i. In this case lMat is invertible.
151 * Contrariwise, if diagonalIsOne is false, we do not know anything
152 * about the diagonal entries. (Note that they may still all be
153 * 1.)<br>
154 * In general, there are two cases:<br>
155 * 1) The matrix is not invertible. Then the method returns false,
156 *    and &iMat remains unchanged.<br>
157 * 2) The matrix is invertible. Then the method returns true, and
158 *    the content of iMat is the inverse of lMat.
159 *
160 * @return true iff lMat is invertible, false otherwise
161 **/
162bool lowerLeftTriangleInverse(
163       const matrix lMat, /**< [in]  (n x n)-matrix in lower left
164                                     triangular form               */
165       matrix &iMat,      /**< [out] inverse of lMat if invertible */
166       bool diagonalIsOne /**< [in]  if true, then all diagonal
167                                     entries of lMat are 1         */
168                              );
169
170/**
171 * Computes the inverse of an (n x n)-matrix which is given by its LU-
172 * decomposition.
173 *
174 * With A denoting the matrix to be inverted, the method expects the
175 * LU-decomposition of A, that is, pMat * A = lMat * uMat, where
176 * the argument matrices have the appropriate proteries; see method
177 * 'luDecomp(const matrix aMat, matrix &pMat, matrix &lMat,
178 * matrix &uMat)'.<br>
179 * Furthermore, uMat is expected to be an (n x n)-matrix. Then A^(-1)
180 * is computed according to A^(-1) = uMat^(-1) * lMat^(-1) * pMat,
181 * since pMat is self-inverse. This will work if and only if uMat is
182 * invertible, because lMat and pMat are in any case invertible.<br>
183 * In general, there are two cases:<br>
184 * 1) uMat and hence A is not invertible. Then the method returns
185 *    false, and &iMat remains unchanged.<br>
186 * 2) uMat and hence A is invertible. Then the method returns true,
187 *    and the content of iMat is the inverse of A.
188 *
189 * @return true if A is invertible, false otherwise
190 * @sa luInverse(const matrix aMat, matrix &iMat)
191 **/
192bool luInverseFromLUDecomp(
193       const matrix pMat, /**< [in]  permutation matrix of an LU-
194                                     decomposition                */
195       const matrix lMat, /**< [in]  lower left matrix of an LU-
196                                     decomposition                */
197       const matrix uMat, /**< [in]  upper right matrix of an LU-
198                                     decomposition                */
199       matrix &iMat       /**< [out] inverse of A if invertible   */
200                          );
201
202/**
203 * Computes the rank of a given (m x n)-matrix.
204 *
205 * The matrix may already be given in row echelon form, e.g., when
206 * the user has before called luDecomp and passes the upper triangular
207 * matrix U to luRank. In this case, the second argument can be set to
208 * true in order to pass this piece of information to the method.
209 * Otherwise, luDecomp will be called first to compute the matrix U.
210 * The rank is then read off the matrix U.
211 *
212 * @return the rank as an int
213 * @sa luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)
214 **/
215int luRank(
216       const matrix aMat,      /**< [in]  input matrix */
217       const bool isRowEchelon /**< [in]  if true then aMat is assumed to be
218                                          already in row echelon form */
219          );
220
221/**
222 * Solves the linear system A * x = b, where A is an (m x n)-matrix
223 * which is given by its LU-decomposition.
224 *
225 * The method expects the LU-decomposition of A, that is,
226 * pMat * A = lMat * uMat, where the argument matrices have the
227 * appropriate proteries; see method
228 * 'luDecomp(const matrix aMat, matrix &pMat, matrix &lMat,
229 * matrix &uMat)'.<br>
230 * Instead of trying to invert A and return x = A^(-1)*b, this
231 * method
232 * 1) computes b' = pMat * b,
233 * 2) solves the simple system L * y = b', and then
234 * 3) solves the simple system U * x = y.
235 * Note that steps 1) and 2) will always work, as L is in any case
236 * invertible. Moreover, y is uniquely determined. Step 3) will only
237 * work if and only if y is in the column span of U. In that case,
238 * there may be infinitely many solutions.
239 * Thus, there are three cases:<br>
240 * 1) There is no solution. Then the method returns false, and &xVec
241 *    as well as &H remain unchanged.<br>
242 * 2) There is a unique solution. Then the method returns true and
243 *    H is the 1x1 matrix with zero entry.<br>
244 * 3) There are infinitely many solutions. Then the method returns
245 *    true and some solution of the given original linear system.
246 *    Furthermore, the columns of H span the solution space of the
247 *    homogeneous linear system. The dimension of the solution
248 *    space is then the number of columns of H.
249 *
250 * @return true if there is at least one solution, false otherwise
251 **/
252bool luSolveViaLUDecomp(
253       const matrix pMat, /**< [in]  permutation matrix of an LU-
254                                     decomposition                */
255       const matrix lMat, /**< [in]  lower left matrix of an LU-
256                                     decomposition                */
257       const matrix uMat, /**< [in]  upper right matrix of an LU-
258                                     decomposition                */
259       const matrix bVec, /**< [in]  right-hand side of the linear
260                                     system to be solved          */
261       matrix &xVec,      /**< [out] solution of A*x = b          */
262       matrix &H          /**< [out] matrix with columns spanning
263                                     homogeneous solution space   */
264                          );
265
266/**
267 * Solves the linear system A * x = b, where A is an (m x n)-matrix
268 * which is given by its LDU-decomposition.
269 *
270 * The method expects the LDU-decomposition of A, that is,
271 * pMat * A = lMat * dMat^(-1) * uMat, where the argument matrices have
272 * the appropriate proteries; see method
273 * 'lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat,
274 * matrix &dMat, matrix &uMat, poly &l, poly &u, poly &lTimesU)'.<br>
275 * Instead of trying to invert A and return x = A^(-1)*b, this
276 * method
277 * 1) computes b' = l * pMat * b,
278 * 2) solves the simple system L * y = b',
279 * 3) computes y' = u * dMat * y,
280 * 4) solves the simple system U * y'' = y',
281 * 5) computes x = 1/(lTimesU) * y''.
282 * Note that steps 1), 2) and 3) will always work, as L is in any case
283 * invertible. Moreover, y and thus y' are uniquely determined.
284 * Step 4) will only work if and only if y' is in the column span of U.
285 * In that case, there may be infinitely many solutions.
286 * In contrast to luSolveViaLUDecomp, this algorithm guarantees that
287 * all divisions which have to be performed in steps 2) and 4) will
288 * work without remainder. This is due to properties of the given LDU-
289 * decomposition. Only the final step 5) performs a division of a vector
290 * by a member of the ground field. (Suppose the ground field is Q or
291 * some rational function field. Then, if A contains only integers or
292 * polynomials, respectively, then steps 1) - 4) will keep all data
293 * integer or polynomial, respectively. This may speed up computations
294 * considerably.)
295 * For the outcome, there are three cases:<br>
296 * 1) There is no solution. Then the method returns false, and &xVec
297 *    as well as &H remain unchanged.<br>
298 * 2) There is a unique solution. Then the method returns true and
299 *    H is the 1x1 matrix with zero entry.<br>
300 * 3) There are infinitely many solutions. Then the method returns
301 *    true and some solution of the given original linear system.
302 *    Furthermore, the columns of H span the solution space of the
303 *    homogeneous linear system. The dimension of the solution
304 *    space is then the number of columns of H.
305 *
306 * @return true if there is at least one solution, false otherwise
307 **/
308bool luSolveViaLDUDecomp(
309       const matrix pMat,  /**< [in]  permutation matrix of an LDU-
310                                      decomposition                     */
311       const matrix lMat,  /**< [in]  lower left matrix of an LDU-
312                                      decomposition                     */
313       const matrix dMat,  /**< [in]  diagonal matrix of an LDU-
314                                      decomposition                     */
315       const matrix uMat,  /**< [in]  upper right matrix of an LDU-
316                                      decomposition                     */
317       const poly l,       /**< [in]  pivot product l of an LDU decomp. */
318       const poly u,       /**< [in]  pivot product u of an LDU decomp. */
319       const poly lTimesU, /**< [in]  product of l and u                */
320       const matrix bVec,  /**< [in]  right-hand side of the linear
321                                      system to be solved               */
322       matrix &xVec,       /**< [out] solution of A*x = b               */
323       matrix &H           /**< [out] matrix with columns spanning
324                                      homogeneous solution space        */
325                          );
326
327/**
328 * Creates a new matrix which is the (nxn) unit matrix, and returns true
329 * in case of success.
330 *
331 * The method will be successful whenever n >= 1. In this case and only then
332 * true will be returned and the new (nxn) unit matrix will reside inside
333 * the second argument.
334 *
335 * @return true iff the (nxn) unit matrix could be build
336 **/
337bool unitMatrix(
338       const int n,     /**< [in]  size of the matrix */
339       matrix &unitMat  /**< [out] the new (nxn) unit matrix */
340               );
341
342/**
343 * Creates a new matrix which is a submatrix of the first argument, and
344 * returns true in case of success.
345 *
346 * The method will be successful whenever rowIndex1 <= rowIndex2 and
347 * colIndex1 <= colIndex2. In this case and only then true will be
348 * returned and the last argument will afterwards contain a copy of the
349 * respective submatrix of the first argument.
350 * Note that this method may also be used to copy an entire matrix.
351 *
352 * @return true iff the submatrix could be build
353 **/
354bool subMatrix(
355       const matrix aMat,    /**< [in]  the input matrix */
356       const int rowIndex1,  /**< [in]  lower row index */
357       const int rowIndex2,  /**< [in]  higher row index */
358       const int colIndex1,  /**< [in]  lower column index */
359       const int colIndex2,  /**< [in]  higher column index */
360       matrix &subMat        /**< [out] the new submatrix */
361              );
362
363/**
364 * Swaps two rows of a given matrix and thereby modifies it.
365 *
366 * The method expects two valid, distinct row indices, i.e. in
367 * [1..n], where n is the number of rows in aMat.
368 **/
369void swapRows(
370       int row1,     /**< [in]     index of first row to swap */
371       int row2,     /**< [in]     index of second row to swap */
372       matrix& aMat  /**< [in/out] matrix subject to swapping */
373             );
374
375/**
376 * Swaps two columns of a given matrix and thereby modifies it.
377 *
378 * The method expects two valid, distinct column indices, i.e. in
379 * [1..n], where n is the number of columns in aMat.
380 **/
381void swapColumns(
382       int column1,  /**< [in]     index of first column to swap */
383       int column2,  /**< [in]     index of second column to swap */
384       matrix& aMat  /**< [in/out] matrix subject to swapping */
385             );
386
387/**
388 * Creates a new matrix which contains the first argument in the top left
389 * corner, and the second in the bottom right.
390 *
391 * All other entries of the resulting matrix which will be created in the
392 * third argument, are zero.
393 **/
394void matrixBlock(
395       const matrix aMat,    /**< [in]  the top left input matrix */
396       const matrix bMat,    /**< [in]  the bottom right input matrix */
397       matrix &block         /**< [out] the new block matrix */
398                );
399
400/**
401 * Computes the characteristic polynomial from a quadratic (2x2) matrix
402 * and returns true in case of success.
403 *
404 * The method will be successful whenever the input matrix is a (2x2) matrix.
405 * In this case, the resulting polynomial will be a univariate polynomial in
406 * the first ring variable of degree 2 and it will reside in the second
407 * argument.
408 * The method assumes that the matrix entries are all numbers, i.e., elements
409 * from the ground field/ring.
410 *
411 * @return true iff the input matrix was (2x2)
412 **/
413bool charPoly(
414       const matrix aMat,    /**< [in]  the input matrix */
415       poly &charPoly        /**< [out] the characteristic polynomial */
416             );
417
418/**
419 * Computes the square root of a non-negative real number and returns
420 * it as a new number.
421 *
422 * The method assumes that the current ground field is the complex
423 * numbers, and that the argument has imaginary part zero.
424 * If the real part is negative, then false is returned, otherwise true.
425 * The square root will be computed in the last argument by Heron's
426 * iteration formula with the first argument as the starting value. The
427 * iteration will stop as soon as two successive values (in the resulting
428 * sequence) differ by no more than the given tolerance, which is assumed
429 * to be a positive real number.
430 *
431 * @return the square root of the non-negative argument number
432 **/
433bool realSqrt(
434       const number n,          /**< [in]  the input number */
435       const number tolerance,  /**< [in]  accuracy of iteration */
436       number &root             /**< [out] the root of n */
437             );
438
439/**
440 * Computes the Hessenberg form of a given square matrix.
441 *
442 * The method assumes that all matrix entries are numbers coming from some
443 * subfield of the reals..
444 * Afterwards, the following conditions will hold:
445 * 1) hessenbergMat = pMat * aMat * pMat^{-1},
446 * 2) hessenbergMat is in Hessenberg form.
447 * During the algorithm, pMat will be constructed as the product of self-
448 * inverse matrices.
449 * The algorithm relies on computing square roots of floating point numbers.
450 * These will be combuted by Heron's iteration formula, with iteration
451 * stopping when two successive approximations of the square root differ by
452 * no more than the given tolerance, which is assumed to be a positve real
453 * number.
454 **/
455void hessenberg(
456       const matrix aMat,      /**< [in]  the square input matrix */
457       matrix &pMat,           /**< [out] the transformation matrix */
458       matrix &hessenbergMat,  /**< [out] the Hessenberg form of aMat */
459       const number tolerance  /**< [in]  accuracy */
460               );
461
462/**
463 * Returns all solutions of a quadratic polynomial relation with real
464 * coefficients.
465 *
466 * The method assumes that the polynomial is univariate in the first
467 * ring variable, and that the current ground field is the complex numbers.
468 * The polynomial has degree <= 2. Thus, there may be
469 * a) infinitely many zeros, when p == 0; then -1 is returned;
470 * b) no zero, when p == constant <> 0; then 0 is returned;
471 * c) one zero, when p is linear; then 1 is returned;
472 * d) one double zero; then 2 is returned;
473 * e) two distinct zeros; then 3 is returned;
474 * In the case e), s1 and s2 will contain the two distinct solutions.
475 * In cases c) and d) s1 will contain the single/double solution.
476 *
477 * @return the number of distinct zeros
478 **/
479int quadraticSolve(
480       const poly p,           /**< [in]  the polynomial        */
481       number &s1,             /**< [out] first zero, if any    */
482       number &s2,             /**< [out] second zero, if any   */
483       const number tolerance  /**< [in] accuracy               */
484                  );
485
486/**
487 * Computes all eigenvalues of a given real quadratic matrix with
488 * multiplicites.
489 *
490 * The method assumes that the current ground field is the complex numbers.
491 * Computations are based on the QR double shift algorithm involving
492 * Hessenberg form and householder transformations.
493 * If the algorithm works, then it returns a list with two entries which
494 * are again lists of the same size:
495 * _[1][i] is the i-th mutually distinct eigenvalue that was found,
496 * _[2][i] is the (int) multiplicity of _[1][i].
497 * If the algorithm does not work (due to an ill-posed matrix), a list with
498 * the single entry (int)0 is returned.
499 * 'tol1' is used for detection of deflation in the actual qr double shift
500 * algorithm.
501 * 'tol2' is used for ending Heron's iteration whenever square roots
502 * are being computed.
503 * 'tol3' is used to distinguish between distinct eigenvalues: When
504 * the Euclidean distance between two computed eigenvalues is less then
505 * tol3, then they will be regarded equal, resulting in a higher
506 * multiplicity of the corresponding eigenvalue.
507 *
508 * @return a list with one entry (int)0, or two entries which are again lists
509 **/
510lists qrDoubleShift(
511       const matrix A,     /**< [in]  the quadratic matrix         */
512       const number tol1,  /**< [in]  tolerance for deflation      */
513       const number tol2,  /**< [in]  tolerance for square roots   */
514       const number tol3   /**< [in]  tolerance for distinguishing
515                                      eigenvalues                  */
516                   );
517
518/**
519 * Computes a factorization of a polynomial h(x, y) in K[[x]][y] up to a
520 * certain degree in x, whenever a factorization of h(0, y) is given.
521 *
522 * The algorithm is based on Hensel's lemma: Let h(x, y) denote a monic
523 * polynomial in y of degree m + n with coefficients in K[[x]]. Suppose there
524 * are two monic factors f_0(y) (of degree n) and g_0(y) of degree (m) such
525 * that h(0, y) = f_0(y) * g_0(y) and <f_0, g_0> = K[y]. Fix an integer d >= 0.
526 * Then there are monic polynomials in y with coefficients in K[[x]], namely
527 * f(x, y) of degree n and g(x, y) of degree m such that
528 *    h(x, y) = f(x, y) * g(x, y) modulo <x^(d+1)>   (*).
529 *
530 * This implementation expects h, f0, g0, and d as described and computes the
531 * factors f and g. Effectively, h will be given as an element of K[x, y] since
532 * all terms of h with x-degree larger than d can be ignored due to (*).
533 * The method expects the ground ring to contain at least two variables; then
534 * x is the first ring variable, specified by the input index xIndex, and y the
535 * second one, specified by yIndex.
536 *
537 * This code was placed here since the algorithm works by successively solving
538 * d linear equation systems. It is hence an application of other methods
539 * defined in this h-file and its corresponding cc-file.
540 *
541 **/
542void henselFactors(
543       const int xIndex,  /**< [in]  index of x in {1, ..., nvars(basering)} */
544       const int yIndex,  /**< [in]  index of y in {1, ..., nvars(basering)} */
545       const poly h,      /**< [in]  the polynomial h(x, y) as above         */
546       const poly f0,     /**< [in]  the first univariate factor of h(0, y)  */
547       const poly g0,     /**< [in]  the second univariate factor of h(0, y) */
548       const int d,       /**< [in]  the degree bound, d >= 0                */
549       poly &f,           /**< [out] the first factor of h(x, y)             */
550       poly &g            /**< [out] the second factor of h(x, y)            */
551                              );
552
553/**
554 * LU-decomposition of a given (m x n)-matrix with performing only those
555 * divisions that yield zero remainders.
556 *
557 * Given an (m x n) matrix A, the method computes its LDU-decomposition,
558 * that is, it computes matrices P, L, D, and U such that<br>
559 * - P * A = L * D^(-1) * U,<br>
560 * - P is an (m x m) permutation matrix, i.e., its row/columns form the
561 *   standard basis of K^m,<br>
562 * - L is an (m x m) matrix in lower triangular form of full rank,<br>
563 * - D is an (m x m) diagonal matrix with full rank, and<br>
564 * - U is an (m x n) matrix in upper row echelon form.<br>
565 * From these conditions, it follows immediately that also
566 * A = P * L * D^(-1) * U, since P is self-inverse.<br>
567 *
568 * In contrast to luDecomp, this method only performs those divisions which
569 * yield zero remainders. Hence, when e.g. computing over a rational function
570 * field and starting with polynomial entries only (or over Q and starting
571 * with integer entries only), then any newly computed matrix entry will again
572 * be polynomial (or integer).
573 *
574 * The method will modify all argument matrices except aMat, so that
575 * they will finally contain the matrices P, L, D, and U as specified
576 * above. Moreover, in order to further speed up subsequent calls of
577 * luSolveViaLDUDecomp, the two denominators and their product will also be
578 * returned.
579 **/
580void lduDecomp(
581       const matrix aMat, /**< [in]  the initial matrix A,          */
582       matrix &pMat,      /**< [out] the row permutation matrix P   */
583       matrix &lMat,      /**< [out] the lower triangular matrix L  */
584       matrix &dMat,      /**< [out] the diagonal matrix D          */
585       matrix &uMat,      /**< [out] the upper row echelon matrix U */
586       poly &l,           /**< [out] the product of pivots of L     */
587       poly &u,           /**< [out] the product of pivots of U     */
588       poly &lTimesU      /**< [out] the product of l and u         */
589             );
590
591#endif
592/* LINEAR_ALGEBRA_H */
Note: See TracBrowser for help on using the repository browser.