source: git/kernel/linearAlgebra.h @ 57592b3

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